# Feasible Adjustments of Statistical Inference in High-Dimensional
Generalized Linear Models

Kazuma Sawaya†, Yoshimasa Uematsu‡, Masaaki Imaizumi†♮ †The University of
Tokyo,
‡Hitotsubashi University,
♮RIKEN Center for Advanced Intelligence Project

###### Abstract.

We developed a statistical inference method applicable to a broad range of
generalized linear models (GLMs) in high-dimensional settings, where the
number of unknown coefficients scales proportionally with the sample size.
Although a pioneering inference method has been developed for logistic
regression, which is a specific instance of GLMs, it is not feasible to apply
this method directly to other GLMs because of unknown hyper-parameters. In
this study, we addressed this limitation by developing a new inference method
designed for a certain class of GLMs. Our method is based on the adjustment of
asymptotic normality in high dimensions and is feasible in the sense that it
is possible even with unknown hyper-parameters. Specifically, we introduce a
novel convex loss-based estimator and its associated system, which are
essential components of inference. Next, we devise a methodology for
identifying the system parameters required by the method. Consequently, we
construct confidence intervals for GLMs in a high-dimensional regime. We prove
that our proposed method has desirable theoretical properties, such as strong
consistency and exact coverage probability. Finally, we experimentally
confirmed its validity.

Date: , Contact<EMAIL_ADDRESS>
##  1\. Introduction

For a pair $(\bm{X},Y)$ of $p$-dimensional random features $\bm{X}$ and random
responses $Y$, we consider the following generalized linear model (GLM):

$\displaystyle\mathbb{E}[Y\mid\bm{X}=\bm{x}]=g\left(\bm{x}^{\top}\bm{\beta}\right),~{}~{}~{}\forall\bm{x}\in\mathbb{R}^{p},$
(1)

where $g:\mathbb{R}\to\mathbb{R}$ is an inverse link function that
monotonically increases, and
$\bm{\beta}=(\beta_{1},\dots,\beta_{p})^{\top}\in\mathbb{R}^{p}$ is an unknown
deterministic coefficient vector. The GLM is one of the most popular classes
of statistical models, including linear regression, logistic regression,
Poisson regression, and many others.

Recently, a high-dimensional regime has been studied. We are particularly
interested in the proportional limit of the coefficient dimension $p\equiv
p(n)$ and sample size $n$:

$\displaystyle
n,p(n)\to\infty\mbox{~{}~{}~{}and~{}~{}~{}}p(n)/n\to\exists\kappa\in(0,\infty).$
(2)

In this setting, the GLM and its special cases have been studied by Bayati and
Montanari [4], Rangan [37], El Karoui et al. [19], Thrampoulidis et al. [45],
Sur and Candès [41], Barbier et al. [3], Bellec [5], to name a few.

Our interest lies in the statistical inference of the unknown coefficient
$\bm{\beta}\equiv\bm{\beta}(n)$ in the regime of (2). The statistical
inference of $\bm{\beta}$ is essential in many applied studies and has been
discussed by many authors such as McCullagh [33], Cordeiro [13], Van de Geer
et al. [47], Tian and Taylor [46], Sur et al. [42], Sur and Candès [41], Fei
and Li [20], Cai et al. [10]. Particularly, a seminal study, Sur and Candès
[41] developed an inferential methodology for logistic regression by adjusting
the asymptotic normality of a maximum likelihood estimator using a state
evolution (SE) system, which is induced by an optimization algorithm for
estimation. This system characterizes some parameters related to the bias
caused by high dimensionality. Zhao et al. [51] extends the approach to
handling features with a more general covariance structure. Salehi et al. [39]
develops a regularized estimator for the logistic regression. Yadlowsky et al.
[50] proposes a useful estimator for a parameter needed for the SE system. We
review these aspects in Section 2.2.

Despite the development of inferential frameworks for logistic regression, it
is not feasible to directly apply SE-based inferences to other GLMs. This is
because of the following reasons: (i) an explicit form of the SE system is not
fully obtained in the general GLM and (ii) the inference requires some unknown
parameters in the SE system, but they are estimable only in the case of
logistic regressions.

In this study, we develop a new statistical method of inference for GLMs with
asymmetric inverse link functions in the high-dimensional regime (2). Our
method is based on the adjustment method for the asymptotic normality proposed
by Sur and Candès [41]. In addition, our method is always feasible because it
is applicable in the presence of unknown parameters in the SE system. Our
approach consists of two steps. First, by employing a new convex loss function
for the estimator of GLMs, we derived a new SE system using a corresponding
algorithm for the loss. Second, we develop novel estimators for the unknown
parameters of the SE system. Consequently, by plugging into the estimator, we
achieve a valid statistical inference.

Our contributions are summarized as follows:

  * •

We develop the statistical inference method for GLMs with an asymmetric
inverse link function in a high-dimensional setting. The method is an
extension of the adjustment method for asymptotic normality only for logistic
regression and is feasible for a wider class of GLMs.

  * •

We derive the new SE system for GLMs and also the method for estimating
unknown parameters. Our estimation method for the SE parameters is new since
it is based on a different idea from those for logistic regression in the
literature.

  * •

We prove the asymptotic validity of our method. Specifically, we show that the
proposed confidence interval asymptotically achieves the exact coverage
probability at a preassigned confidence level. This theoretical result is
supported by our experiments.

###  1.1. Related Works

Statistical estimation and inferential problems in the high-dimensional regime
(2) have been studied. For example, linear models are considered by Bayati and
Montanari [4], El Karoui et al. [19], Thrampoulidis et al. [44, 45], Mei and
Montanari [34], El Karoui [18], Hastie et al. [24], Miolane and Montanari
[35]; logistic regression models are by Sur et al. [42], Sur and Candès [41],
and others [43, 5, 30]. In particular, Barbier et al. [3] presented
theoretical results on GLM estimation errors and learnability in the regime
(2). Its Bayesian model formulation is essentially the same as ours, but the
scope of that paper is different from ours; we suppose that the
hyperparameters are unknown and propose an estimation method for them. We also
extended this to statistical inference, which is feasible for real data
analyses.

For the practical use of inference methods in the regime (2), their
feasibility is also an important concern. To address this issue, many studies
have estimated the unknown hyper-parameters and coordinates necessary for
inference. For the hyper-parameter estimation problem for SE-based inference,
ProbeFrontier by Sur and Candès [41] was the first heuristic approach in the
logistic regression case. This method estimates the hyper-parameters of SE
systems based on the phase transition property of the maximum likelihood
estimators for logistic regression. SLOE by Yadlowsky et al. [50] is a
consistent and useful estimator based on another representation of the SE
system and leave-one-out technique [19]. Bellec [5] provides a more widely
applicable inference framework that allows valid inference even if the link
function is misspecified. For a classification problem, Mai and Liao [31], Mai
et al. [32], Seddik et al. [40] established a feasible method applicable for
inference using prior knowledge of the covariance matrix of the covariate.

Various theoretical tools have been proposed for analyzing statistical models
in the regime (2), including (i) approximate message passing (AMP) algorithms
[16, 8, 4], (ii) convex Gaussian minimax theorem (CGMT) [44, 45], (iii) leave-
one-out (LOO) techniques [19], (iv) second-order Poincaré inequalities [12],
and (v) second-order Stein’s formulae [7, 6]. In addition to (i), Rangan [37]
proposed generalized AMP (GAMP) applicable to GLMs. We use this to
characterize the asymptotic behavior of an estimator of GLMs.

In a high-dimensional regime different from (2), statistical inferences under
sparse conditions have also been studied. Guo and Chen [23], Janková and Van
De Geer [26], Salehi et al. [39], Cai et al. [9], Li et al. [29] considered a
statistical inference for a coefficient vector of high-dimensional GLMs.
Janková et al. [25] proposed a goodness-of-fit testing for the GLM models.

###  1.2. Notation

Define $\mathbb{R}_{+}=(0,\infty)$. For a vector
$\bm{b}=(b_{1},\dots,b_{p})\in\mathbb{R}^{p}$, $b_{j}$ denotes a $j$-th
element of $\bm{b}$ for $j=1,\dots,p$. The cumulative distribution and density
functions of the standard normal distribution a denoted by
$\Phi:\mathbb{R}\to\mathbb{R}$ and $\phi:\mathbb{R}\to\mathbb{R}$,
respectively. For function $f:\mathbb{R}\to\mathbb{R}$, $f^{\prime}(\cdot)$
denotes the derivative of $f(\cdot)$. We say $f(\cdot)$ is $C$-Lipschitz with
$C>0$ if $f(\cdot)$ is Lipschitz continuous and its Lipschitz constant is $C$.
For any convex function $f:\mathbb{R}\to\mathbb{R}$ and constant $\eta>0$, we
define the proximal operator $\mathrm{prox}_{\eta f}:\mathbb{R}\to\mathbb{R}$
as $\mathrm{prox}_{\eta
f}(x)=\operatorname*{arg\,min}_{z\in\mathbb{R}}\left\\{\eta
f(z)+(x-z)^{2}/2\right\\}$.

##  2\. Preliminary

In this section, we discuss the basic concept of GLM and its problem in
moderately high dimensions. We also review the statistical inference framework
for logistic regression based on a state evolution (SE) system.

###  2.1. Generalized Linear Model

Suppose that we observe i.i.d. $n$ pairs $\\{(\bm{X}_{i},Y_{i})\\}_{i=1}^{n}$
of a feature vector $\bm{X}_{i}\equiv\bm{X}_{i}(n)\in\mathbb{R}^{p}$ and a
target variable $Y_{i}\equiv Y_{i}(n)\in\mathcal{Y}$ that follow the GLM (1),
where $\mathcal{Y}$ is a response space, such as
$\mathbb{R},\mathbb{R}_{+},\\{0,1\\},\\{0,1,2,\ldots\\}$, and so on.
Hereafter, we drop the dependence on $n$ whenever it is clear from context. We
assume that the feature vector is generated independently from
$\bm{X}_{i}\sim\mathcal{N}_{p}(\bm{0},\bm{\Sigma})$ using the covariance
matrix $\bm{\Sigma}\equiv\bm{\Sigma}(n)\in\mathbb{R}^{p\times p}$. The GLM can
represent several models by specifying the inverse link function $g(\cdot)$
and distribution of $Y_{i}$ for a given $\bm{X}_{i}$. For example, if
$Y_{i}\mid\bm{X}_{i}\sim\mathrm{Ber}(g(\bm{X}_{i}^{\top}\bm{\beta}))$ with
$g(t)=1/(1+\exp(-t))$, the GLM reduces to a logistic regression model. If
$Y_{i}\mid\bm{X}_{i}\sim\mathrm{Pois}(g(\bm{X}_{i}^{\top}\bm{\beta}))$ with
$g(t)=\exp(t)$, then the GLM represents a Poisson regression model.

If a random variable $Y_{i}\mid\bm{X}_{i}$ has a density function (continuous
case) or a mass function (discrete case) $f(\cdot\mid\bm{X}_{i})$, then the
maximum likelihood estimator (MLE) of $\bm{\beta}$ is defined as

$\displaystyle\widehat{\bm{\beta}}^{\mathrm{MLE}}(n)=\operatorname*{arg\,max}_{b\in\mathbb{R}^{p}}\sum_{i=1}^{n}\log
f(Y_{i}\mid\bm{X}_{i}).$ (3)

In a low-dimensional setting where $n\to\infty$ with fixed $p<\infty$, the MLE
$\widehat{\bm{\beta}}^{\mathrm{MLE}}$ asymptotically obeys the normal
distribution,
$\sqrt{n}(\widehat{\bm{\beta}}^{\mathrm{MLE}}-\bm{\beta})\overset{d}{\to}\mathcal{N}_{p}(\bm{0},\mathcal{I}_{\bm{\beta}}^{-1})$,
where $\mathcal{I}_{\bm{\beta}}$ is the Fisher information matrix at the true
coefficient $\bm{\beta}$. Once the variance is estimated consistently, we can
consider the inference of $\bm{\beta}$ conventionally. By contrast, this
classical approach is no longer valid in a proportionally high-dimensional
setting (2), as presented in Sur and Candès [41].

###  2.2. State Evolution (SE)-based Inference: Logistic Regression Case

For a logistic regression model, Sur and Candès [41], Zhao et al. [51]
developed an inferential framework under the high-dimensional regime (2).
Specifically, the studies showed that a properly corrected MLE is
asymptotically normally distributed. That is, as $n,p(n)\to\infty$ with
$p(n)/n\to\kappa\in(0,1)$, we have the following for any $j=1,\dots,p$:

$\displaystyle\frac{\sqrt{n}(\widehat{\beta}_{j}^{\mathrm{MLE}}-\mu\beta_{j})}{\sigma{/\tau_{j}}}\overset{d}{\to}\mathcal{N}(0,1),$
(4)

where $\tau_{j}^{2}:=\mathrm{Var}(X_{ij}\mid\bm{X}_{i\setminus j})$ is the
conditional variance of the $j$-th element of $X_{i}$ with given
$\bm{X}_{i\setminus j}:=(X_{i,1},\dots,X_{i,j-1},X_{i,j+1},\dots,X_{i,p})$.
Here, $(\mu,\sigma)$ are the state evolution (SE) parameters that satisfy the
following system, called the SE system:

$\displaystyle\begin{cases}\kappa^{2}\sigma^{2}&=\mathbb{E}_{(Q_{1},Q_{2})}\left[2g(Z)\left(\eta
g(D)\right)^{2}\right],\\\ 0&=\mathbb{E}_{(Q_{1},Q_{2})}\left[g(Z)Z(\eta
g(D))\right],\\\ 1-\kappa&=\mathbb{E}_{(Q_{1},Q_{2})}\left[2g(Z)(1+\eta
g^{\prime}(D))^{-1}\right],\end{cases}$ (5)

with $D=\mathrm{prox}_{\eta G}\left(\mu Z+\sqrt{\kappa}\sigma Q_{2}\right)$,
$Z=\gamma Q_{1}$, and $(Q_{1},Q_{2})\sim\mathcal{N}_{2}(\bm{0},\bm{I}_{2})$,
where

$\gamma^{2}=\lim_{n\to\infty}\mathrm{Var}(\bm{X}_{i}(n)^{\top}\bm{\beta}(n))$

is a signal strength (SS) parameter. The covariance structure of $\bm{X}_{i}$
is completely captured by $\gamma^{2}$ in this system. In this logistic
regression case, the inverse link and the integrated inverse link functions
are defined as $g(t)=1/(1+\exp(-t))$ and $G(t):=\log(1+\exp(t))$,
respectively. Here, $\mu,\sigma,\eta\in\mathbb{R}$ are the SE parameters.
Remarkably, SE-based inference reduces the inference problem on high-
dimensional coefficients to an SE system with a small number of scalar SE
parameters.

The SE system is derived as a fixed-point equation of an iterative algorithm
for an optimization problem on estimators. In the logistic regression case, we
derive the SE system (5) from the optimization problem (3). The derivation
process is based on the framework of approximate message passing, the details
of which are provided in Feng et al. [21].

Based on the asymptotic normality of the MLE in (4), we may construct a
confidence interval for each $\beta_{j}$ for $j=1,\ldots,p$, provided that the
SE parameters are available. Specifically, given $(\mu,\sigma)$, the MLE-
centric confidence interval with confidence level $1-\alpha\in(0,1)$ is formed
as

$\displaystyle\mathrm{CI}_{(1-\alpha)}:=\left[\frac{\widehat{\beta}^{\mathrm{MLE}}_{j}}{\mu}-z_{(1-\alpha/2)}\frac{\sigma{/\tau_{j}}}{\sqrt{n}\mu},\frac{\widehat{\beta}^{\mathrm{MLE}}_{j}}{\mu}+z_{(1-\alpha/2)}\frac{\sigma{/\tau_{j}}}{\sqrt{n}\mu}\right],$
(6)

where $z_{\alpha/2}$ is the $(\alpha/2)$-quantile of a standard Gaussian
distribution. By properly specifying of the parameters, we obtain
$\mathrm{Pr}(\beta_{j}\in\mathrm{CI}_{(1-\alpha)})\to 1-\alpha$ as
$n,p(n)\to\infty$ in the sense of (2).

Importantly, we must know the SS parameter $\gamma^{2}$ in advance because the
SE parameters $(\mu,\sigma)$ are understood as a function of $\gamma^{2}$. Sur
and Candès [41] proposed an estimation procedure of $\gamma^{2}$ based on the
phase transitions of their existence of the MLE [11]. Yadlowsky et al. [50]
developed another approach for this estimation problem.

##  3\. Estimator and SE System for GLMs

In this section, we derive a generic estimator for the GLM and associated SE
system. It is necessary to perform SE-based inference as in Sur and Candès
[41].

###  3.1. Surrogate Estimator for Coefficient Vector $\bm{\beta}$

First, we develop an estimator for the true coefficient vector $\bm{\beta}$ by
using a surrogate loss function associated with Auer et al. [2], Agarwal et
al. [1]. We introduce the loss function
$\ell:\mathbb{R}^{p}\times\mathbb{R}^{p}\times\mathcal{Y}\to\mathbb{R}$
defined as

$\displaystyle\ell(\bm{b};\bm{x},y):=G\left(\bm{x}^{\top}\bm{b}\right)-y\bm{x}^{\top}\bm{b},$
(7)

where $G(\cdot)$ is the function satisfying $G^{\prime}=g$. Using the loss, we
define the surrogate estimator

$\displaystyle\widehat{\bm{\beta}}(n):=\operatorname*{arg\,min}_{\bm{b}\in\mathbb{R}^{p}}\sum_{i=1}^{n}\ell(\bm{b};\bm{X}_{i},Y_{i}).$
(8)

The use of a surrogate estimator is justified in two ways. First, if
$g(\cdot)$ monotonically increases, the surrogate loss is convex in $\bm{b}$.
This property provides computational and statistical advantages for a broader
class of GLMs and plays a crucial role in the derivation of SE systems; see
Remark 3.1 below. Second, the minimizer of surrogate risk coincides with the
true coefficient $\bm{\beta}$, which is expressed as follows:

######  Proposition 3.1 (Lemma 1 in Agarwal et al. [1]).

Consider an $\mathbb{R}^{p}\times\mathcal{Y}$-valued random element
$(\bm{X},Y)$ that follows GLM (1). If $g(\cdot)$ is increasing, then the
coefficient vector $\bm{\beta}$ in (1) satisfies

$\displaystyle\bm{\beta}=\operatorname*{arg\,min}_{\bm{b}\in\mathbb{R}^{p}}\mathbb{E}\left[\ell(\bm{b};\bm{X},Y)\mid\bm{X}\right].$
(9)

These results demonstrate that the surrogate loss estimator is a reasonable
extension of MLE for GLMs.

######  Remark 3.1 (Convexity of surrogate loss).

The convex property of the surrogate loss plays a critical role in SE-based
inference. Specifically, SE systems such as (5) can only be constructed if the
estimator is a minimizer of a convex function; for details, see Feng et al.
[21]. Therefore, the squared loss $(y-g(\bm{x}^{\top}\bm{b}))^{2}$ is often
inapplicable.

######  Remark 3.2 (Connection to exponential family).

If we model the distribution of $Y_{i}\mid\bm{X}_{i}$ to follow an exponential
dispersion family [27] and $g(\cdot)$ is chosen as the inverse of the
canonical link function, then the surrogate estimator (8) is equivalent to the
MLE (3). This covers a logistic and a Poisson regression model.

###  3.2. SE System for GLMs

We developed an SE system associated with the surrogate estimator (8). As this
estimator is given by the minimization problem of the convex loss function, we
can derive the SE system for the estimator, as described in Section 2.2.
Recall that $G(\cdot)$ is a function satisfying $G^{\prime}=g$, and
$\gamma^{2}=\lim_{n\to\infty}\mathrm{Var}(\bm{X}_{i}(n)^{\top}\bm{\beta}(n))$
is the signal strength (SS) parameter.

######  Proposition 3.2.

We consider the surrogate estimator in (8) to minimize the surrogate loss
function. If $\kappa\in(0,1)$ and the estimator is bounded, then the
corresponding SE system is given by

$\displaystyle\begin{cases}\kappa^{2}\sigma^{2}&=\eta^{2}\mathbb{E}_{(Q_{1},Q_{2},\bar{Y})}\left[\left(\bar{Y}-g\left(D\right)\right)^{2}\right],\\\
0&=\mathbb{E}_{(Q_{1},Q_{2},\bar{Y})}\left[Z\left(\bar{Y}-g\left(D\right)\right)\right],\\\
1-\kappa&=\mathbb{E}_{(Q_{1},Q_{2},\bar{Y})}\left[(1+\eta
g^{\prime}\left(D\right))^{-1}\right],\end{cases}$ (10)

where $D=\mathrm{prox}_{\eta G}(\mu Z+\sqrt{\kappa}\sigma Q_{2}+\eta\bar{Y})$,
$Z=\gamma Q_{1}$, $(Q_{1},Q_{2})\sim\mathcal{N}_{2}(\bm{0},\bm{I}_{2})$, and
$\bar{Y}\sim P_{\bar{Y}}$ with $P_{\bar{Y}}$ being identical to the
distribution of $\lim_{n\to\infty}Y_{1}(n)$ whose mean is $g(Z)$ where
$Z\overset{\rm d}{=}\lim_{n\to\infty}\bm{X}_{1}(n)^{\top}\bm{\beta}(n)$.

We demonstrate the solutions to the system for Poisson regression in Figure 1
and piecewise regression in Figure 2.

In view of (10), the SE parameters $(\mu,\sigma^{2},\eta)$ depend on three
components: (i) $\kappa=\lim_{n\to\infty}p(n)/n$, (ii) the SS parameter
$\gamma^{2}=\lim_{n\to\infty}\mathrm{Var}(\bm{X}_{i}(n)^{\top}\bm{\beta}(n))$,
and (iii) the distribution of $\bar{Y}$ determined by the model. In practice,
we can set $\kappa=p(n)/n$ and consider it as known. However, the estimation
of $\gamma^{2}$ and other parameters (if present) of the distribution of
$\bar{Y}$ for each GLM is highly nontrivial. In the next section, we propose a
simple method for estimating them.

Figure 1. Computed bias $\mu$ (left) and variance $\sigma^{2}$ (right) for
Poisson regression model
$Y\mid\bm{X}\sim\mathrm{Pois}(\exp(\bm{X}^{\top}\bm{\beta}))$ as the solutions
to the SE system over $100$ simulations with their mean and $95\%$ bootstrap
confidence intervals.

Figure 2. Computed bias $\mu$ (left) and variance $\sigma^{2}$ (right) for
piecewise regression model
$Y\mid\bm{X}\sim\mathcal{N}(g(\bm{X}^{\top}\bm{\beta}),0.2)$ with
$g(t)=\min(5t,0.1t)$ as the solutions to the SE system over $100$ simulations
with their mean and $95\%$ bootstrap confidence intervals.

######  Remark 3.3 (Specification of $\bar{Y}$).

If the distribution of $\bar{Y}$ has a single parameter (e.g., Bernoulli,
Poisson, and exponential distribution), then the model assumption
$\mathbb{E}[Y\mid\bm{X}]=g(\bm{X}^{\top}\bm{\beta})$ implies that $\gamma^{2}$
fully characterizes the inferential parameters $(\mu,\sigma^{2},\eta)$. Thus,
in this case, the SS parameter $\gamma^{2}$ is the only unknown parameter to
be estimated. Meanwhile, estimating the distribution of $\bar{Y}$ is difficult
if there are multiple parameters. Even in such cases, we propose a method for
their estimation, focusing on the case of $\bar{Y}$ following a normal
distribution; see Section 4.2.2.

##  4\. Estimation for Parameters in SE System

In this section, we propose a novel estimator for the SS parameter
$\gamma^{2}$ and other parameters of the distribution of $\bar{Y}$ to access
the SE parameters $(\mu,\sigma^{2})$ for inference. Once we estimate
$\gamma^{2}$ and the parameters of $\bar{Y}$, we can solve system (10) by
substituting the estimators to obtain the SE parameters.

###  4.1. Additive Form of GLM

As the first step, we converted the GLM (1) into another representation. That
is, we convert it into the following additive model with observation
$(Y_{i},\bm{X}_{i})$ as

$\displaystyle Y_{i}=g(\bm{X}_{i}^{\top}\bm{\beta})+e_{i}\quad{\rm with}\quad
e_{i}=Y_{i}-\mathbb{E}[Y_{i}\mid\bm{X}_{i}],~{}i=1,...,n,$ (11)

where $e_{i}$ is a noise variable that satisfies
$\mathbb{E}[e_{i}\mid\bm{X}_{i}]=0$. Note that $e_{i}\mid\bm{X}_{i}$ does not
necessarily follow an identical distribution for $i=1,\dots,n$.

Additive model (11) provides a more concrete distribution structure of
$Y_{i}$. While the original GLM (1) highlights the conditional mean, the
additive model provides a more analyzable form of the distribution by
introducing the variable $e_{i}$. This structure is essential for the
estimation method described below.

###  4.2. Estimation of Parameters

We provide an estimation procedure for $\gamma^{2}$ and moments of $\bar{Y}$.

####  4.2.1. Estimator of SS Parameter $\gamma^{2}$

We consider the estimation of the SS parameter $\gamma^{2}$. By the additive
model of GLMs (11), the strong law of large numbers implies the following as
$n\to\infty$:

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}Y_{i}(n)=\frac{1}{n}\sum_{i=1}^{n}\left\\{g(\bm{X}_{i}(n)^{\top}\bm{\beta}(n))+e_{i}\right\\}\overset{\rm
a.s.}{\longrightarrow}\mathbb{E}\left[g(Z)+e_{1}\right]=\mathbb{E}[g(Z)].$

The last equality follows from the model’s assumption that
$\mathbb{E}[e_{1}]=\mathbb{E}[\mathbb{E}[e_{1}\mid\bm{X}_{1}]]=0$. Given that
$g(\cdot)$ is known, we can simulate $\mathbb{E}[g({Z_{\varsigma}})]$ using
${Z_{\varsigma}}\sim\mathcal{N}(0,\varsigma^{2})$ for each $\varsigma>0$.
Thus, we obtain an estimator of $\gamma$ with the given observation
$\\{Y_{i}\\}_{i=1}^{n}$ as a solution to the linear equation:

$\displaystyle\widehat{\gamma}^{2}:=\left\\{\varsigma^{2}:\frac{1}{n}\sum_{i=1}^{n}Y_{i}=\widehat{\mathbb{E}}_{\varsigma}[g({Z_{\varsigma}})]\right\\},$
(12)

where we denote as $\widehat{\mathbb{E}}_{\varsigma}[\cdot]$ a simulated mean
by generated Gaussian random variables; more precisely, we independently
generate
$Z_{\varsigma}^{(1)},\dots,Z_{\varsigma}^{(m)}\sim\mathcal{N}(0,\varsigma^{2})$
with a large $m\in\mathbb{N}$ and calculate
$\widehat{\mathbb{E}}_{\varsigma}[g(Z_{\varsigma})]:=m^{-1}\sum_{k=1}^{m}g(Z_{\varsigma}^{(k)})$.

Figure 3. Illustration of the estimation procedure of $\widehat{\gamma}$.
Given $n^{-1}\sum_{i=1}^{n}Y_{i}$ and $g(\cdot)$, we seek the intersection of
the simulated curve $f(\varsigma)=\widehat{\mathbb{E}}[g(Z_{\varsigma})]$ (red
curve) and $\tilde{f}(\varsigma)=n^{-1}\sum_{i=1}^{n}Y_{i}$ (blue dotted
line). Then, the horizontal axis of the intersection is the estimator defined
in (12).

This estimation is infeasible if $g(\cdot)$ is an odd function. This is
because $\mathbb{E}[g(Z_{\varsigma})]$ is zero for all $\varsigma>0$, because
the Gaussian density is even, violating the uniqueness of a solution to the
minimization problem in (12).

####  4.2.2. Estimation of Additional Parameter in $\bar{Y}$

For some types of GLMs, the distribution of $\bar{Y}$ in the SE system has an
additional unknown parameter that should be estimated. This estimator is
required unless the distribution of $\bar{Y}\mid Z$ is completely
characterized by its mean $g(Z)$; for example, when $\bar{Y}\mid
Z\sim\mathcal{N}(g(Z),\sigma_{{e}}^{2})$, the variance parameter
$\sigma_{{e}}^{2}>0$ is an additional unknown parameter. By contrast, if we
consider a binary choice model $\bar{Y}\mid Z\sim\mathrm{Ber}(g(Z))$ with
$Z\sim\mathcal{N}(0,\gamma^{2})$, including a logistic regression model, the
distribution of $\bar{Y}$ is fully determined by $\gamma^{2}$. Hence the SE
parameters in the system (10) become available once $\widehat{\gamma}^{2}$ in
(12) is obtained.

We consider the case in which the $e_{i}$ in (11) independently follows
$\mathcal{N}(0,\sigma_{e}^{2})$ with unknown $\sigma_{e}>0$. From the additive
model representation in (11), we have, as $n\to\infty$,

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}Y_{i}(n)^{2}=\frac{1}{n}\sum_{i=1}^{n}\left\\{g(\bm{X}_{i}(n)^{\top}\bm{\beta}(n))+e_{i}\right\\}^{2}\overset{\rm
a.s.}{\longrightarrow}\mathbb{E}[g(Z)^{2}]+\sigma_{{e}}^{2},$

where the last expression follows from the independence of $e_{i}$ and
$\bm{X}_{i}(n)$. Given that $n^{-1}\sum_{i=1}^{n}Y_{i}^{2}$ is provided and
$\mathbb{E}[g(Z)^{2}]$ can be estimated by plugging into the estimator
$\widehat{\gamma}^{2}$, we estimate $\sigma_{{e}}^{2}$ by

$\displaystyle\widehat{\sigma}_{{e}}^{2}:=\left(\frac{1}{n}\sum_{i\in
I_{n}}Y_{i}(n)^{2}\right)-\widehat{\mathbb{E}}_{\widehat{\gamma}}\left[g\left(Z_{\widehat{\gamma}_{I^{c}}}\right)^{2}\right],$
(13)

where $I_{n}\subset\\{1,\dots,n\\}$ such that $\lim_{n\to\infty}|I_{n}|/n$ is
a constant (e.g., $I_{n}=\\{1,2,\ldots,\lfloor n/2\rfloor\\}$),
$I_{n}^{c}=\\{1,\dots,n\\}\backslash I_{n}$, and
$\widehat{\gamma}_{I_{n}^{c}}$ denotes the estimator (12) made only from data
in $I_{n}^{c}$. The data split prevents bias caused by data reuse.

##  5\. Statistical Inference for GLMs

Along with the ingredients obtained thus far, we constructed a confidence
interval with reasonable coverage, even for high dimensions. We also provide
theoretical guarantees for estimators and confidence intervals.

###  5.1. Confidence Interval

We construct a confidence interval for $\bm{\beta}$ by using the estimators in
Sections 4 and $\widehat{\bm{\beta}}\equiv\widehat{\bm{\beta}}(n)$ in (8). Let
$(\widehat{\mu},\widehat{\sigma}^{2},\widehat{\eta})$ be the solutions to the
SE system (10), with $\gamma^{2}$ replaced by $\widehat{\gamma}^{2}$. We also
introduce an estimator $\widehat{\tau}_{j}$ of the conditional variance
$\tau_{j}^{2}$, the details of which are provided in Section A.1. Then, our
confidence interval with a confidence level $1-\alpha\in(0,1)$ is defined as:

$\displaystyle\mathcal{CI}_{{1-\alpha},j}:=\left[\frac{\widehat{\beta}_{j}}{\widehat{\mu}}-z_{(1-\alpha/2)}\frac{\widehat{\sigma}}{\sqrt{n}\widehat{\mu}\widehat{\tau}_{j}},\frac{\widehat{\beta}_{j}}{\widehat{\mu}}+z_{(1-\alpha/2)}\frac{\widehat{\sigma}}{\sqrt{n}\widehat{\mu}\widehat{\tau}_{j}}\right],~{}j=1,\dots,p.$
(14)

###  5.2. Theoretical Validation

We derive the theoretical results of the proposed methods. The following
assumptions were made:

######  Assumption 1.

We consider the following conditions:

  * (A1)

$\bm{X}(n)$ is generated as
$\bm{X}(n)\sim\mathcal{N}_{p}(\bm{0},\bm{\Sigma}(n))$, and
$\mathrm{Var}(\bm{X}(n)^{\top}\bm{\beta}(n))$ has a convergence limit
$\gamma^{2}<\infty$ as $n\to\infty$.

  * (A2)

An inverse link function $g:\mathbb{R}\to\mathbb{R}$ is monotonically
increasing and $L$-smooth (i.e., the derivative of $g(\cdot)$ is $L$-Lipschitz
continuous).

  * (A3)

We can write a GLM as $Y=h(\bm{X}^{\top}\bm{\beta},\varepsilon)$ with a non-
random function $h:\mathbb{R}^{2}\to\mathbb{R}$ and an $\mathbb{R}$-valued
random variable $\varepsilon$ independent of $\bm{X}$, such that $\varepsilon$
has a finite second moment and $h(\cdot,\cdot)$ is Lipschitz continuous with
respect to the first argument.

  * (A4)

There exists a positive solution to SE system (10).

The convergence of
$\mathrm{Var}(\bm{X}^{\top}\bm{\beta})=\bm{\beta}^{\top}\bm{\Sigma}\bm{\beta}$
in (A1) is a constraint on $\bm{\beta}$ and $\bm{\Sigma}$. (A3) is a technical
requirement for the Lipschitzness of SE parameters with respect to
$\gamma^{2}$. Examples of $h(\cdot,\cdot)$ for each model are provided in
Section A.5.

Depending on the model chosen, the MLE can diverge in the current high-
dimensional setting. For instance, the MLE for logistic regression models can
diverge to infinity if $\kappa>0.5$ [11]. However, it is difficult to impose
restrictions that guarantee uniform existence on all GLMs because such
conditions generally vary depending on the model. Therefore, we focus only on
the situation in which the estimator (8) exists asymptotically almost surely;
that is, $\lim_{n\to\infty}\|\widehat{\bm{\beta}}(n)\|<\infty$ a.s. If the
estimator does not exist, we may consider a regularization approach to
guarantee its existence; see Section A.4 for this point.

First, we obtain the asymptotic normality of the adjusted test statistics for
each coordinate of $\widehat{\bm{\beta}}$ with the oracle SE parameters $\mu$
and $\sigma^{2}$.

######  Proposition 5.1.

Suppose that we know a solution to the system of nonlinear equations (10), and
the estimator (8) almost surely exists asymptotically. Under (A1), (A2), and
(A4) in Assumption 1, as $n,p(n)\to\infty$ and $p(n)/n\to\kappa\in(0,1)$, we
have the following for $j=1,...,p(n)$:

$\displaystyle\frac{\sqrt{n}(\widehat{\beta}_{j}-\mu\beta_{j})}{\sigma/\tau_{j}}\overset{{d}}{\to}\mathcal{N}(0,1).$
(15)

This is the first result to demonstrate the marginal asymptotic normality of
estimators for GLMs in a high-dimensional scheme, including MLEs with
canonical links. This differs from the classical result
$\sqrt{n}(\widehat{\bm{\beta}}^{\rm
MLE}-\bm{\beta})\overset{{d}}{\to}\mathcal{N}(\bm{0},\mathcal{I}_{\bm{\beta}}^{-1})$
in a setup with $p<\infty$. We regard the marginal convergence result as an
extension of that of Zhao et al. [51] on MLEs for logistic regression to more
general models. Although this is an important result alone, a feasible
construction of the asymptotic normality is provided in Proposition 5.3.

####  5.2.1. Consistency of Proposed Estimator for SE Parameters

We show the strong consistency of the estimators $\widehat{\gamma}^{2}$ and
$\widehat{\sigma}_{e}^{2}$, as one of our main theoretical contributions. To
identify these estimators, we impose an additional assumption on $g(\cdot)$.

######  Assumption 2.

The odd part $g_{0}(x):=(g(x)+g(-x))/2$ of the inverse link function
$g(\cdot)$ is strictly monotonic for $\mathbb{R}_{+}$.

Intuitively, this implies that the form of $g(\cdot)$ is not point-symmetric
around point $(0,g(0))$, which is a generalization of a non-odd function. Even
if $g(\cdot)$ does not satisfy Assumption 2, we can define an approximator
that is sufficiently close to $g(\cdot)$ and that satisfies Assumption 2. We
will describe its details in Section A.3.

Assumption 2 implies that $\mathbb{E}[g(Z_{\varsigma})]$ is a strictly
monotone function of $\varsigma>0$.

######  Lemma 5.1.

Let $\mathbb{E}[g(Z_{\varsigma})]<\infty$ for $\varsigma>0$ and $g_{0}$ be
strictly monotone on $\mathbb{R}_{+}$. Then, the map
$\varsigma\mapsto\mathbb{E}[g(Z_{\varsigma})]$ is strictly monotonic in
$\varsigma>0$.

This ensures the uniqueness of the estimator $\widehat{\gamma}^{2}$. We then
obtain the following results for the strong consistency of the estimators for
the SE parameters:

######  Theorem 5.1.

Under Assumption 1 (A1) and Assumption 2, $\widehat{\gamma}^{2}$ defined in
(12) satisfies, as $n,p(n)\to\infty$ and $p(n)/n\to\kappa\in(0,\infty)$,

$\displaystyle\widehat{\gamma}^{2}\overset{\mathrm{a.s.}}{\longrightarrow}\gamma^{2}.$

In addition, we discuss the specific case in which $e_{i}$ in (11) is an
independent Gaussian, which requires an additional analysis of its variance
$\sigma_{e}^{2}$. By using the consistency of the signal strength estimator
$\widehat{\gamma}^{2}$, we can obtain the asymptotic validity of the variance
estimator $\widehat{\sigma}_{e}^{2}$.

######  Proposition 5.2.

Suppose that the random variable $e_{i}$ in model (11) follows
$e_{i}\sim\mathcal{N}(0,\sigma_{e}^{2})$ for $i=1,\ldots,n$ independently of
$\bm{X}_{i}$. Then, under (A1) in Assumption 1 and Assumption 2,
$\widehat{\sigma}_{{e}}^{2}$ defined in (13) satisfies the following as
$n,p(n)\to\infty$ and $p(n)/n\to\kappa\in(0,\infty)$:

$\displaystyle\widehat{\sigma}_{{e}}^{2}\overset{\mathrm{a.s.}}{\longrightarrow}\sigma_{e}^{2}.$

####  5.2.2. Asymptotic Validity of Proposed Statistical Inference

As the next theoretical contribution, we demonstrate the asymptotic validity
of our feasible adjustment for statistical inference. To observe this, we used
the asymptotic normality with the estimated SE parameters.

######  Proposition 5.3.

Suppose that all the conditions in Assumption 1 and Assumption 2 hold, the
estimator (8) almost surely exists asymptotically, and
$\widehat{\tau}_{j}^{2}$ is a consistent estimator of the conditional variance
$\tau_{j}^{2}$. Then, for any confidence level $(1-\alpha)\in(0,1)$, we obtain
the following for every $j=1,\dots,p(n)$ satisfying
$\sqrt{n}\tau_{j}\beta_{j}=O(1)$:

$\displaystyle\frac{\sqrt{n}(\widehat{\beta}_{j}-\widehat{\mu}\beta_{j})}{\widehat{\sigma}/\widehat{\tau}_{j}}\overset{d}{\to}\mathcal{N}(0,1),$

as $n,p(n)\to\infty$, where $p(n)/n\to\kappa\in(0,1)$.

Consequently, the proposed confidence interval asymptotically achieved a
preassigned confidence level $(1-\alpha)$.

######  Theorem 5.2.

Assume the settings of Proposition 5.3. Then, for any confidence level
$(1-\alpha)\in(0,1)$, we obtain the following for every $j=1,\dots,p(n)$
satisfying $\sqrt{n}\tau_{j}\beta_{j}=O(1)$:

$\displaystyle\mathrm{Pr}\left(\beta_{j}\in\mathcal{CI}_{{1-\alpha},j}\right)\to
1-\alpha,$

as $n,p(n)\to\infty$, where $p(n)/n\to\kappa\in(0,1)$.

##  6\. Experiment

###  6.1. Empirical Performance of $\widehat{\gamma}$ and
$\widehat{\sigma}_{e}^{2}$

We evaluated the performance of $\widehat{\gamma}$ in (12) using the simulated
data. We fix $n=4000$ and vary $p=\kappa n$ over
$\kappa\in\\{0.1,0.2,0.3,0.4,0.5\\}$. The SS parameter is also fixed at
$\gamma^{2}=1$ and generates $n$ realizations of the feature vector
$\bm{X}\in\mathbb{R}^{p}$ from the $p$-variate normal distribution
$\mathcal{N}_{p}(0,\bm{\Sigma})$ with $\bm{\Sigma}_{ij}=0.5^{|i-j|}$. The true
regression coefficients $\beta_{j}$ for $j=1,\ldots,p$ are drawn independently
from a normal distribution with the variance determined by $\gamma^{2}$. We
calculate $g(\bm{X}^{\top}\bm{\beta})$ conditional on $\bm{X}$ and
$\bm{\beta}$ and finally draw the response $Y$ from a given distribution. We
consider two cases in this section: (i) Poisson regression case
$g(t)=\mathrm{exp}(t)$,
$Y\mid\bm{X}\sim\mathrm{Pois}(g(\bm{X}^{\top}\bm{\beta}))$, (ii) piecewise
regression case $g(t)=\min(5t,0.1t)$,
$Y\mid\bm{X}\sim\mathcal{N}(g(\bm{X}^{\top}\bm{\beta}),0.04)$, and
complementary log-log (cloglog) regression case $g(t)=1-\exp(-\exp(t))$,
$Y\mid\bm{X}\sim\mathrm{Bern}(g(\bm{X}^{\top}\bm{\beta}))$. The results in
Figure 4 demonstrate that the estimator $\widehat{\gamma}$ is consistent for
any $\kappa$.

Figure 4. Numerical performance of $\widehat{\gamma}$ in (12) to estimate
$\gamma=1$ over 1000 replications.

Using the same data-generating process for $\bm{X}$ and $\bm{\beta}$ as
described in the previous section, we study the performance of
$\widehat{\sigma}_{e}^{2}$ (13). Here, we consider two situations: (i)
piecewise regression case $g(t)=\min(5t,0.1t)$,
$Y\mid\bm{X}\sim\mathcal{N}(g(\bm{X}^{\top}\bm{\beta}),0.04)$ and (ii) squared
regression case $g(t)=t^{2}$,
$Y\mid\bm{X}\sim\mathcal{N}(g(\bm{X}^{\top}\bm{\beta}),0.04)$. The results in
Figure 5 demonstrate numerical consistency.

Figure 5. Numerical performance of $\widehat{\sigma}_{e}^{2}$ over 1000
simulations. The blue dotted line represents the true value
$\sigma^{2}_{\varepsilon}=0.04$.

###  6.2. Coverage Proportion of Proposed Confidence Interval

We compared the empirical coverage proportion of the proposed confidence
interval (14) with that of the classical confidence interval with $n=1000$,
$\kappa\in\\{0.1,0.2,0.3,0.4,0.5\\}$, and $\gamma\in\\{1,2\\}$. The proposed
CI is calculated by solving a system of nonlinear equations using the
estimated signal strength. The classical CI for $\beta_{j},j=1,\ldots,p$ is
constructed as

$\displaystyle\left[\widehat{\beta}_{j}-z_{(1-\alpha/2)}({\widehat{\mathcal{I}}_{jj}^{-1}(\widehat{\bm{\beta}})/n})^{1/2},\widehat{\beta}_{j}+z_{(1-\alpha/2)}({\widehat{\mathcal{I}}_{jj}^{-1}(\widehat{\bm{\beta}})/n})^{1/2}\right]$

where
$\widehat{\mathcal{I}}(\widehat{\bm{\beta}})=n^{-1}\sum_{i=1}^{n}g^{\prime}(\bm{X}_{i}^{\top}\widehat{\bm{\beta}})\bm{X}_{i}\bm{X}_{i}^{\top}$
denotes the empirical Fisher information matrix evaluated at
$\widehat{\bm{\beta}}$. We analyze the Poisson regression model, as described
in the previous sections. The results in Figure 6 demonstrate that our
proposed interval achieves theoretical coverage in all cases, whereas the
conventional method decreases the coverage as $\kappa$ increases.

Figure 6. Coverage proportion of true coefficient with $90\%$ classical
(white) and proposed (gray) CIs over $100$ simulations with $n=1000$.

###  6.3. Real Data Application

We consider an application for the Cleveland Clinic Heart Disease dataset [14]
from the UCI Machine Learning Repository [17]. This included $303$
observations of $14$ variables. The target is the presence of heart disease,
which is integer-valued from 0 (no presence) to 4. To realize a setting with a
high-dimensional regime, for each $i=1,\ldots,303$, we generate an independent
random vector from $\mathcal{N}(\bm{0},\bm{I}_{198})$ and concatenate it with
the original $14$ variable to obtain the covariate vector $\bm{X}_{i}$. Hence,
we have $\kappa=(14+198)/303\approx 0.7$. We constructed CIs with $\alpha=0.1$
for $\bm{X}_{i}^{\top}\bm{\beta}$ by our proposed CI and a classical CI by the
MLE for Poisson regression (See Section A.6 for the construction of the
classical CI). We approximate the coverage proportion of CIs to
$\bm{X}_{i}^{\top}\bm{\beta}$ using $Y_{i}$, and preclude samples with
responses that take zero because they cannot be covered.

As a result, the coverage proportion of classical CIs is $139/139=100\%$ and
that of the proposed CIs is $126/139\approx 90.6\%$. Given that the ideal
coverage proportion is $1-\alpha=0.9$, the proposed CI provides more
appropriate coverage, whereas the classical CI is too conservative. This
result implies that classical CIs overfit the samples and underestimate the
uncertainty of the estimation, whereas the proposed CIs control the coverage
rate nearly at the preassigned level and can evaluate the uncertainty much
more accurately.

##  7\. Conclusion

We developed a statistical inference method for GLMs with asymmetric inverse
link functions in high dimensions. Our approach extends the SE-based inference
method designed for a logistic regression model of GLMs. Specifically, we
developed a surrogate estimator for GLMs, an associated SE system, and a
method to estimate the necessary parameters in the system. Our methodology
works well in terms of both theory and experimentation. One limitation of our
method is that it requires Gaussianity of the covariate; however, this
limitation can be relaxed by applying studies on the universality of
methodologies in terms of data, for example, Montanari and Saeed [36],
Venkataramanan et al. [49].

##  Appendix A Supportive Information

###  A.1. Construction of $\widehat{\tau}_{j}^{2}$

To construct the corrected confidence interval (14) with correlated features,
we have to estimate the conditional variance parameter,
$\tau_{j}^{2}=\mathrm{Var}(\bm{X}_{ij}|\bm{X}_{i\setminus j})$, for each
$j=1,\ldots,p$. In this article, we follow Zhao et al. [51]. They consider the
residual sum of squares $\mathrm{RSS}_{j}$ obtained by regression of
$\bm{X}_{j}=(\bm{X}_{1j},\ldots,\bm{X}_{nj})^{\top}\in\mathbb{R}^{n}$ onto
$\bm{X}_{\setminus
j}=(\bm{X}_{1},\ldots,\bm{X}_{j-1},\bm{X}_{j+1},\ldots,\bm{X}_{p})\in\mathbb{R}^{n\times(p-1)}$.
Owing to the Gaussianity of $\bm{X}_{i}$, it satisfies

$\displaystyle\mathrm{RSS}_{j}=\bm{X}_{j}^{\top}\bm{P}_{\bm{X}_{\setminus
j}}^{\perp}\bm{X}_{j}\sim\tau_{j}^{2}\chi_{n-p+1}^{2},$

where $\bm{P}_{\bm{X}_{\setminus j}}^{\perp}$ is the orthogonal projection
matrix onto the orthogonal complement of column space spanned by
$\bm{X}_{\setminus j}$. Then, we immediately obtain an unbiased estimator of
$\tau_{j}^{2}$:

$\displaystyle\widehat{\tau}_{j}^{2}=\frac{1}{n-p+1}\mathrm{RSS}_{j}.$

###  A.2. Master Theorem of Generalized Approximate Message Passing (GAMP)

We present the GAMP algorithm by Rangan [37] and its associated theoretical
result, which are key tools to inspect the limiting distributional behavior of
the estimator $\widehat{\bm{\beta}}$.

First, we provide the GAMP algorithm, which generate a sequence of parameters
$\tilde{\bm{\beta}}^{k}$ for an index $k\in\mathbb{N}\cup\\{0\\}$ based on the
minimization problem in (8) and a limit of the sequence corresponds to the
estimator $\widehat{\bm{\beta}}$. Let
$\bm{X}=(\bm{X}_{1},\ldots,\bm{X}_{n})^{\top}\in\mathbb{R}^{n\times p}$,
$\breve{\bm{X}}=\bm{X}/\sqrt{n}$, and
${\bm{Y}}=(Y_{1},\ldots,Y_{n})\in\mathbb{R}^{n}$. Given initial values
$\bar{\eta}_{0}>0,\tilde{\bm{\beta}}^{0}$, and
$\tilde{\bm{\xi}}^{0}=\breve{\bm{X}}\tilde{\bm{\beta}}^{0}$, a GAMP recursion
takes the form, for each $k\in\mathbb{N}\cup\\{0\\}$,

$\displaystyle\tilde{\bm{\beta}}^{k+1}$
$\displaystyle=\frac{\bar{\eta}_{k+1}}{\kappa}\breve{\bm{X}}^{\top}\left\\{\bm{Y}-g\left(\mathrm{prox}_{\bar{\eta}_{k}G}(\tilde{\bm{\xi}}^{k}+\bar{\eta}_{k}\bm{Y})\right)\right\\}+\frac{\bar{\eta}_{k+1}}{\bar{\eta}_{k}}\tilde{\bm{\beta}}^{k},$
(16) $\displaystyle\tilde{\bm{\xi}}^{k+1}$
$\displaystyle=\breve{\bm{X}}\tilde{\bm{\beta}}^{k+1}-\bar{\eta}_{k+1}\left\\{\bm{Y}-g\left(\mathrm{prox}_{\bar{\eta}_{k}G}(\tilde{\bm{\xi}}^{k}+\bar{\eta}_{k}\bm{Y})\right)\right\\}.$
(17)

Here, $\bar{\eta}_{k}$ is updated with $\bar{\mu}_{k}$ and
$\bar{\sigma}^{2}_{k}$ as following:

$\displaystyle\bar{\eta}_{k+1}$
$\displaystyle=\kappa\bar{\eta}_{k}\left(1-\mathbb{E}_{(Q_{1},Q_{2},U)}\left[\frac{1}{1+\bar{\eta}_{k}g^{\prime}\left(d_{k}\right)}\right]\right)^{-1},$
(18) $\displaystyle\bar{\mu}_{k+1}$
$\displaystyle=\frac{\bar{\eta}_{k+1}}{\gamma^{2}}\mathbb{E}_{(Q_{1},Q_{2},U)}\left[\gamma
Q_{1}\left\\{h(\gamma Q_{1},U)-g\left(d_{k}\right)\right\\}\right],$ (19)
$\displaystyle\bar{\sigma}_{k+1}^{2}$
$\displaystyle=\frac{\bar{\eta}_{k+1}^{2}}{\kappa^{2}}\mathbb{E}_{(Q_{1},Q_{2},U)}\left[\left\\{h(\gamma
Q_{1},U)-g\left(d_{k}\right)\right\\}^{2}\right],$ (20)

where

$\displaystyle
d_{k}=\mathrm{prox}_{\bar{\eta}_{k}(\gamma)G}(\bar{\mu}_{k}\gamma
Q_{1}+\sqrt{\kappa}\bar{\sigma}_{k}Q_{2}+\eta_{k}h(\gamma Q_{1},U)),$
$\displaystyle(Q_{1},Q_{2})\sim\mathcal{N}_{2}(\bm{0},\bm{I}_{2}),\qquad
U\perp\\!\\!\\!\perp(Q_{1},Q_{2}).$

Here, $U\in\mathbb{R}$ is a random variable independent of $(Q_{1},Q_{2})$.
This algorithm is proposed by Rangan [37]. It is closely related to the
linearized alternating direction method of multiplier [38].

Second, we provide the following result for the estimator
$\widehat{\bm{\beta}}$, which is regarded as a limit of the parameter sequence
generated by the GAMP algorithm. The result follows the master theorem for the
GAMP algorithm presented by Feng et al. [21].

######  Lemma A.1.

Assume that $\\{(\bm{X}_{i},Y_{i})\\}_{i=1}^{n}$ is an i.i.d. sample, and
$\bm{X}_{i}\sim\mathcal{N}_{p}(\bm{0},\bm{I}_{p})$ independent of $\bm{\beta}$
and $\bm{\varepsilon}=(\varepsilon_{1},\ldots,\varepsilon_{n})$. For
$r\in[2,\infty)$, suppose that the empirical distributions
$p^{-1}\sum_{j=1}^{p}\delta_{\sqrt{n}\beta_{j}}$ and
$n^{-1}\sum_{i=1}^{n}\delta_{\varepsilon_{i}}$ converge in the $r$-Wasserstein
sense to the distributions of $\bar{\beta}$ and $\bar{\varepsilon}$ with
finite $r$-th order moment, respectively. Let $Z\sim\mathcal{N}(0,\gamma^{2})$
be independent of $\bar{\varepsilon}$, and $G,\tilde{G}\sim\mathcal{N}(0,1)$
be independent of $\bar{\beta}$ and $Z$. Then, for finite
$\widehat{\bm{\beta}}$ defined in (8) and any pseudo-Lipschitz function
${\psi}:\mathbb{R}^{2}\to\mathbb{R}$ and
$\tilde{\psi}:\mathbb{R}^{3}\to\mathbb{R}$ of order $r$111A function
$\psi:\mathbb{R}^{m}\to\mathbb{R}$ is said to be pseudo-Lipschitz of order $r$
if there exists a constant $L>0$ such that for any
$t_{0},t_{1}\in\mathbb{R}^{m}$,
$\left\lVert\psi(t_{0})-\psi(t_{1})\right\rVert_{2}\leq L(1+\left\lVert
t_{0}\right\rVert_{2}^{r-1}+\left\lVert
t_{1}\right\rVert_{2}^{r-1})\left\lVert t_{0}-t_{1}\right\rVert$., we have

$\displaystyle\frac{1}{p}\sum_{j=1}^{p}{\psi}\left(\sqrt{n}\widehat{\beta}_{j},\sqrt{n}\beta_{j}\right)$
$\displaystyle\overset{\rm
a.s.}{\longrightarrow}\mathbb{E}\left[{\psi}\left(\mu\bar{\beta}+\sigma
G,\bar{\beta}\right)\right],$
$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\tilde{\psi}\left(\bm{X}_{i}^{\top}\widehat{\bm{\beta}},\bm{X}_{i}^{\top}\bm{\beta},\varepsilon_{i}\right)$
$\displaystyle\overset{\rm
a.s.}{\longrightarrow}\mathbb{E}\left[\tilde{\psi}\left(\mathrm{prox}_{\eta\ell}\left(\mu_{Z}Z+\sigma_{Z}\tilde{G}\right),Z,\bar{\varepsilon}\right)\right],$

as $n,p(n)\to\infty$ with $p(n)/n\to\kappa\in(0,1)$.

The convergence of pseudo-Lipschitz functions of order $1$ is equivalent to
the convergence of the $1$-Wasserstein distance between the empirical
distributions of each coordinate, as established through Kantorovich-
Rubinstein duality [28].

###### Proof of Lemma A.1.

We prove this lemma by the general proof strategy discussed in Section 4.4 in
Feng et al. [21]. It consists of three steps: (i) find a fixed point of the
GAMP recursion; (ii) consider the stationary version of the GAMP recursion;
(iii) show that the stationary version of the GAMP iterate converges to the
estimator $\widehat{\bm{\beta}}$. Following the strategy, the step (i) and the
step (ii) are simply achieved by the convexity and smoothness of the surrogate
loss function by the assumption. Given the fixed point $\bar{\eta}_{*}$ of
$\bar{\eta}_{k}$, and initial values
$\bm{\beta}^{0},\bm{\xi}^{0}=\breve{\bm{X}}\bm{\beta}^{0}$, the stationary
version of the GAMP algorithm takes the form, for each
$k\in\mathbb{N}\cup\\{0\\}$,

$\displaystyle\bm{\beta}^{k+1}$
$\displaystyle=\frac{\bar{\eta}_{*}}{\kappa}\breve{\bm{X}}^{\top}\left\\{\bm{Y}-g\left(\mathrm{prox}_{\bar{\eta}_{*}G}(\bm{\xi}^{k}+\bar{\eta}_{*}\bm{Y})\right)\right\\}+\bm{\beta}^{k}$
(21) $\displaystyle\bm{\xi}^{k+1}$
$\displaystyle=\breve{\bm{X}}\bm{\beta}^{k+1}-\bar{\eta}_{*}\left\\{\bm{Y}-g\left(\mathrm{prox}_{\bar{\eta}_{*}G}(\bm{\xi}^{k}+\bar{\eta}_{*}\bm{Y})\right)\right\\}.$
(22)

In the Step (iii), our goal is to ensure the algorithmic convergence of scaled
GAMP iterates $\widehat{\bm{\beta}}^{k}=\bm{\beta}^{k}/\sqrt{n}$ to the
estimator $\widehat{\bm{\beta}}$ as $k\to\infty$. Denote
$\ell(\bm{b})\equiv\sum_{i=1}^{n}\ell(\bm{b};\bm{X}_{i},Y_{i})$ for any
$\bm{b}\in\mathbb{R}^{p}$. By Taylor’s theorem, we have

$\displaystyle\ell(\widehat{\bm{\beta}})=\ell(\widehat{\bm{\beta}}^{k})+\left(\widehat{\bm{\beta}}-\widehat{\bm{\beta}}^{k}\right)^{\top}\nabla\ell(\widehat{\bm{\beta}}^{k})+\frac{1}{2}\left(\widehat{\bm{\beta}}-\widehat{\bm{\beta}}^{k}\right)^{\top}\nabla^{2}\ell\left(t\widehat{\bm{\beta}}+(1-t)\widehat{\bm{\beta}}^{k}\right)\left(\widehat{\bm{\beta}}-\widehat{\bm{\beta}}^{k}\right),$

for some $t\in(0,1)$. Thus, Lemma C.1 implies, for some non-increasing
positive function $\omega:\mathbb{R}_{+}\to\mathbb{R}_{+}$,

$\displaystyle\ell(\widehat{\bm{\beta}})\geq\ell(\widehat{\bm{\beta}}^{k})+\left(\widehat{\bm{\beta}}-\widehat{\bm{\beta}}^{k}\right)^{\top}\nabla\ell(\widehat{\bm{\beta}}^{k})+\frac{1}{2}n\omega\left(\max\left\\{\|\widehat{\bm{\beta}}\|,\|\widehat{\bm{\beta}}^{k}\|\right\\}\right)\|\widehat{\bm{\beta}}-\widehat{\bm{\beta}}^{k}\|^{2},$

with probability at least $1-c_{1}\exp(-c_{2}n)$ where $c_{1},c_{2}>0$ are
some positive constants. Here, we use
$\omega(\|t\widehat{\bm{\beta}}+(1-t)\widehat{\bm{\beta}}^{k}\|)\geq\omega(\max(\|\widehat{\bm{\beta}}\|,\|\widehat{\bm{\beta}}^{k}\|))$
since $\omega(\cdot)$ is non-increasing. Using optimality of the estimator
$\ell(\widehat{\bm{\beta}}^{k})\geq\ell(\widehat{\bm{\beta}})$ and the Cauchy-
Schwarz inequality as
$(\widehat{\bm{\beta}}-\widehat{\bm{\beta}}^{k})^{\top}\nabla\ell(\widehat{\bm{\beta}}^{k})\geq-\|{\nabla\ell(\widehat{\bm{\beta}}^{k})\|}\|\widehat{\bm{\beta}}-\widehat{\bm{\beta}}^{k}\|$
yields, with probability at least $1-c_{1}\exp(-c_{2}n)$,

$\displaystyle\|\widehat{\bm{\beta}}-\widehat{\bm{\beta}}^{k}\|\leq\frac{2}{\omega\left(\max\left\\{\|\widehat{\bm{\beta}}\|,\|\widehat{\bm{\beta}}^{k}\|\right\\}\right)}\left\lVert\frac{1}{n}\nabla\ell(\widehat{\bm{\beta}}^{k})\right\rVert\leq\frac{2}{\omega(\|\widehat{\bm{\beta}}\|)\omega(\|\widehat{\bm{\beta}}^{k}\|)}\left\lVert\frac{1}{n}\nabla\ell(\widehat{\bm{\beta}}^{k})\right\rVert,$

where the last inequality follows from the fact that $0<\omega(\cdot)<1$ and
$\omega$ is non-increasing. Next, we consider controlling
$\|\nabla\ell(\widehat{\bm{\beta}}^{k})\|$. We have

$\displaystyle\mathrm{prox}_{\bar{\eta}_{*}G}\left(\bm{\xi}^{k-1}+\bar{\eta}_{*}\bm{Y}\right)$
$\displaystyle=\bm{\xi}^{k-1}+\bar{\eta}_{*}\bm{Y}-\bar{\eta}_{*}g\left(\mathrm{prox}_{\bar{\eta}_{*}G}\left(\bm{\xi}^{k-1}+\bar{\eta}_{*}\bm{Y}\right)\right)$
(23) $\displaystyle=\bm{\xi}^{k-1}-\bm{\xi}^{k}+\breve{\bm{X}}\bm{\beta}^{k}$
(24)
$\displaystyle=\bm{\xi}^{k-1}-\bm{\xi}^{k}+\bm{X}\widehat{\bm{\beta}}^{k},$
(25)

by the definition of the proximal operator and (22). Thus,

$\displaystyle\bm{\beta}^{k}-\bm{\beta}^{k-1}$
$\displaystyle=\frac{\bar{\eta}_{*}}{\kappa}\breve{\bm{X}}^{\top}\left\\{\bm{Y}-g\left(\mathrm{prox}_{\bar{\eta}_{*}G}(\bm{\xi}^{k-1}+\bar{\eta}_{*}\bm{Y})\right)\right\\}$
$\displaystyle=\frac{\bar{\eta}_{*}}{\kappa}\breve{\bm{X}}^{\top}\left\\{\bm{Y}-g\left(\bm{X}\widehat{\bm{\beta}}^{k}+\bm{\xi}^{k-1}-\bm{\xi}^{k}\right)\right\\},$

by (21). Using this and triangle inequalities give

$\displaystyle\left\lVert\nabla\ell(\widehat{\bm{\beta}}^{k})\right\rVert$
$\displaystyle=\left\lVert\bm{X}^{\top}\left\\{g(\bm{X}^{\top}\widehat{\bm{\beta}}^{k})-\bm{Y}\right\\}\right\rVert$
$\displaystyle\leq\left\lVert\bm{X}^{\top}\left\\{\bm{Y}-g(\bm{X}\widehat{\bm{\beta}}^{k}+\bm{\xi}^{k-1}-\bm{\xi}^{k})\right\\}\right\rVert$
$\displaystyle\quad+\left\lVert\bm{X}^{\top}\left\\{g(\bm{X}\widehat{\bm{\beta}}^{k}+\bm{\xi}^{k-1}-\bm{\xi}^{k})-g(\bm{X}^{\top}\widehat{\bm{\beta}}^{k})\right\\}\right\rVert$
$\displaystyle\leq\frac{p}{\bar{\eta}_{*}}\left\lVert\bm{\beta}^{k}-\bm{\beta}^{k-1}\right\rVert+\left\lVert\bm{X}\right\rVert_{\rm
op}\left\lVert
g(\bm{X}\widehat{\bm{\beta}}^{k}+\bm{\xi}^{k-1}-\bm{\xi}^{k})-g(\bm{X}^{\top}\widehat{\bm{\beta}}^{k})\right\rVert$
$\displaystyle\leq\frac{p}{\bar{\eta}_{*}}\left\lVert\bm{\beta}^{k}-\bm{\beta}^{k-1}\right\rVert+L_{g}\left\lVert\bm{X}\right\rVert_{\rm
op}\left\lVert\bm{\xi}^{k-1}-\bm{\xi}^{k}\right\rVert,$

where $L_{g}=\sup_{z}g^{\prime}(z)$. This establishes, with probability at
least $1-c_{1}\exp(-c_{2}n)$,

$\displaystyle\|\widehat{\bm{\beta}}-\widehat{\bm{\beta}}^{k}\|\leq
c\left\\{\frac{\kappa}{\bar{\eta}_{*}}\left\lVert\widehat{\bm{\beta}}^{k}-\widehat{\bm{\beta}}^{k-1}\right\rVert+\frac{C_{g}}{n}\left\lVert\bm{X}\right\rVert_{\rm
op}\left\lVert\bm{\xi}^{k-1}-\bm{\xi}^{k}\right\rVert\right\\},$

with
$c={2}/\left(\omega(\|\widehat{\bm{\beta}}\|)\omega(\|\widehat{\bm{\beta}}^{k}\|)\right)$.
Finally, Lemma C.2 and the Borel-Cantelli lemma implies

$\displaystyle\lim_{k\to\infty}\lim_{n\to\infty}\left\lVert\widehat{\bm{\beta}}-\widehat{\bm{\beta}}^{k}\right\rVert=_{\rm
a.s.}0.$

This completes the proof. ∎

Lemma A.1 provides a statement on the convergence of the distance between
1-dimensional distributions, which is computed by averaging over the
coordinates of $\widehat{\bm{\beta}},\bm{\beta}\in\mathbb{R}^{p}$. However,
for statistical inference purposes, it is necessary to investigate the
limiting behavior of the marginal distributions of each coordinate. To achieve
this, we utilize a property of the surrogate loss function.

######  Lemma A.2.

For any invertible matrix $\bm{L}\in\mathbb{R}^{p\times p}$ satisfying
$\bm{\Sigma}=\bm{L}\bm{L}^{\top}$, $\bm{L}^{\top}\widehat{\bm{\beta}}$
minimizes the surrogate loss function in (8) for the true coefficient
$\bm{L}^{\top}\bm{\beta}$ and the covariate
$\bm{L}^{-1}\bm{X}\sim\mathcal{N}_{p}(\bm{0},\bm{I}_{p})$.

###### Proof of Lemma A.2.

Since the surrogate loss depends on $\bm{X}\in\mathbb{R}^{p}$ and
$\bm{b}\in\mathbb{R}^{p}$ only through their inner product
$\bm{X}^{\top}\bm{b}$, we have
$\ell(\bm{b};\bm{X},Y)=\ell(\bm{L}^{\top}\bm{b};\bm{L}^{-1}\bm{X},Y)$. ∎

Using this in reverse, once we show Lemma A.1 for a design with identity
covariance, we can rewrite the estimator corresponding to the unit covariance
in $\bm{L}^{\top}\widehat{\bm{\beta}}$ to obtain the general covariance
result.

####  A.2.1. Limit of Estimation and Classification Errors

As the consequence of the master theorem of GAMP (Lemma A.1), we obtain the
convergence limits of the mean squared error (MSE) and cosine similarity-like
classification error. If we set $\psi(s,t)=(s-t)^{2}$, then Lemma A.1 implies
the MSE limit

$\displaystyle\frac{1}{p}\|\widehat{\bm{\beta}}-\bm{\beta}\|_{2}^{2}\overset{\mathrm{a.s.}}{\longrightarrow}(\mu-1)^{2}\mathbb{E}[\bar{\beta}^{2}]+\sigma^{2}.$

By setting $\psi(s,t)=(s-\mu t)^{2}$, we have a corrected MSE limit

$\displaystyle\frac{1}{p}\|\widehat{\bm{\beta}}-\mu\bm{\beta}\|_{2}^{2}\overset{\mathrm{a.s.}}{\longrightarrow}\sigma^{2}.$

For classification errors, by taking the ratio of $\psi(s,t)=st$ and
$\psi(s,t)=t^{2}$, we can obtain the convergence limit of a cosine similarity-
like measure

$\displaystyle\frac{\widehat{\bm{\beta}}^{\top}\bm{\beta}}{\|\bm{\beta}\|^{2}}\overset{\mathrm{a.s.}}{\longrightarrow}\mu.$

###  A.3. Approximation of Inverse Link Functions

In this section, we explore a way to approximate an inverse link function that
does not satisfy Assumption 1 with an inverse link function that falls within
a manageable range.

Firstly, when dealing with a function $g(\cdot)$ that lacks Lipschitzness or
smoothness, we consider the approach of clipping either its value or
derivative by a suitably large threshold. This enables us to handle functions
that are otherwise difficult to work with. As an example, we approximate the
exponential inverse link function $g(t)=e^{t}$ using the following
formulation:

$\displaystyle\tilde{g}(t)=e^{t}1\\{t\leq\log(50)\\}+50(t+1-\log(50))1\\{t>\log(50)\\}.$

Here, $\tilde{g}$ satisfies $L^{\prime}$-Lipschitz continuity and
$L^{\prime}$-smoothness with $L^{\prime}=50$. In our numerical experiments
conducted in Section 6 on Poisson regression, we constructed an estimate for
the data generated from the inverse link function $g(t)=e^{t}$, assuming that
$\tilde{g}(\cdot)$ above serves as an approximation of the inverse link
function. This construction yielded satisfactory results.

Next, we consider the case where $g(\cdot)$ does not satisfy the symmetry
assumption, which is essential for the existence and uniqueness of our
estimation for the SS parameter $\gamma^{2}$. We consider a non-odd
approximation of an odd inverse link. Note that our approximation here
directly applies to the data-generating process. For example, if $g(\cdot)$ is
not exactly linear but follows the form $g(t)=\min(1.1t,0.9t)$, we can
consistently estimate the SS parameter for this model. The same holds true
when $g(\cdot)$ is not precisely sigmoid but instead follows
$g(t)=1/(1+\exp(-t+0.01))$.

###  A.4. Regularized Estimator

In this section, we consider the following regularized estimator:

$\displaystyle\widehat{\bm{\beta}}_{\lambda}\in\operatorname*{arg\,min}_{\bm{b}\in\mathbb{R}^{p}}\left\\{\sum_{i=1}^{n}\ell(\bm{b};\bm{X}_{i},Y_{i})+\lambda\sum_{j=1}^{p}J(b_{j})\right\\},$

where $J:\mathbb{R}\to\mathbb{R}$ is some regularization function and
$\lambda>0$ is a tuning parameter. It is important to consider the regularized
estimator $\widehat{\bm{\beta}}_{\lambda}$, because in some cases such as
logistic regression with $n<2p$, the unregularized estimator
$\widehat{\bm{\beta}}$ does not exist (see, for example, Candès and Sur [11]).
The regularization imposes constraints on the estimator within specific
regions around the origin. Consequently, the issue of non-existence can be
mitigated by employing a suitable regularized estimator. Note that even when
we employ the regularized estimator $\widehat{\bm{\beta}}_{\lambda}$, the
method of estimating $\gamma^{2}$ and $\sigma_{e}^{2}$ is not affected by the
method of estimating $\bm{\beta}$.

We display the system of nonlinear equations which characterizes the state
evolution parameters in $L_{2}$ penalized cases. If we set $J(t)=t^{2}$, we
have

$\displaystyle\begin{cases}\kappa^{2}\sigma^{2}&=\eta^{2}\mathbb{E}_{(Q_{1},Q_{2},\bar{Y})}\left[\left(\bar{Y}-g\left(L\right)\right)^{2}\right],\\\
2\gamma^{2}\lambda\mu&=\mathbb{E}_{(Q_{1},Q_{2},\bar{Y})}\left[Z\left(\bar{Y}-g\left(L\right)\right)\right],\\\
1-\kappa+2\lambda\eta&=\mathbb{E}_{(Q_{1},Q_{2},\bar{Y})}\left[({\displaystyle
1+\eta g^{\prime}\left(L\right)})^{-1}\right],\end{cases}$ (26)

where $L=\mathrm{prox}_{\eta G}(\mu Z+\sqrt{\kappa}\sigma Q_{2}+\eta\bar{Y})$,
$Z=\gamma Q_{1}$, $(Q_{1},Q_{2})\sim\mathcal{N}_{2}(0,\bm{I}_{2})$ and
$\bar{Y}=h(Z,\bar{\varepsilon})$. We can see that (26) admits the case $n<p$
since the left-hand side of the last equation can be positive with
sufficiently large $\lambda>0$ while the right-hand side is always positive.

Since regularized estimators are biased, the distributional characterization
of the limit of $\widehat{\bm{\beta}}_{\lambda}$ is somewhat different from
the unregularized case. Actually, we have an extension of Lemma A.1 as

$\displaystyle\frac{1}{p}\sum_{j=1}^{p}\psi\left(\sqrt{n}\widehat{\beta}_{j},\sqrt{n}\beta_{j}\right)\overset{\mathrm{a.s.}}{\longrightarrow}\mathbb{E}[\psi\left(\mathrm{prox}_{\eta
J}(\mu\bar{\beta}+\sigma G),\bar{\beta}\right)],$

under the settings of Lemma A.1. For the unregularized case
$J\equiv\mathrm{const.}$, this is reduced to Lemma A.1. Suppose that
$J(\cdot)$ is differentiable and strongly convex, and
$\bm{X}\sim\mathcal{N}_{p}(\bm{0},\bm{I}_{p})$ here. Using the fact that
$\mathrm{prox}_{\eta J}(x)=x-\eta J^{\prime}(\mathrm{prox}_{\eta J}(x))$ for
$x\in\mathbb{R}$ and $\eta>0$ by the definition of the proximal operator, we
finally obtain

$\displaystyle\sqrt{n}\frac{\widehat{\beta}_{j}^{\mathrm{(d)}}-\mu\beta_{j}}{\sigma}\overset{{d}}{\to}\mathcal{N}(0,1),$

where $\widehat{\bm{\beta}}^{\mathrm{(d)}}$ is a debiased estimator

$\displaystyle\widehat{\bm{\beta}}^{\mathrm{(d)}}=\widehat{\bm{\beta}}_{\lambda}+\eta
J^{\prime}(\widehat{\bm{\beta}}_{\lambda}).$

###  A.5. GLM Form with $h(Z,\bar{\varepsilon})$

In this section, we discuss the other form of GLM used in (A3) in Assumption 1
as follows:

$\displaystyle Y_{i}$
$\displaystyle=h\left(\bm{X}_{i}^{\top}\bm{\beta},\varepsilon_{i}\right),$
(27)

where $h:\mathbb{R}\times\mathbb{R}\to\mathbb{R}$ is a deterministic function
determined by the distribution of $Y\mid\bm{X}$, and
$\varepsilon_{i}\in\mathbb{R}$ is an error variable independent of $\bm{X}$.
While the GLM 1 only specifies the conditional mean, (27) completely
determines the distributional behavior of $Y$ depending on $\bm{X}$. The major
difference from the additive model (11) used in our estimation is that the
random variable $\varepsilon_{i}$ is independent to $\bm{X}_{i}^{\top}$.

This model (27) is flexible in design and allows for an intuitive
representation of GLMs. We see several examples of $h(Z,\bar{\varepsilon})$
introduced in (10).

$\bullet$ Bernoulli case (binary choice model). When we use a model
$\bar{Y}|Z\sim\mathrm{Ber}(g(Z))$ with some $g:\mathbb{R}\to[0,1]$, the
inverse transformation method yields

$\displaystyle\bar{Y}=1\\{g(Z)\leq\bar{\varepsilon}\\},\quad\bar{\varepsilon}\sim\mathrm{Unif}[0,1].$

$\bullet$ Exponential case. If $\bar{Y}|Z\sim\mathrm{Exp}(g(Z))$ with some
$g:\mathbb{R}\to(0,\infty)$, the inverse transformation method yields

$\displaystyle\bar{Y}=-\frac{1}{g(Z)}\log(\bar{\varepsilon}),\quad\bar{\varepsilon}\sim\mathrm{Unif}[0,1].$

$\bullet$ Poisson case. If $\bar{Y}|Z\sim\mathrm{Pois}(g(Z))$ with some
$g:\mathbb{R}\to(0,\infty)$, we have

$\displaystyle\bar{Y}=\min\left\\{k\in\mathbb{N}\cup\\{0\\}\left|\sum_{l=1}^{k+1}\bar{\varepsilon}_{l}>g(Z)\right.\right\\},\quad\bar{\varepsilon}_{l}\overset{\rm
iid}{\sim}\mathrm{Exp}(1).$

$\bullet$ Gaussian case. If
$\bar{Y}|Z\sim\mathcal{N}(g(Z),\sigma_{\bar{\varepsilon}}^{2})$ with some
$g:\mathbb{R}\to\mathbb{R}$, we have

$\displaystyle\bar{Y}=g(Z)+\bar{\varepsilon},\quad\bar{\varepsilon}\sim\mathcal{N}(0,\sigma_{\bar{\varepsilon}}^{2}).$

As evident from the aforementioned examples, in the case of modeling using a
one-parameter distribution, the distribution of $\bar{\varepsilon}$ can be
fully determined without any additional parameters.

###  A.6. Construction of CI in Real Data Application

This section specifies the construction of the classical and proposed CI in
Section 6.3. To begin with, we compute the MLE $\widehat{\bm{\beta}}$ for the
Poisson regression model. A classical MLE theory implies, for any
$\bm{x}\in\mathbb{R}^{p}$,

$\displaystyle\sqrt{n}(\bm{x}^{\top}\widehat{\bm{\beta}}-\bm{x}^{\top}\bm{\beta})\overset{{d}}{\to}\mathcal{N}(0,\bm{x}^{\top}\mathcal{I}_{\bm{\beta}}^{-1}\bm{x}),$

as $n\to\infty$ with fixed $p$. Using this, we construct the classical CI of
each $\bm{X}_{i}^{\top}\bm{\beta}$ with a preassigned level $(1-\alpha)$ as

$\displaystyle\left[\bm{X}_{i}^{\top}\widehat{\bm{\beta}}-z_{(1-\alpha/2)}\sqrt{\frac{{\bm{X}_{i}^{\top}{\widehat{\mathcal{I}}}(\widehat{\bm{\beta}})\bm{X}_{i}}}{n}},\bm{X}_{i}^{\top}\widehat{\bm{\beta}}+z_{(1-\alpha/2)}\sqrt{\frac{{\bm{X}_{i}^{\top}{\widehat{\mathcal{I}}}(\widehat{\bm{\beta}})\bm{X}_{i}}}{n}}\right].$

We have proposed a CI for $\beta_{j},j=1,\ldots,p$ in (14). Using a similar
technique, we can construct a valid CI for
$\bm{X}_{i}^{\top}\bm{\beta},i=1,\ldots,n$ in proportionally high dimensions.
Actually, we obtain

######  Proposition A.1.

Under the setting of Theorem 5.2, we have, for each $i=1,\ldots,n$,

$\displaystyle\frac{\bm{X}_{i}^{\top}\widehat{\bm{\beta}}+\eta\ell^{\prime}(\bm{X}_{i}^{\top}\widehat{\bm{\beta}})-\mu_{Z}\bm{X}_{i}^{\top}\bm{\beta}}{\sigma_{Z}}\overset{{d}}{\to}\mathcal{N}(0,1),$

as $n,p(n)\to\infty$ with $p(n)/n\to\kappa$. Here,
$\ell^{\prime}(\bm{X}_{i}^{\top}\widehat{\bm{\beta}})=g(\bm{X}_{i}^{\top}\widehat{\bm{\beta}})-Y_{i}$.

This proposition yields the asymptotically level $(1-\alpha)$ confidence
interval of $\bm{X}_{i}^{\top}\bm{\beta}$,

$\displaystyle\frac{1}{\mu_{Z}}\left[\bm{X}_{i}^{\top}\widehat{\bm{\beta}}+\eta\ell^{\prime}(\bm{X}_{i}^{\top}\widehat{\bm{\beta}})-\sigma_{Z}z_{(1-\alpha/2)},\bm{X}_{i}^{\top}\widehat{\bm{\beta}}+\eta\ell^{\prime}(\bm{X}_{i}^{\top}\widehat{\bm{\beta}})+\sigma_{Z}z_{(1-\alpha/2)}\right].$

###### Proof of Proposition A.1.

Since
$\bm{X}_{i}^{\top}\widehat{\bm{\beta}}=(\bm{L}^{-1}\bm{X}_{i})^{\top}(\bm{L}^{\top}\widehat{\bm{\beta}})$
with $\bm{L}^{-1}\bm{X}_{i}\sim\mathcal{N}_{p}(\bm{0},\bm{I}_{p})$, we can
repeat the arguments in Proposition 5.1 by replacing $\bm{\theta}$ and
$\widehat{\bm{\theta}}$ with
$(\bm{X}_{1}^{\top}\bm{\beta},\ldots,\bm{X}_{n}^{\top}\bm{\beta})^{\top}$ and
$(\bm{X}_{1}^{\top}\widehat{\bm{\beta}}+\eta\ell^{\prime}(\bm{X}_{1}^{\top}\widehat{\bm{\beta}}),\ldots,\bm{X}_{n}^{\top}\widehat{\bm{\beta}}+\eta\ell^{\prime}(\bm{X}_{n}^{\top}\widehat{\bm{\beta}}))$,
respectively. ∎

##  Appendix B Proofs of Main Results

###### Proof of Proposition 3.2.

Since the surrogate loss has an equivalent form of the negative log-likelihood
function for logistic regression, it can be derived according to the
discussion of Feng et al. [21] in Section 4.7. Note that the construction of
$h(\cdot,\cdot)$ and the inverse link function $g(\cdot)$ are generalized in
our case. ∎

###### Proof of Proposition 5.1.

Lemma A.2 implies that, for any $j=1,\ldots,p$,

$\displaystyle\frac{\widehat{\beta}_{j}-\mu\beta_{j}}{{\sigma}/\tau_{j}}=\frac{\widehat{\theta}_{j}-\mu\theta_{j}}{{\sigma}},$

where we define

$\displaystyle\bm{\theta}=\bm{L}^{\top}\bm{\beta},\quad\widehat{\bm{\theta}}=\bm{L}^{\top}\widehat{\bm{\beta}},$
(28)

by a Cholesky factorization $\bm{\Sigma}=\bm{L}\bm{L}^{\top}$. We have
$\widehat{\theta}_{j}=\tau_{j}\widehat{\beta}_{j}$ and
$\theta_{j}=\tau_{j}\beta_{j}$ by the rearrangement of indices. Here, define

$\displaystyle\mu_{n}=\frac{\widehat{\bm{\theta}}^{\top}\bm{\theta}}{\left\lVert\bm{\theta}\right\rVert^{2}},\qquad\mbox{and}\qquad\sigma_{n}^{2}=\frac{1}{\kappa}\|\widehat{\bm{\theta}}-\mu_{n}\bm{\theta}\|^{2}.$
(29)

Then, we have

$\displaystyle\sqrt{n}\frac{\widehat{\theta}_{j}-\mu\theta_{j}}{\sigma}=\sqrt{n}\frac{\widehat{\theta}_{j}-\mu_{n}\theta_{j}}{\sigma_{n}}\frac{\sigma_{n}}{\sigma}+\sqrt{n}\frac{(\mu_{n}-\mu)\theta_{j}}{\sigma}.$

Lemma C.4 gives us

$\displaystyle\frac{\widehat{\bm{\theta}}-\mu_{n}\bm{\theta}}{\sigma_{n}}\overset{\rm
d}{=}\frac{\bm{P}_{\bm{\theta}}^{\perp}\bm{Z}}{\|\bm{P}_{\bm{\theta}}^{\perp}\bm{Z}\|},$
(30)

where $\bm{Z}=(Z_{1},\ldots,Z_{p})\sim\mathcal{N}_{p}(\bm{0},\bm{I}_{p})$.
Triangle inequalities yield that

$\displaystyle\frac{\left\lVert\bm{Z}\right\rVert}{\sqrt{p}}-\frac{|\bm{\theta}^{\top}\bm{Z}|}{\sqrt{p}\left\lVert\bm{\theta}\right\rVert}\leq\frac{\|\bm{P}_{\bm{\theta}}^{\perp}\bm{Z}\|}{\sqrt{p}}\leq\frac{\left\lVert\bm{Z}\right\rVert}{\sqrt{p}}+\frac{|\bm{\theta}^{\top}\bm{Z}|}{\sqrt{p}\left\lVert\bm{\theta}\right\rVert}.$

Since
$|\bm{\theta}^{\top}\bm{Z}|/(\sqrt{p}\left\lVert\bm{\theta}\right\rVert)\overset{\mathrm{a.s.}}{\longrightarrow}0$
and
$\left\lVert\bm{Z}\right\rVert/\sqrt{p}\overset{\mathrm{a.s.}}{\longrightarrow}1$,
we have
$\|\bm{P}_{\bm{\theta}}^{\perp}\bm{Z}\|/\sqrt{p}\overset{\mathrm{a.s.}}{\longrightarrow}1$.
Then, this fact and (30) imply that

$\displaystyle\sqrt{n}\frac{\widehat{\theta}_{j}-\mu_{n}\theta_{j}}{\sigma_{n}}\overset{\rm
d}{=}\frac{1}{\sqrt{\kappa}}\sigma_{j}Q+o_{p}(1),\qquad\sigma_{j}^{2}=1-\frac{\theta_{j}^{2}}{\left\lVert\bm{\theta}\right\rVert^{2}},$

where $Q\sim\mathcal{N}(0,1)$. Here we use the fact that the covariance matrix
of $\bm{P}_{\bm{\theta}}^{\perp}\bm{Z}$ is
$\bm{P}_{\bm{\theta}}^{\perp}\bm{P}_{\bm{\theta}}^{\perp}=\bm{I}_{p}-\bm{\theta}\bm{\theta}^{\top}/\|\bm{\theta}\|^{2}$.
Thus, the facts that $\mu_{n}\overset{\mathrm{a.s.}}{\longrightarrow}\mu$ and
$\sigma_{n}^{2}\overset{\mathrm{a.s.}}{\longrightarrow}\sigma^{2}$ by Lemma
C.5 conclude the proof. ∎

###### Proof of Lemma 5.1.

Note that a density function of Gaussian distribution with mean zero is an
even function. Since a product of even functions is even and a product of an
even function and an odd function is odd, we have

$\displaystyle\mathbb{E}[g(Z_{\varsigma})]=\frac{1}{\varsigma}\int_{-\infty}^{\infty}g(x)\phi\left(\frac{x}{\varsigma}\right)dx=\int_{-\infty}^{\infty}g(\varsigma
y)\phi\left(y\right)dy=2\int_{0}^{\infty}g_{0}(\varsigma
y)\phi\left(y\right)dy,$

where the second identity is from a change of variables $y=x/\varsigma$. Thus,
For any $\varsigma_{1}>\varsigma_{0}>0$, we can say that

$\displaystyle\mathbb{E}[g(Z_{\varsigma_{1}})]-\mathbb{E}[g(Z_{\varsigma_{0}})]=2\int_{0}^{\infty}\left\\{g_{0}(\varsigma_{1}y)-g_{0}(\varsigma_{0}y)\right\\}\phi\left(y\right)dy.$

is strictly positive or negative by $\phi(\cdot)>0$. Therefore, if the even
part $g_{0}(\cdot)$ of $g(\cdot)$ is monotone on $\mathbb{R}_{+}$, then
$\mathbb{E}[g(Z_{\varsigma})]$ is monotone in $\varsigma>0$. ∎

###### Proof of Theorem 5.1.

We assume that the simulated mean satisfies
$\widehat{\mathbb{E}}[\cdot]=\mathbb{E}[\cdot]$. The law of large numbers
yields

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}Y_{i}=\frac{1}{n}\sum_{i=1}^{n}\left\\{g(\bm{X}_{i}^{\top}\bm{\beta})+e_{i}\right\\}\overset{\rm
a.s.}{\longrightarrow}\mathbb{E}[g(Z)],\qquad Z\sim\mathcal{N}(0,\gamma^{2}),$
(31)

as $n,p(n)\to\infty$ with $p(n)/n\to\kappa\in(0,\infty)$. Here, we use the
fact that $\mathbb{E}[e_{i}]=\mathbb{E}[\mathbb{E}[e_{i}|\bm{X}_{i}]]=0$. Let
$\Psi_{n}(\varsigma^{2})=n^{-1}\sum_{i=1}^{n}\left\\{Y_{i}-{\mathbb{E}}[g(Z_{\varsigma})]\right\\}$
and $\Psi(\varsigma^{2})=\mathbb{E}[\Psi_{n}(\varsigma^{2})]$. Then, we have

$\displaystyle\sup_{\varsigma^{2}>0}\left\lvert\Psi_{n}(\varsigma^{2})-\Psi(\varsigma^{2})\right\rvert=\left\lvert\frac{1}{n}\sum_{i=1}^{n}Y_{i}-\mathbb{E}[Y_{1}]\right\rvert\overset{\rm
a.s.}{\longrightarrow}0,$ (32)

by the strong law of large numbers. Since $\mathbb{E}[g(Z_{\varsigma})]$ is a
strictly monotone function of $\varsigma>0$ by the assumption of this theorem,
for every $\epsilon>0$, we have

$\displaystyle\inf_{\varsigma^{2}:|\varsigma^{2}-\gamma^{2}|>\epsilon}\left\lvert\Psi(\varsigma^{2})\right\rvert>0.$
(33)

Also, (31) and (32) directly implies

$\displaystyle\left\lvert\Psi(\gamma)\right\rvert=\left\lvert\mathbb{E}[Y_{1}]-\mathbb{E}[g(Z_{\gamma})]\right\rvert\leq\left\lvert\mathbb{E}[Y_{1}]-n^{-1}\sum
Y_{i}\right\rvert+\left\lvert n^{-1}\sum
Y_{i}-\mathbb{E}[g(Z)]\right\rvert\overset{\mathrm{a.s.}}{\longrightarrow}0,$
(34)

as $n\to\infty$. Since $\Psi(\gamma)$ does not depend on $n$, one can see that

$\displaystyle\Psi(\gamma)=0.$ (35)

Finally, (32), (33), and (35) together imply $\widehat{\gamma}^{2}\overset{\rm
a.s.}{\longrightarrow}\gamma^{2}$ by Theorem 5.9 in van der Vaart [48]. ∎

###### Proof of Proposition 5.2.

We assume that the simulated mean satisfies
$\widehat{\mathbb{E}}[\cdot]=\mathbb{E}[\cdot]$. The law of large numbers
implies

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}Y_{i}^{2}=\frac{1}{n}\sum_{i=1}^{n}\left\\{g(\bm{X}_{i}^{\top}\bm{\beta})+e_{i}\right\\}^{2}\overset{\rm
a.s.}{\longrightarrow}\mathbb{E}[g(Z)^{2}]+\sigma_{e}^{2},$

as $n,p(n)\to\infty$ with $p(n)/n\to\kappa\in(0,\infty)$. Here we use the
independence between $Z$ and $\bar{\varepsilon}$. By Theorem 5.1 and the
continuous mapping theorem, $g(Z_{\widehat{\gamma}})\overset{\rm
a.s.}{\longrightarrow}g(Z)$ by the fact that $Z=\gamma Z_{1}$. Thus,
$\mathbb{E}[g(Z_{\widehat{\gamma}})]\to\mathbb{E}[g(Z)]$ by Theorem 2 (b) in
Ferguson [22]. At last, we obtain

$\displaystyle\widehat{\sigma}_{e}^{2}=\left(\frac{1}{n}\sum_{i=1}^{n}Y_{i}^{2}\right)-{\mathbb{E}}\left[g\left(Z_{\widehat{\gamma}}\right)^{2}\right]\overset{\rm
a.s.}{\longrightarrow}\mathbb{E}[g(Z)^{2}]+\sigma_{e}^{2}-\mathbb{E}[g(Z)^{2}]=\sigma_{e}^{2}.$

∎

###### Proof of Theorem 5.2.

By Proposition 5.3, $\mathcal{CI}_{\alpha,j}$ is clearly asymptotically level
$(1-\alpha)$ confidence interval, since we have

$\displaystyle\mathrm{Pr}\left(\frac{\widehat{\beta}_{j}}{\widehat{\mu}}+z_{\alpha/2}\frac{\widehat{\sigma}}{\sqrt{n}\widehat{\mu}\widehat{\tau}_{j}}\leq\beta_{j}\leq\frac{\widehat{\beta}_{j}}{\widehat{\mu}}-z_{\alpha/2}\frac{\widehat{\sigma}}{\sqrt{n}\widehat{\mu}\widehat{\tau}_{j}}\right)$
$\displaystyle=\mathrm{Pr}\left(z_{\alpha/2}\leq\frac{\sqrt{n}(\widehat{\beta}_{j}-\widehat{\mu}\beta_{j})}{\widehat{\sigma}/\widehat{\tau}_{j}}\leq
z_{(1-\alpha/2)}\right)$ $\displaystyle\to(1-\alpha).$

∎

###### Proof of Proposition 5.3.

We have

$\displaystyle\sqrt{n}\widehat{\tau}_{j}\frac{\widehat{\beta}_{j}-\widehat{\mu}\beta_{j}}{\widehat{\sigma}}$
$\displaystyle=\sqrt{n}\widehat{\tau}_{j}\frac{\widehat{\beta}_{j}-\mu\beta_{j}}{{\sigma}}\frac{\sigma}{\widehat{\sigma}}-\sqrt{n}\widehat{\tau}_{j}\frac{(\mu-\widehat{\mu})\beta_{j}}{\widehat{\sigma}}$
$\displaystyle=\sqrt{n}\widehat{\tau}_{j}\frac{\widehat{\beta}_{j}-\mu\beta_{j}}{{\sigma}}-\sqrt{n}\widehat{\tau}_{j}\frac{\widehat{\beta}_{j}-\mu\beta_{j}}{{\sigma}}\frac{\widehat{\sigma}-\sigma}{\widehat{\sigma}}-\sqrt{n}\widehat{\tau}_{j}\frac{(\mu-\widehat{\mu})\beta_{j}}{\widehat{\sigma}}$
$\displaystyle=\sqrt{n}\tau_{j}\frac{\widehat{\beta}_{j}-\mu\beta_{j}}{{\sigma}}-o_{p}(1),$

where the last identity follows from Lemma C.3,
$\sqrt{n}\tau_{j}\beta_{j}=O(1)$, and the assumption (A1). Finally,
Proposition 5.1 concludes the proof. ∎

##  Appendix C Technical Lemmas

For a matrix $A\in\mathbb{R}^{m\times m}$, $A\succeq 0$ is defined to mean
that $A$ is positive semi-definite.

######  Lemma C.1.

Suppose that $\bm{X}_{i}\sim\mathcal{N}_{p}(\bm{0},\bm{I}_{p})$. Under the
assumption (A2), for some constant $\epsilon_{0}$ such that
$0\leq\epsilon\leq\epsilon_{0}$,

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\nabla^{2}\ell(\bm{b};\bm{X}_{i},Y_{i})\succeq\left(\inf_{z:|z|\leq\frac{3\left\lVert\bm{b}\right\rVert}{\sqrt{\epsilon}}}g^{\prime}(z)\right)\left(\sqrt{1-\epsilon}-\sqrt{\kappa}-2\sqrt{\frac{H(\epsilon)}{1-\epsilon}}\right)^{2}\bm{I}_{p}$

with $H(\epsilon)=-\epsilon\log\epsilon-(1-\epsilon)\log(1-\epsilon)$ holds
for any $\bm{b}\in\mathbb{R}^{p}$ with probability at least
$1-2\exp(-nH(\epsilon))-2\exp(-n/2)$.

###### Proof of Lemma C.1.

Since the proof of Lemma 3 in Sur et al. [42] only uses the specific structure
$\nabla^{2}\ell(\bm{b};\bm{X}_{i},Y_{i})=\rho^{\prime\prime}(\bm{X}_{i}^{\top}\bm{b})\bm{X}_{i}\bm{X}_{i}^{\top}$
with $\rho^{\prime}(t)=1/(1+\exp(-t))$, we can immediately extend it to our
general structure
$\nabla^{2}\ell(\bm{b};\bm{X}_{i},Y_{i})=g^{\prime}(\bm{X}_{i}^{\top}\bm{b})\bm{X}_{i}\bm{X}_{i}^{\top}$.
∎

######  Remark C.1.

This lemma implies that, for sufficiently small $\epsilon>0$,

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\nabla^{2}\ell(\bm{b};\bm{X}_{i},Y_{i})\succeq\omega(\left\lVert\bm{b}\right\rVert)\bm{I}_{p},$

for some non-increasing positive function
$\omega:\mathbb{R}_{+}\to\mathbb{R}_{+}$ independent of $n$.

######  Lemma C.2.

Consider the setting of Lemma A.1. For the generalized approximate message
passing recursion (21)-(22), we have the Cauchy property of the recursion:

$\displaystyle\lim_{k\to\infty}\lim_{n\to\infty}\|\bm{\beta}^{k+1}-\bm{\beta}^{k}\|^{2}=_{\rm
a.s.}0,$
$\displaystyle\lim_{k\to\infty}\lim_{n\to\infty}\|\bm{\xi}^{k+1}-\bm{\xi}^{k}\|^{2}=_{\rm
a.s.}0.$

###### Proof of Lemma C.2.

For the stationary version of the GAMP recursion, we have

$\displaystyle\bar{\eta}_{*}\left\\{\bm{Y}-g\left(\mathrm{prox}_{\bar{\eta}_{*}G}(\bm{\xi}^{k}+\bar{\eta}_{*}\bm{Y})\right)\right\\}=\mathrm{prox}_{\bar{\eta}_{*}\ell}(\bm{\xi}^{k})-\bm{\xi}^{k},$

where $\mathrm{prox}_{\bar{\eta}_{*}\ell}(\cdot)$ depends on the first input
of the surrogate loss function $\ell$, by the fact that $\mathrm{prox}_{\eta
G}(x)=x-g(\mathrm{prox}_{\eta G}(x))$ for any $x\in\mathbb{R},\eta>0$ by the
definition of the proximal operator. Thus, taking
$\Psi(z;b)=z-\mathrm{prox}_{b\ell}(z)$, we can straightforwardly repeat the
argument of Lemma 6.9 in Donoho and Montanari [15] and complete the proof. ∎

Recall that $\widehat{\eta}$ is a solution of the SE system as defined in
Section 5.1.

######  Lemma C.3.

Under the settings in Proposition 5.3, we have

$\displaystyle\widehat{\eta}\overset{{p}}{\to}\eta,\quad\widehat{\mu}\overset{{p}}{\to}\mu,\quad\widehat{\sigma}^{2}\overset{{p}}{\to}\sigma^{2},$

as $n,p(n)\to\infty$ with $p(n)/n\to\kappa$.

###### Proof of Lemma C.3.

In this proof, we denote $L_{f}$ as the Lipschitz constant of a function
$f(\cdot)$. To emphasize the dependence on the signal strength, we denote
$\eta(\widehat{\gamma})\equiv\widehat{\eta},\mu(\widehat{\gamma})\equiv\widehat{\mu}$,
$\sigma(\widehat{\gamma})\equiv\widehat{\sigma}$, $\eta(\gamma)\equiv\eta$,
$\mu(\gamma)\equiv\mu$, and $\sigma(\gamma)\equiv\sigma$. To begin with, note
that, as discussed in Section 4.4 in Feng et al. [21], the solutions to the
system of nonlinear equations (10) can be rewritten as fixed points of the
following recursions:

$\displaystyle\eta_{k+1}(\gamma)$
$\displaystyle=\kappa\eta_{k}(\gamma)\left(1-\mathbb{E}_{(Q_{1},Q_{2},U)}\left[\frac{1}{1+\eta_{k}(\gamma)g^{\prime}\left(d_{k}(\gamma)\right)}\right]\right)^{-1}$
(36) $\displaystyle\mu_{k+1}(\gamma)$
$\displaystyle=\frac{\eta_{k+1}(\gamma)}{\gamma^{2}}\mathbb{E}_{(Q_{1},Q_{2},U)}\left[\gamma
Q_{1}\left\\{h(\gamma Q_{1},U)-g\left(d_{k}(\gamma)\right)\right\\}\right]$
(37) $\displaystyle\sigma_{k+1}^{2}(\gamma)$
$\displaystyle=\frac{\eta_{k+1}^{2}(\gamma)}{\kappa^{2}}\mathbb{E}_{(Q_{1},Q_{2},U)}\left[\left\\{h(\gamma
Q_{1},U)-g\left(d_{k}(\gamma)\right)\right\\}^{2}\right],$ (38)

where

$\displaystyle
d_{k}(\gamma)=\mathrm{prox}_{\eta_{k}(\gamma)G}(\mu_{k}(\gamma)\gamma
Q_{1}+\sqrt{\kappa}\sigma_{k}(\gamma)Q_{2}+\eta_{k}(\gamma)h(\gamma
Q_{1},U)),$
$\displaystyle(Q_{1},Q_{2})\sim\mathcal{N}_{2}(\bm{0},\bm{I}_{2}),\qquad
U\perp\\!\\!\\!\perp(Q_{1},Q_{2}),$

and $h(\gamma Q_{1},U)$ is designed to have the same distribution as $Y$.

Step 1. In this step, we inductively show that in each step $k\in\mathbb{N}$,
$(\eta_{k}(\widehat{\gamma}),\mu_{k}(\widehat{\gamma}),\sigma_{k}(\widehat{\gamma}))$
of the recursion converges to
$(\eta_{k}({\gamma}),\mu_{k}({\gamma}),\sigma_{k}({\gamma}))$ using the fact
that $\widehat{\gamma}^{2}\overset{\rm a.s.}{\longrightarrow}\gamma^{2}$ by
Theorem 5.1. Let $(\eta_{0},\mu_{0},\sigma_{0})$ be a given triplet of the
initializers, and assume $\eta_{0}>0$ and
$\max\\{\left\lvert\eta_{k-1}(\widehat{\gamma})-\eta_{k-1}(\gamma)\right\rvert,\left\lvert\mu_{k-1}(\widehat{\gamma})-\mu_{k-1}(\gamma)\right\rvert,\left\lvert\sigma_{k-1}(\widehat{\gamma})-\sigma_{k-1}(\gamma)\right\rvert\\}=o_{p}(1)$.

$\bullet$ Bound for $d_{0}$. We have

$\displaystyle|d_{0}(\widehat{\gamma})-d_{0}(\gamma)|$
$\displaystyle=\left\lvert\mathrm{prox}_{\eta_{0}G}(\mu_{0}\widehat{\gamma}Q_{1}+\sqrt{\kappa}\sigma_{0}Q_{2}+\eta_{0}h(\widehat{\gamma}Q_{1},U))-\mathrm{prox}_{\eta_{0}G}(\mu_{0}\gamma
Q_{1}+\sqrt{\kappa}\sigma_{0}Q_{2}+\eta_{0}h(\gamma Q_{1},U))\right\rvert$
$\displaystyle\leq\left\lvert\mu_{0}(\widehat{\gamma}-\gamma)Q_{1}+\eta_{0}h(\widehat{\gamma}Q_{1},U)-\eta_{0}h(\gamma
Q_{1},U)\right\rvert$
$\displaystyle\leq\left((\mu_{0}+L_{h}\eta_{0})Q_{1}\right)\left\lvert\widehat{\gamma}-\gamma\right\rvert:=L_{d_{0}}\left\lvert\widehat{\gamma}-\gamma\right\rvert,$

where the first inequality follows from the Lipschitz continuous of the
proximal operator by Lemma C.6 (i), and the last inequality is from the
Lipschitz condition on $h(\cdot,\cdot)$ and the triangle inequality. Thus,
$\left\lvert d_{0}(\widehat{\gamma})-d_{0}(\gamma)\right\rvert=o_{p}(1)$

$\bullet$ Bound for $\eta_{1}$. We use the fact that
$\left\lvert\frac{1}{1-a}-\frac{1}{1-b}\right\rvert\leq\frac{1}{(1-a)(1-b)}\left\lvert
a-b\right\rvert$ for $0<a,b<1$, and
$\left\lvert\frac{1}{1+a}-\frac{1}{1+b}\right\rvert\leq\left\lvert
a-b\right\rvert$ for $0<a,b$. Define a constant

$\displaystyle
C_{\eta_{1}}=\left(1-\mathbb{E}\left[\frac{1}{1+\eta_{0}g^{\prime}(d_{0}(\widehat{\gamma}))}\right]\right)^{-1}\left(1-\mathbb{E}\left[\frac{1}{1+\eta_{0}g^{\prime}(d_{0}({\gamma}))}\right]\right)^{-1}>0.$

Since $1+\eta_{0}g^{\prime}(\cdot)>1$ by the monotonically increasing property
of $g(\cdot)$, we have

$\displaystyle\left\lvert\eta_{1}(\widehat{\gamma})-\eta_{1}(\gamma)\right\rvert$
$\displaystyle\leq\kappa\eta_{0}C_{\eta_{1}}\mathbb{E}\left[\left\lvert\frac{1}{1+\eta_{0}g^{\prime}\left(d_{0}(\widehat{\gamma})\right)}-\frac{1}{1+\eta_{0}g^{\prime}\left(d_{0}({\gamma})\right)}\right\rvert\right]$
$\displaystyle\leq\kappa\eta_{0}^{2}C_{\eta_{1}}\mathbb{E}\left[\left\lvert
g^{\prime}(d_{0}(\widehat{\gamma}))-g^{\prime}(d_{0}({\gamma}))\right\rvert\right]$
$\displaystyle\leq\kappa\eta_{0}^{2}C_{\eta_{1}}L_{g^{\prime}}L_{d_{0}}\left\lvert\widehat{\gamma}-\gamma\right\rvert,$

where the last inequality uses the $L_{g^{\prime}}$-smoothness of the inverse
link function $g(\cdot)$ and the Lipschitz continuity of $d_{0}(\gamma)$.
Then,
$\left\lvert\eta_{1}(\widehat{\gamma})-\eta_{1}(\gamma)\right\rvert=o_{p}(1).$

$\bullet$ Bound for $\mu_{1}$. By the triangle inequality, we have

$\displaystyle\left\lvert\mu_{1}(\widehat{\gamma})-\mu(\gamma)\right\rvert$
$\displaystyle\leq\left\lvert\left(\frac{\eta_{1}(\widehat{\gamma})}{\widehat{\gamma}}-\frac{\eta_{1}(\gamma)}{\gamma}\right)\mathbb{E}\left[Q_{1}\left\\{h(\widehat{\gamma}Q_{1},U)-g(d_{0}(\widehat{\gamma}))\right\\}\right]\right\rvert$
$\displaystyle\quad+\left\lvert\frac{\eta_{1}(\gamma)}{\gamma}\mathbb{E}\left[Q_{1}\left\\{h(\widehat{\gamma}Q_{1},U)-h(\gamma
Q_{1},U)-g(d_{0}(\widehat{\gamma}))+g(d_{0}(\gamma))\right\\}\right]\right\rvert.$

Using the fact that $\frac{c}{a+b}=\frac{c}{a}-\frac{cb}{a(a+b)}$ for any
$a\neq 0,b\neq-a,c\in\mathbb{R}$, we have

$\displaystyle\frac{\eta_{1}(\widehat{\gamma})}{\widehat{\gamma}}-\frac{\eta_{1}(\gamma)}{\gamma}=\frac{\eta_{1}(\widehat{\gamma})}{{\gamma}}-\frac{\eta_{1}(\widehat{\gamma})\left(\widehat{\gamma}-\gamma\right)}{\gamma\widehat{\gamma}}-\frac{\eta_{1}(\gamma)}{\gamma}=\frac{\eta_{1}(\widehat{\gamma})-\eta_{1}(\gamma)}{\gamma}-o_{p}(1)=o_{p}(1).$

Thus, by the Lipschitz continuity of $h(\cdot,\cdot)$, $g(\cdot)$, and
$d_{0}$, we have

$\displaystyle\left\lvert\mu_{1}(\widehat{\gamma})-\mu_{1}(\gamma)\right\rvert\leq\left\lvert\frac{\eta_{1}(\gamma)}{\gamma}\left(\widehat{\gamma}-\gamma\right)\mathbb{E}\left[Q_{1}(L_{h}Q_{1}-L_{g}L_{d_{0}})\right]\right\rvert+o_{p}(1)=o_{p}(1).$

$\bullet$ Bound for $\sigma_{1}$. By the triangle inequality and Jensen’s
inequality, we have

$\displaystyle\left\lvert\sigma_{1}^{2}(\widehat{\gamma})-\sigma_{1}^{2}(\gamma)\right\rvert$
$\displaystyle\leq\frac{1}{\kappa^{2}}\left\lvert\left(\eta_{1}^{2}(\widehat{\gamma})-\eta_{1}^{2}(\gamma)\right)\mathbb{E}\left[\left\\{h(\widehat{\gamma}Q_{1},U)-g\left(d_{0}(\widehat{\gamma})\right)\right\\}^{2}\right]\right\rvert$
$\displaystyle\quad+\frac{\eta_{1}^{2}(\gamma)}{\kappa^{2}}\mathbb{E}\left\lvert\left\\{h(\widehat{\gamma}Q_{1},U)-g\left(d_{0}(\widehat{\gamma})\right)\right\\}^{2}-\left\\{h({\gamma}Q_{1},U)-g\left(d_{0}({\gamma})\right)\right\\}^{2}\right\rvert$
$\displaystyle\leq\frac{1}{\kappa^{2}}\left\lvert\widehat{\gamma}-\gamma\right\rvert{\left(\eta_{1}(\widehat{\gamma})+\eta_{1}(\gamma)\right)\mathbb{E}\left[\left\\{h(\widehat{\gamma}Q_{1},U)-g\left(d_{0}(\widehat{\gamma})\right)\right\\}^{2}\right]}$
$\displaystyle\quad+\frac{\eta_{1}^{2}(\gamma)}{\kappa^{2}}\mathbb{E}\left[\left\lvert\left\\{h(\widehat{\gamma}Q_{1},U)-g\left(d_{0}(\widehat{\gamma})\right)\right\\}-\left\\{h({\gamma}Q_{1},U)-g\left(d_{0}({\gamma})\right)\right\\}\right\rvert\right.$
$\displaystyle\quad\cdot\left.\left\\{h(\widehat{\gamma}Q_{1},U)-g\left(d_{0}(\widehat{\gamma})\right)+h({\gamma}Q_{1},U)-g\left(d_{0}({\gamma})\right)\right\\}\right]$
$\displaystyle\leq
o_{p}(1)+\frac{\eta_{1}^{2}(\gamma)}{\kappa^{2}}\left\lvert\widehat{\gamma}-\gamma\right\rvert\mathbb{E}\left\lvert(L_{h}Q_{1}-L_{g}L_{d_{0}})\left(h(\widehat{\gamma}Q_{1},U)-g\left(d_{0}(\widehat{\gamma})\right)+h({\gamma}Q_{1},U)-g\left(d_{0}({\gamma})\right)\right)\right\rvert$
$\displaystyle=o_{p}(1).$

Also, since $|a-b|\leq\sqrt{|a^{2}-b^{2}|}$ for $a,b>0$, we have
$\left\lvert\sigma_{1}(\widehat{\gamma})-\sigma_{1}(\gamma)\right\rvert=o_{p}(1)$.

$\bullet$ Bound for $d_{k-1}$. By the triangle inequality, we have

$\displaystyle|d_{k-1}(\widehat{\gamma})-d_{k-1}(\gamma)|$
$\displaystyle=\left|\mathrm{prox}_{\eta_{k-1}(\widehat{\gamma})G}(\mu_{k-1}(\widehat{\gamma})\widehat{\gamma}Q_{1}+\sqrt{\kappa}\sigma_{k-1}(\widehat{\gamma})Q_{2}+\eta_{k-1}(\widehat{\gamma})h(\widehat{\gamma}Q_{1},U))\right.$
$\displaystyle\quad\left.-\mathrm{prox}_{\eta_{k-1}(\gamma)G}(\mu_{k-1}(\gamma)\gamma
Q_{1}+\sqrt{\kappa}\sigma_{k-1}(\gamma)Q_{2}+\eta_{k-1}(\gamma)h(\gamma
Q_{1},U))\right|$
$\displaystyle\leq\left|\mathrm{prox}_{\eta_{k-1}(\widehat{\gamma})G}\left(\mu_{k-1}(\widehat{\gamma})\widehat{\gamma}Q_{1}+\sqrt{\kappa}\sigma_{k-1}(\widehat{\gamma})Q_{2}+\eta_{k-1}(\widehat{\gamma})h(\widehat{\gamma}Q_{1},U)\right)\right.$
$\displaystyle\quad\left.-\mathrm{prox}_{\eta_{k-1}(\gamma)G}\left(\mu_{k-1}(\widehat{\gamma})\widehat{\gamma}Q_{1}+\sqrt{\kappa}\sigma_{k-1}(\widehat{\gamma})Q_{2}+\eta_{k-1}(\widehat{\gamma})h(\widehat{\gamma}Q_{1},U)\right)\right|$
$\displaystyle\quad+\left|\mathrm{prox}_{\eta_{k-1}(\gamma)G}\left(\mu_{k-1}(\widehat{\gamma})\widehat{\gamma}Q_{1}+\sqrt{\kappa}\sigma_{k-1}(\widehat{\gamma})Q_{2}+\eta_{k-1}(\widehat{\gamma})h(\widehat{\gamma}Q_{1},U)\right)\right.$
$\displaystyle\quad\left.-\mathrm{prox}_{\eta_{k-1}(\gamma)G}\left(\mu_{k-1}(\gamma)\gamma
Q_{1}+\sqrt{\kappa}\sigma_{k-1}(\gamma)Q_{2}+\eta_{k-1}(\gamma)h(\gamma
Q_{1},U)\right)\right|$ $\displaystyle\leq
C\left\lvert\eta_{k-1}(\widehat{\gamma})-\eta_{k-1}(\gamma)\right\rvert+\left\lvert
Q_{1}\left(\mu_{k-1}(\widehat{\gamma})\widehat{\gamma}-\mu_{k-1}(\gamma)\gamma\right)\right\rvert+\sqrt{\kappa}\left\lvert
Q_{2}(\sigma_{k-1}\left(\widehat{\gamma})-\sigma_{k-1}(\gamma)\right)\right\rvert$
$\displaystyle\quad+\left\lvert\eta_{k-1}(\widehat{\gamma})h(\widehat{\gamma}Q_{1},U)-\eta_{k-1}({\gamma})h({\gamma}Q_{1},U)\right\rvert,$

where $C$ is some positive constant, and the last inequality follows from
Lemma C.7 and Lemma C.6 (i). Using the triangle inequality again, we obtain

$\displaystyle|d_{k-1}(\widehat{\gamma})-d_{k-1}(\gamma)|$
$\displaystyle\leq\left\lvert
Q_{1}\mu_{k-1}(\widehat{\gamma})(\widehat{\gamma}-\gamma)\right\rvert+\left\lvert
Q_{1}\left(\mu_{k-1}(\widehat{\gamma})-\mu_{k-1}(\gamma)\right)\right\rvert\gamma$
$\displaystyle\quad+\left\lvert\eta_{k-1}(\widehat{\gamma})Q_{1}L_{h}\left(\widehat{\gamma}-\gamma\right)\right\rvert+\left\lvert\left(\eta_{k-1}(\widehat{\gamma})-\eta_{k-1}(\gamma)\right)h(\gamma
Q_{1},U)\right\rvert+o_{p}(1)=o_{p}(1).$

$\bullet$ Bound for $\eta_{k}$. Define a constant

$\displaystyle
C_{\eta_{k}}=\left(1-\mathbb{E}\left[\frac{1}{1+\eta_{k-1}(\widehat{\gamma})g^{\prime}(d_{k-1}(\widehat{\gamma}))}\right]\right)^{-1}\left(1-\mathbb{E}\left[\frac{1}{1+\eta_{k-1}(\gamma)g^{\prime}(d_{k-1}({\gamma}))}\right]\right)^{-1}>0.$

By the triangle inequality and the technique that we use for bounding
$\eta_{1}$, we have

$\displaystyle\left\lvert\eta_{k}(\widehat{\gamma})-\eta_{k}(\gamma)\right\rvert$
$\displaystyle\leq\kappa\eta_{k-1}(\gamma)C_{\eta_{k}}\mathbb{E}\left[\left\lvert\eta_{k-1}(\widehat{\gamma})g^{\prime}\left(d_{k-1}(\widehat{\gamma})\right)-\eta_{k-1}({\gamma})g^{\prime}\left(d_{k-1}({\gamma})\right)\right\rvert\right]+O_{p}(\left\lvert\eta_{k-1}(\widehat{\gamma})-\eta_{k-1}(\gamma)\right\rvert)$
$\displaystyle\leq\kappa\eta_{k-1}^{2}(\gamma)C_{\eta_{k}}L_{g^{\prime}}\mathbb{E}[O_{p}(\left\lvert
d_{k-1}(\widehat{\gamma})-d_{k-1}(\gamma)\right\rvert)]+O_{p}(\left\lvert\eta_{k-1}(\widehat{\gamma})-\eta_{k-1}(\gamma)\right\rvert)$
$\displaystyle=o_{p}(1).$

$\bullet$ Bound for $\mu_{k}$ and $\sigma_{k}$. By the same way to show the
bound for $\mu_{1}$ and $\sigma_{1}$, we have

$\displaystyle\left\lvert\mu_{k}(\widehat{\gamma})-\mu_{k}(\gamma)\right\rvert=\left\lvert\sigma_{k}(\widehat{\gamma})-\sigma_{k}(\gamma)\right\rvert=o_{p}(1).$

Thus, we obtain

$\displaystyle\left\lvert\eta_{k}(\widehat{\gamma})-\eta_{k}(\gamma)\right\rvert=\left\lvert\mu_{k}(\widehat{\gamma})-\mu_{k}(\gamma)\right\rvert=\left\lvert\sigma_{k}(\widehat{\gamma})-\sigma_{k}(\gamma)\right\rvert=o_{p}(1),$
(39)

for any $k\in\mathbb{N}$ by induction.

Step 2. In this step, we get the conclusion from the results in Step 1. First,
we obtain

$\displaystyle\left\lvert\eta(\widehat{\gamma})-\eta(\gamma)\right\rvert=\left\lvert\mu(\widehat{\gamma})-\mu(\gamma)\right\rvert=\left\lvert\sigma(\widehat{\gamma})-\sigma(\gamma)\right\rvert=o_{p}(1).$
(40)

This follows from the fact that

$\displaystyle\left\lvert\eta(\widehat{\gamma})-\eta(\gamma)\right\rvert\leq\left\lvert\eta(\widehat{\gamma})-\eta_{k}(\widehat{\gamma})\right\rvert+\left\lvert\eta({\gamma})-\eta_{k}(\gamma)\right\rvert+o_{p}(1),$

by (39) and the first two terms on the right-hand side converge to zero for a
large $k$ limit. The results for $\mu$ and $\sigma$ also follow in the same
manner. ∎

######  Lemma C.4.

Let $\widehat{\bm{\theta}}$ be the estimator, i.e. surrogate loss minimizer,
in a GLM with a true coefficient vector $\bm{\theta}$ and features drawn
i.i.d. from $\mathcal{N}_{p}(\bm{0},\bm{I}_{p})$. Define
$(\mu_{n},\sigma_{n})$ as in (29). Then,

$\displaystyle\frac{\widehat{\bm{\theta}}-\mu_{n}\bm{\theta}}{\sigma_{n}}$

is uniformly distributed on the unit sphere lying in $\bm{\theta}^{\perp}$.

###### Proof of Lemma C.4.

Define an orthogonal projection matrix
$\bm{P}_{\bm{\theta}}=\bm{\theta}\bm{\theta}^{\top}/\left\lVert\bm{\theta}\right\rVert^{2}$
onto the line including $\bm{\theta}$, and an orthogonal projection matrix
$\bm{P}_{\bm{\theta}}^{\perp}=\bm{I}_{p}-\bm{P}_{\bm{\theta}}$ onto the
orthogonal complement of the line including $\bm{\theta}$. Let
$\bm{U}\in\mathbb{R}^{p\times p}$ be any orthogonal matrix obeying
$\bm{U}\bm{\theta}=\bm{\theta}$, i.e. any rotation operator about
$\bm{\theta}$. Then, since
$\widehat{\bm{\theta}}=\bm{P}_{\bm{\theta}}\widehat{\bm{\theta}}+\bm{P}_{\bm{\theta}}^{\perp}\widehat{\bm{\theta}}$,
we have

$\displaystyle\bm{U}\widehat{\bm{\theta}}=\bm{U}\bm{P}_{\bm{\theta}}\widehat{\bm{\theta}}+\bm{U}\bm{P}_{\bm{\theta}}^{\perp}\widehat{\bm{\theta}}=\bm{P}_{\bm{\theta}}\widehat{\bm{\theta}}+\bm{U}\bm{P}_{\bm{\theta}}^{\perp}\widehat{\bm{\theta}}.$

Using this, we obtain

$\displaystyle\frac{\bm{U}\bm{P}_{\bm{\theta}}^{\perp}\widehat{\bm{\theta}}}{\|\bm{P}_{\bm{\theta}}^{\perp}\widehat{\bm{\theta}}\|}\overset{\rm
d}{=}\frac{\bm{P}_{\bm{\theta}}^{\perp}\widehat{\bm{\theta}}}{\|\bm{P}_{\bm{\theta}}^{\perp}\widehat{\bm{\theta}}\|}=\frac{\widehat{\bm{\theta}}-\mu_{n}\bm{\theta}}{\sigma_{n}},$
(41)

where the first identity follows from the fact that
$\bm{U}\widehat{\bm{\theta}}\overset{\rm d}{=}\widehat{\bm{\theta}}$ since
$\bm{U}\widehat{\bm{\theta}}$ is the estimator with a true coefficient
$\bm{U}\bm{\theta}=\bm{\theta}$ and features drawn i.i.d. from
$\mathcal{N}(\bm{0},\bm{I}_{p})$. (41) reveals that
$(\widehat{\bm{\theta}}-\mu_{n}\bm{\theta})/\sigma_{n}$ is rotationally
invariant about $\bm{\theta}$, lies in $\bm{\theta}^{\perp}$, and has a unit
norm. These conclude the proof. ∎

######  Lemma C.5.

Suppose that
$0<\|\bm{L}^{\top}\widehat{\bm{\beta}}\|,\|\bm{L}^{\top}\bm{\beta}\|\leq C$
for some constant $C>0$. Let $\bm{\theta},\widehat{\bm{\theta}}$ be defined in
(28). If we define

$\displaystyle\mu_{n}=\frac{\widehat{\bm{\theta}}^{\top}\bm{\theta}}{\left\lVert\bm{\theta}\right\rVert^{2}},\qquad\sigma_{n}^{2}=\frac{1}{\kappa}\|\widehat{\bm{\theta}}-\mu_{n}\bm{\theta}\|^{2},$
(42)

then, under assumptions (A1) and (A2), we have

$\displaystyle\mu_{n}\overset{\mathrm{a.s.}}{\longrightarrow}\mu,\qquad\sigma_{n}^{2}\overset{\mathrm{a.s.}}{\longrightarrow}\sigma^{2},$

as $n,p(n)\to\infty$ with $p(n)/n\to\kappa$.

###### Proof of Lemma C.5.

Since $\|\bm{\theta}\|$ is bounded by an assumption,
$\bm{U}\bm{\theta}=(\left\lVert\bm{\theta}\right\rVert,\ldots,\left\lVert\bm{\theta}\right\rVert)/\sqrt{p}$
with some orthogonal matrix $\bm{U}\in\mathbb{R}^{p\times p}$ satisfies an
assumption of Lemma A.1 with $\bar{\beta}\sim\delta_{\gamma/\sqrt{p}}$ by the
fact that
$\left\lVert\bm{U}\bm{\theta}\right\rVert^{2}=\left\lVert\bm{\theta}\right\rVert^{2}=\bm{\beta}^{\top}\bm{\Sigma}\bm{\beta}=\gamma^{2}$.
Then applying Lemma A.1 to $(\bm{U}\bm{\theta},\bm{U}\widehat{\bm{\theta}})$
with $\psi(s,t)=st,t^{2}$ and considering their ratio gives

$\displaystyle\frac{\langle\bm{U}\bm{\theta},\bm{U}\widehat{\bm{\theta}}\rangle}{\left\lVert\bm{U}\bm{\theta}\right\rVert^{2}}=\mu_{n}\overset{\mathrm{a.s.}}{\longrightarrow}\mu.$

Also, setting $\psi(s,t)=(s-\mu t)^{2}$ yields

$\displaystyle\frac{1}{\kappa}\|\bm{U}\widehat{\bm{\theta}}-\mu\bm{U}\bm{\theta}\|^{2}=\frac{1}{\kappa}\|\widehat{\bm{\theta}}-\mu\bm{\theta}\|^{2}\overset{\mathrm{a.s.}}{\longrightarrow}\sigma^{2}.$

Thus, we obtain

$\displaystyle\sigma_{n}^{2}=\frac{1}{\kappa}\|\widehat{\bm{\theta}}-\mu\bm{\theta}\|^{2}+2(\mu-\mu_{n})\bm{\theta}^{\top}\widehat{\bm{\theta}}-(\mu^{2}-\mu_{n}^{2})\left\lVert\bm{\theta}\right\rVert^{2}\overset{\mathrm{a.s.}}{\longrightarrow}\sigma^{2}.$

∎

######  Lemma C.6.

Define $\left\lVert f\right\rVert_{\infty}:=\sup_{t\in\mathbb{R}}|f(t)|$ for a
function $f:\mathbb{R}\to\mathbb{R}$. A proximal operator
$\mathrm{prox}(x,f,b)\equiv\mathrm{prox}_{bF}(x)$ with $x\in\mathbb{R}$,
$b>0$, and a monotone continuous function $f=F^{\prime}$, is Lipschitz
continuous with respect to

  1. (i)

$x\in\mathbb{R}$ with constant $1$,

  2. (ii)

a monotone continuous function $f:\mathbb{R}\rightarrow\mathbb{R}$ in terms of
$L^{\infty}$ norm $\left\lVert\cdot\right\rVert_{\infty}$ with constant $b>0$,

  3. (iii)

and $b>0$ with constant $\left\lVert f\right\rVert_{\infty}$.

###### Proof of Lemma C.6.

For any $x,y\in\mathbb{R}$ and fixed monotone $f(\cdot)$, suppose that
$u=\text{prox}(x,f,b)$ and $v=\text{prox}(y,f,b)$ and $u>v$ without loss of
generality. Note that $F(\cdot)$ is convex since its derivative $f(\cdot)$ is
a monotonically increasing function. By the first-order condition, we have
$x=u+bf(u)$ and $y=v+bf(v)$. Since $f(\cdot)$ is monotonically increasing,
$0\leq(u-v)(f(u)-f(v))$. Using this,

$\displaystyle 0\leq(u-v)(bf(u)-bf(v))=(u-v)(x-u-y+v)=(u-v)(x-y)-(u-v)^{2},$

by $b>0.$ This immediately implies $|u-v|\leq|x-y|$, i.e.,
$\text{prox}(x,f,b)$ is 1-Lipschitz continuous with respect to the first
argument.

Next, set any two continuous monotone functions $f(\cdot)$ and $g(\cdot)$.
Suppose that $\text{prox}(x,f,b)\geq\text{prox}(x,g,b)$ with fixed
$x\in\mathbb{R}$ without loss of generality. Then,

$\displaystyle\left\lvert{\text{prox}(x,f,b)-\text{prox}(x,g,b)}\right\rvert$
$\displaystyle={\text{prox}(x,f,b)-\text{prox}(x,g,b)}$
$\displaystyle=b\left(g(\text{prox}(x,g,b))-f(\text{prox}(x,f,b))\right)$
$\displaystyle\leq b\left(g(\text{prox}(x,f,b))-f(\text{prox}(x,f,b))\right)$
$\displaystyle\leq b\left\lVert f-g\right\rVert_{\infty},$ (43)

where the second identity follows from the fact that
$\text{prox}(x,f)=x-bf(\text{prox}(x,f,b))$, and the first inequality is from
monotinicity of $g(\cdot)$. This means $\text{prox}(x,f,b)$ is Lipschitz
continuous for the second argument with constant $b>0$ in terms of
$\left\lVert\cdot\right\rVert_{\infty}$.

At last, using the fact that $\mathrm{prox}(x,g,b)=\mathrm{prox}(x,f,\alpha
b)$ for $g(x)=\alpha f(x),\alpha>0$ by the definition of the proximal
operator, we have, for any $b,b^{\prime}>0$,

$\displaystyle\left\lvert\mathrm{prox}(x,f,b)-\mathrm{prox}_{b^{\prime}}(x,f,b^{\prime})\right\rvert$
$\displaystyle=\left\lvert\mathrm{prox}_{b}(x,f,b)-\mathrm{prox}\left(x,\frac{b^{\prime}}{b}f,b\right)\right\rvert$
$\displaystyle\leq b\left\lVert f-\frac{b^{\prime}}{b}f\right\rVert_{\infty}$
$\displaystyle=\left\lVert bf-b^{\prime}f\right\rVert_{\infty}$
$\displaystyle=\left\lvert b-b^{\prime}\right\rvert\left\lVert
f\right\rVert_{\infty},$

where the inequality follows from (C). ∎

######  Lemma C.7.

Let $F:\mathbb{R}\to\mathbb{R}$ be a strictly convex function and
$f=F^{\prime}$. Suppose that $f(\cdot)$ takes bounded values on bounded
domains. For any bounded $x\in\mathbb{R}$ and $c,c^{\prime}>0$, we have

$\displaystyle\left\lvert\mathrm{prox}_{cF}(x)-\mathrm{prox}_{c^{\prime}F}(x)\right\rvert\leq
C\left\lvert c-c^{\prime}\right\rvert,$

for some positive constant C.

###### Proof of Lemma C.7.

By the proofs of Lemma C.6 (ii) and (iii), we can improve Lemma C.6 (iii) as

$\displaystyle\left\lvert\mathrm{prox}_{cF}(x)-\mathrm{prox}_{c^{\prime}F}(x)\right\rvert\leq\left\lvert
f(\mathrm{prox}_{cF}(x))\right\rvert\left\lvert c-c^{\prime}\right\rvert.$

Thus, we can complete the proof by showing that $\mathrm{prox}_{cF}(x)$ is
bounded. Remind that the definition of the proximal operator is

$\displaystyle\mathrm{prox}_{cF}(x)=\operatorname*{arg\,min}_{z\in\mathbb{R}}\left\\{cF(z)+\frac{1}{2}(z-x)^{2}\right\\}.$

We denote the objective function as
$H(z):=cF(z)+\frac{1}{2}(z-x)^{2}=H_{1}(z)+H_{2}(z)$. Obviously, $H_{2}$ has
the minimum value at $x$.

(i) For the case that a minimizer $\tilde{z}$ of $H_{1}$ is bounded. Note that
the minimizer is unique by the strict convexity of $H_{1}$. Without loss of
generality, suppose that $\tilde{z}\leq z$. In this case, we have
$H_{1}^{\prime}(z^{\prime})<0$ and $H_{2}^{\prime}(z^{\prime})<0$ for
$z^{\prime}<\tilde{z}$, and also have $H_{1}^{\prime}(z^{\prime})>0$ and
$H_{2}^{\prime}(z^{\prime})>0$ for $z<z^{\prime}$. Hence,
$H^{\prime}(z^{\prime})<0$ holds for $z^{\prime}<\tilde{z}$ and
$H^{\prime}(z^{\prime})>0$ holds for $z<z^{\prime}$, thus
$\mathrm{prox}_{cF}(x)\notin[-\infty,\tilde{z})\cup(z,\infty]$. In contrast,
$H^{\prime}(z^{\prime})$ may be both positive or negative for
$z^{\prime}\in[\tilde{z},z]$. Hence, $\mathrm{prox}_{cF}(x)\in[\tilde{z},z]$
holds and thus bounded itself.

(ii) For the case when $H_{1}$ has a minimum at an unbounded point, such as
$H_{1}(z)=ce^{z}$, we can also show that $\mathrm{prox}_{cF}(x)$ is bounded.
In this case, we can assume that $H_{1}$ is monotonically increasing without
loss of generality. For $z>x$, we have $H^{\prime}(z)>0$ by
$H_{1}^{\prime}(z)>0$ and $H_{2}(z)>0$. Since $H_{1}^{\prime}(z)$ is positive
and monotonically increasing by the monotonicity and the convexity of
$H_{1}(z)$, there exists a constant $\tilde{C}<x$ such that
$H_{1}^{\prime}(z)=cf(z)<-H_{2}^{\prime}(z)=-(z-x)$ for any $z<\tilde{C}$.
Putting together the results, we have $\mathrm{prox}_{cF}(x)\in[\tilde{C},x]$.
∎

## References

  * Agarwal et al. [2014] Agarwal, A., S. Kakade, N. Karampatziakis, L. Song, and G. Valiant (2014).  Least squares revisited: Scalable approaches for multi-class prediction.  In International Conference on Machine Learning, pp. 541–549. PMLR.
  * Auer et al. [1995] Auer, P., M. Herbster, and M. K. Warmuth (1995).  Exponentially many local minima for single neurons.  Advances in neural information processing systems 8.
  * Barbier et al. [2019] Barbier, J., F. Krzakala, N. Macris, L. Miolane, and L. Zdeborová (2019).  Optimal errors and phase transitions in high-dimensional generalized linear models.  Proceedings of the National Academy of Sciences 116(12), 5451–5460.
  * Bayati and Montanari [2011] Bayati, M. and A. Montanari (2011).  The dynamics of message passing on dense graphs, with applications to compressed sensing.  IEEE Transactions on Information Theory 57(2), 764–785.
  * Bellec [2022] Bellec, P. C. (2022).  Observable adjustments in single-index models for regularized m-estimators.  arXiv preprint arXiv:2204.06990.
  * Bellec and Shen [2022] Bellec, P. C. and Y. Shen (2022).  Derivatives and residual distribution of regularized m-estimators with application to adaptive tuning.  In Conference on Learning Theory, pp. 1912–1947. PMLR.
  * Bellec and Zhang [2021] Bellec, P. C. and C.-H. Zhang (2021).  Second-order stein: Sure for sure and other applications in high-dimensional inference.  The Annals of Statistics 49(4), 1864–1903.
  * Bolthausen [2014] Bolthausen, E. (2014).  An iterative construction of solutions of the tap equations for the sherrington–kirkpatrick model.  Communications in Mathematical Physics 325(1), 333–366.
  * Cai et al. [2021] Cai, T. T., Z. Guo, and R. Ma (2021).  Statistical inference for high-dimensional generalized linear models with binary outcomes.  Journal of the American Statistical Association, 1–14.
  * Cai et al. [2023] Cai, T. T., Z. Guo, and R. Ma (2023).  Statistical inference for high-dimensional generalized linear models with binary outcomes.  Journal of the American Statistical Association 118(542), 1319–1332.
  * Candès and Sur [2020] Candès, E. J. and P. Sur (2020).  The phase transition for the existence of the maximum likelihood estimate in high-dimensional logistic regression.  The Annals of Statistics 48(1), 27–42.
  * Chatterjee [2009] Chatterjee, S. (2009).  Fluctuations of eigenvalues and second order poincaré inequalities.  Probability Theory and Related Fields 143(1-2), 1–40.
  * Cordeiro [1983] Cordeiro, G. M. (1983).  Improved likelihood ratio statistics for generalized linear models.  Journal of the Royal Statistical Society: Series B (Methodological) 45(3), 404–413.
  * Detrano et al. [1989] Detrano, R., A. Janosi, W. Steinbrunn, M. Pfisterer, J.-J. Schmid, S. Sandhu, K. H. Guppy, S. Lee, and V. Froelicher (1989).  International application of a new probability algorithm for the diagnosis of coronary artery disease.  The American Journal of Cardiology 64(5), 304–310.
  * Donoho and Montanari [2016] Donoho, D. and A. Montanari (2016).  High dimensional robust m-estimation: Asymptotic variance via approximate message passing.  Probability Theory and Related Fields 166(3), 935–969.
  * Donoho et al. [2009] Donoho, D. L., A. Maleki, and A. Montanari (2009).  Message-passing algorithms for compressed sensing.  Proceedings of the National Academy of Sciences 106(45), 18914–18919.
  * Dua and Graff [2017] Dua, D. and C. Graff (2017).  UCI machine learning repository.
  * El Karoui [2018] El Karoui, N. (2018).  On the impact of predictor geometry on the performance on high-dimensional ridge-regularized generalized robust regression estimators.  Probability Theory and Related Fields 170, 95–175.
  * El Karoui et al. [2013] El Karoui, N., D. Bean, P. J. Bickel, C. Lim, and B. Yu (2013).  On robust regression with high-dimensional predictors.  Proceedings of the National Academy of Sciences 110(36), 14557–14562.
  * Fei and Li [2021] Fei, Z. and Y. Li (2021).  Estimation and inference for high dimensional generalized linear models: A splitting and smoothing approach.  Journal of Machine Learning Research 22(58), 1–32.
  * Feng et al. [2022] Feng, O. Y., R. Venkataramanan, C. Rush, and R. J. Samworth (2022).  A unifying tutorial on approximate message passing.  Foundations and Trends® in Machine Learning 15(4), 335–536.
  * Ferguson [1996] Ferguson, T. S. (1996).  A Course in Large Sample Theory.  CRC Press.
  * Guo and Chen [2016] Guo, B. and S. X. Chen (2016).  Tests for high dimensional generalized linear models.  Journal of the Royal Statistical Society. Series B (Statistical Methodology), 1079–1102.
  * Hastie et al. [2022] Hastie, T., A. Montanari, S. Rosset, and R. J. Tibshirani (2022).  Surprises in high-dimensional ridgeless least squares interpolation.  The Annals of Statistics 50(2), 949–986.
  * Janková et al. [2020] Janková, J., R. D. Shah, P. Bühlmann, and R. J. Samworth (2020).  Goodness-of-fit testing in high dimensional generalized linear models.  Journal of the Royal Statistical Society Series B: Statistical Methodology 82(3), 773–795.
  * Janková and Van De Geer [2016] Janková, J. and S. Van De Geer (2016).  Confidence regions for high-dimensional generalized linear models under sparsity.  arXiv preprint arXiv:1610.01353.
  * Jørgensen [1987] Jørgensen, B. (1987).  Exponential dispersion models.  Journal of the Royal Statistical Society: Series B (Methodological) 49(2), 127–145.
  * Kantorovich and Rubinshtein [1958] Kantorovich, L. V. and S. Rubinshtein (1958).  On a space of totally additive functions.  Vestnik of the St. Petersburg University: Mathematics 13(7), 52–59.
  * Li et al. [2023] Li, S., L. Zhang, T. T. Cai, and H. Li (2023).  Estimation and inference for high-dimensional generalized linear models with knowledge transfer.  Journal of the American Statistical Association, 1–12.
  * Loureiro et al. [2021] Loureiro, B., C. Gerbelot, H. Cui, S. Goldt, F. Krzakala, M. Mézard, and L. Zdeborová (2021).  Capturing the learning curves of generic features maps for realistic data sets with a teacher-student model.  arXiv preprint arXiv:2102.08127.
  * Mai and Liao [2019] Mai, X. and Z. Liao (2019).  High dimensional classification via regularized and unregularized empirical risk minimization: Precise error and optimal loss.  arXiv preprint arXiv:1905.13742.
  * Mai et al. [2019] Mai, X., Z. Liao, and R. Couillet (2019).  A large scale analysis of logistic regression: Asymptotic performance and new insights.  In ICASSP 2019-2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp. 3357–3361. IEEE.
  * McCullagh [1980] McCullagh, P. (1980).  Regression models for ordinal data.  Journal of the Royal Statistical Society: Series B (Methodological) 42(2), 109–127.
  * Mei and Montanari [2022] Mei, S. and A. Montanari (2022).  The generalization error of random features regression: Precise asymptotics and the double descent curve.  Communications on Pure and Applied Mathematics 75(4), 667–766.
  * Miolane and Montanari [2021] Miolane, L. and A. Montanari (2021).  The distribution of the lasso: Uniform control over sparse balls and adaptive parameter tuning.  The Annals of Statistics 49(4), 2313–2335.
  * Montanari and Saeed [2022] Montanari, A. and B. N. Saeed (2022).  Universality of empirical risk minimization.  In Conference on Learning Theory, pp. 4310–4312. PMLR.
  * Rangan [2011] Rangan, S. (2011).  Generalized approximate message passing for estimation with random linear mixing.  In 2011 IEEE International Symposium on Information Theory Proceedings, pp. 2168–2172. IEEE.
  * Rangan et al. [2016] Rangan, S., P. Schniter, E. Riegler, A. K. Fletcher, and V. Cevher (2016).  Fixed points of generalized approximate message passing with arbitrary matrices.  IEEE Transactions on Information Theory 62(12), 7464–7474.
  * Salehi et al. [2019] Salehi, F., E. Abbasi, and B. Hassibi (2019).  The impact of regularization on high-dimensional logistic regression.  Advances in Neural Information Processing Systems 32.
  * Seddik et al. [2021] Seddik, M. E. A., C. Louart, R. Couillet, and M. Tamaazousti (2021).  The unexpected deterministic and universal behavior of large softmax classifiers.  In International Conference on Artificial Intelligence and Statistics, pp. 1045–1053. PMLR.
  * Sur and Candès [2019] Sur, P. and E. J. Candès (2019).  A modern maximum-likelihood theory for high-dimensional logistic regression.  Proceedings of the National Academy of Sciences 116(29), 14516–14525.
  * Sur et al. [2019] Sur, P., Y. Chen, and E. J. Candès (2019).  The likelihood ratio test in high-dimensional logistic regression is asymptotically a rescaled chi-square.  Probability theory and related fields 175(1), 487–558.
  * Tan and Bellec [2023] Tan, K. and P. C. Bellec (2023).  Multinomial logistic regression: Asymptotic normality on null covariates in high-dimensions.  arXiv preprint arXiv:2305.17825.
  * Thrampoulidis et al. [2015] Thrampoulidis, C., E. Abbasi, and B. Hassibi (2015).  Lasso with non-linear measurements is equivalent to one with linear measurements.  Advances in Neural Information Processing Systems 28.
  * Thrampoulidis et al. [2018] Thrampoulidis, C., E. Abbasi, and B. Hassibi (2018).  Precise error analysis of regularized $m$-estimators in high dimensions.  IEEE Transactions on Information Theory 64(8), 5592–5628.
  * Tian and Taylor [2017] Tian, X. and J. Taylor (2017).  Asymptotics of selective inference.  Scandinavian Journal of Statistics 44(2), 480–499.
  * Van de Geer et al. [2014] Van de Geer, S., P. Bühlmann, Y. Ritov, and R. Dezeure (2014).  On asymptotically optimal confidence regions and tests for high-dimensional models.  The Annals of Statistics, 1166–1202.
  * van der Vaart [2000] van der Vaart, A. W. (2000).  Asymptotic statistics, Volume 3.  Cambridge University Press.
  * Venkataramanan et al. [2022] Venkataramanan, R., K. Kögler, and M. Mondelli (2022).  Estimation in rotationally invariant generalized linear models via approximate message passing.  In International Conference on Machine Learning, pp. 22120–22144. PMLR.
  * Yadlowsky et al. [2021] Yadlowsky, S., T. Yun, C. Y. McLean, and A. D’Amour (2021).  Sloe: A faster method for statistical inference in high-dimensional logistic regression.  Advances in Neural Information Processing Systems 34, 29517–29528.
  * Zhao et al. [2022] Zhao, Q., P. Sur, and E. J. Candes (2022).  The asymptotic distribution of the mle in high-dimensional logistic models: Arbitrary covariance.  Bernoulli 28(3), 1835–1861.