1 Introduction↩︎

This paper studies the performance of bagging estimators trained on subsampled, overlapping datasets in the context robust linear regression and logistic regression.

1.1 M-estimation in the proportional regime↩︎

We consider an M-estimation problem in the proportional regime where sample size \(n\) and dimension \(p\) are of the same order: Throughout the paper, \(\delta > 1\) is a fixed constant and the ratio \[n/p = \delta \label{regime}\tag{1}\] is held fixed as \(n,p\to+\infty\) simultaneously. The practitioner collects data \((y_i,\boldsymbol{x}_i)_{i\in[n]}\) with scalar-valued responses \(y_i\) and feature vectors \(\boldsymbol{x}_i \in \mathbb{R}^p\). For a given subset of observations \(I\subset [n]\), an estimator \(\hat{\boldsymbol{b}}(I)\) is trained on the subset of observations \((y_i,\boldsymbol{x}_i)_{i\in I}\) using an optimization problem of the form \[\hat{\boldsymbol{b}}(I) = \mathop{\mathrm{arg\,min}}_{\boldsymbol{b}\in \mathbb{R}^p} \sum_{i\in I} \ell_i(\boldsymbol{x}_i^\top \boldsymbol{b}) \label{hbbI95intro}\tag{2}\] where for each \(i\in[n]\), the loss \(\ell_i(\cdot)\) is convex and depends implicitly on the response \(y_i\). We will focus on two regression settings: robust linear regression and Generalized Linear Models (GLM), including logistic regression. In robust regression, the response is of the form \[\label{y95i95robust} y_i = \boldsymbol{x}_i^T \boldsymbol{\beta}^* + \varepsilon_i\tag{3}\] for some possibly heavy-tailed noise \(\varepsilon_i\) independent of \(\boldsymbol{x}_i\). In this case the loss \(\ell_i\) in 2 is given by \[\ell_i(u) = \rho(y_i-u) \label{ell95i95robust}\tag{4}\] where \(\rho\) is a deterministic function, for instance the Huber loss \(\rho (u) = \int_0^{|u|}\min(1,t)dt\) or its smooth variants, e.g., \(\rho(u)=\sqrt{1+u^2}\). The asymptotics of the performance of 2 with \(I=\{1,...,n\}\) and the loss 4 in robust regression in the proportional regime 1 are now well understood [1]–[4] as we will review in 2. A typical example of GLM to which our results apply is the case of binary logistic regression, where \(\ell_i\) in 2 is the negative log-likelihood \[\ell_i(u) = \log(1+e^u) - u y_i, \qquad y_i\in \{0,1\} \label{logi95loss}\tag{5}\] which is now also well understood for \(I=[n]\) in 2. Related results will be reviewed in 3. The goal of the present paper is to study the performance of bagging several estimators of the form 2 obtained from several subsampled datasets \(I_1,...,I_M\).

1.2 Bagging estimators trained on subsampled datasets without replacement↩︎

Let \(M>0\) be a fixed integer, held fixed as \(n,p\to+\infty\). The practitioner then samples \(M\) subsets of \([n]\) according to the uniform distribution on all subsets of \([n]\) of size \(qn\) for some \(q\in (0,1]\), that is, \[I_1,...,I_M \sim^{\text{iid}} \text{Unif}\{I\subset [n] : |I|=qn \}. \label{I95195I95M}\tag{6}\] Each \(I_m\) thus samples a subset of \([n]\) of size \(qn\) without replacement, and the set of indices \(I_1,...,I_M\) are all independent. While the set of indices are independent, the corresponding subsampled datasets \[(\boldsymbol{x}_i, y_i)_{i \in I_m}, \quad \text{ and } \quad (\boldsymbol{x}_i, y_i)_{i \in I_{m'}} \label{two95datasets}\tag{7}\] are not independent as soon as there is some overlap in the sense \(I_m\cap I_{m'} \ne \emptyset\).

The goal of the paper is to understand the performance of bagging the corresponding subsampled estimates: with the notation \(\hat{\boldsymbol{b}}(I)\) in 2 and \(I_1,...,I_M\) in 6 , the practitioner construct the bagged estimate \[\label{bagging95intro} \bar\boldsymbol{b}= \frac{1}{M} \sum_{m=1}^M \hat{\boldsymbol{b}}(I_m).\tag{8}\]

Related works↩︎

Bagging as a generally applicable principle was introduced in [5], [6]. Early analysis of bagging in low-dimensional regimes were performed in [7] among others. In the proportional regime 1 , [8] demonstrated the role of bagging as an implicit regularization technique when the base learners \(\hat{\boldsymbol{b}}(I_m)\) are least-squares estimates. Bagging Ridge estimators was studied in [9], [10] who characterized the limit of the squared error of 8 using random matrix theory. The implicit regularization power of bagging in the proportional regime is again seen in [9], [10], where it is shown that the optimal risk among Ridge estimates can also be achieved by bagging Ridgeless estimates and optimally choosing the sub sample size. Estimating the risk of a bagged estimate such as 8 for regularized least-squares estimates is done in [9]–[11]. The risk of bagging random-features estimators, trained on the full dataset but with each base learner having independent weights within the random feature activations, is characterized in [12]. Most recently, [13] studied the limiting equations of several resampling schemes including bootstrap and resampling without replacement, and characterized self-consistent equations for the limiting risk of estimators obtained by minimization of the negative log-likelihood and an additive Ridge penalty.

Organization↩︎

We will first study and state our main results for robust regression in 2. 3 extends the results to logistic regression. Numerical simulations are provided in 2.5 in robust regression and in 3.3 in logistic regression. The main results are proved in 4 simultaneously for robust linear regression and logistic regression. 5 contains several auxiliary lemmas used in the proof in 4.

Notation↩︎

For vectors \(\|\cdot\|\) or \(\|\cdot\|_2\) is the Euclidean norm, while \(\|\cdot\|_{op}\) and \(\|\cdot\|_F\) denote the operator norm and Frobenius norm of matrices. The arrow \(\to^P\) denotes convergence in probability and \(o_P(1)\) denotes any sequence of random variables converging to 0 in probability. The stochastically bounded notation \(O_P(r_n)\) for \(r_n>0\) denotes a sequence of random variables such that for any \(\eta>0\), there exists \(K>0\) with \(\mathop{\mathrm{\mathbb{P}}}(O_P(r_n)>K r_n)\le \eta\).

2 Robust regression↩︎

This section focuses on robust regression in the linear model 3 , where the noise variables \(\varepsilon_i\) are possibly heavy-tailed. Throughout the paper, our working assumption for the robust linear regression setting is the following.

2.1 A review of existing results in robust linear regression↩︎

The seminal works [1]–[3], [14] characterized the performance of robust M-estimation in the proportional regime 1 . For a convex loss \(\rho:\mathbb{R}\to\mathbb{R}\) ans \(\ell_i\) as in 1, these works characterized the limiting squared risk \(\|\hat{\boldsymbol{b}}(\{1,...,n\})-\boldsymbol{\beta}^*\|^2\) of an estimator \(\hat{\boldsymbol{b}}(\{1,...,n\})\), trained on the full dataset, i.e., taking \(I=\{1,...,n\}\) in 2 . In particular, [1]–[4], [14] show that under the design of \((\boldsymbol{x}_i,y_i)\) given in 1, the squared risk of \(\hat{\boldsymbol{b}}(\{1,...,n\})\) converges in probability to a constant, and this constant is found by solving a system of two nonlinear equations with two unknowns. If a subset \(I\subset [n]\) of size \(|I|=qn\) is used to train 2 , simply changing \(\delta=n/p\) to \(\delta q=|I|/p\), these results imply the convergence in probability \(\|\hat{\boldsymbol{b}}(I)-\boldsymbol{\beta}^*\|^2\to^P \sigma^2\) where \((\sigma,\gamma)\) is the solution to the system \[\begin{align} \tag{9} \tfrac{\sigma^2}{\delta q}&= \mathop{\mathrm{\mathbb{E}}}[ (\sigma G_i - \mathop{\mathrm{prox}}[\gamma\ell_i](\sigma G_i))^2 ] \\ 1 - \tfrac{1}{\delta q} &= \sigma^{-1} \mathop{\mathrm{\mathbb{E}}}[G_i \mathop{\mathrm{prox}}[\gamma\ell_i](\sigma G_i) ] \tag{10} \end{align}\] where \(G_i\sim N(0,1)\) is independent of \(\varepsilon_i\). Above, \(\mathop{\mathrm{prox}}[f](x_0)=\mathop{\mathrm{arg\,min}}_{x\in\mathbb{R}}(x_0-x)^2/2 + f(x)\) denotes the proximal operator of a convex function \(f\) for any \(x_0\in \mathbb{R}\). The system 9 10 was predicted in [1] using a heuristic leave-one-out argument. Early rigorous results [2], [3], [14] assumed either \(\rho\) strongly convex [2] or added an additive strongly-convex Ridge penalty to the M-estimation problem [3], [14]; [4] generalized such results without strong convexity.

We now subsample without replacement, obtaining iid subsets \(I_1,...,I_M\) as in 6 . For each \(m=1,...,M\) the theory above applies individually to \(\hat{\boldsymbol{b}}(I_m)\). In particular \(\|\hat{\boldsymbol{b}}(I_m)-\boldsymbol{\beta}^*\|^2\to^P \sigma^2\). By expanding the square, the squared L2 error of the average \(\bar\boldsymbol{b}\) in 8 is given by \[\|\bar\boldsymbol{b}-\boldsymbol{\beta}^*\|^2 = \frac{1}{M^2} \sum_{m=1}^M \|\hat{\boldsymbol{b}}(I_m)-\boldsymbol{\beta}^*\|^2 + \frac{1}{M^2} \sum_{m=1}^M \sum_{m'=1: m'\ne m}^M (\hat{\boldsymbol{b}}(I_m)-\boldsymbol{\beta}^*)^T (\hat{\boldsymbol{b}}(I_{m'})-\boldsymbol{\beta}^*). \label{risk95bagged}\tag{11}\] Since previous works established that \(\|\hat{\boldsymbol{b}}(I_m)-\boldsymbol{\beta}^*\|^2\to^P \sigma^2\), the first term above is clearly \(\sigma^2/M\). The crux of the problem is thus to characterize the limit in probability, if any, of each term \((\hat{\boldsymbol{b}}(I_m)-\boldsymbol{\beta}^*)^T (\hat{\boldsymbol{b}}(I_{m'})-\boldsymbol{\beta}^*)\) in the second term inside the double sum.

2.2 A glance at our results↩︎

Since \(\rho\) in 4 is Lipschitz and differentiable, the system 9 10 admits a unique solution [15]. Let \((\sigma,\gamma )\) be the solution¹ to this system.

The key to understanding the performance of the aforementioned bagging procedure 8 and, for instance, characterizing the limits of \(\| \bar\boldsymbol{b}- \boldsymbol{\beta}^*\|^2\), is the following equation with unknown \(\eta\in[-1,1]\): \[\begin{align} \eta&=\frac{{q^2}\delta}{\sigma^2} \mathop{\mathrm{\mathbb{E}}}\Bigl[ \Bigl(\sigma G_i - \mathop{\mathrm{prox}}[\gamma \ell_i](\sigma G_i)\Bigr) \Bigl(\sigma G_i - \mathop{\mathrm{prox}}[\gamma \ell_i](\sigma \tilde{G}_i)\Bigr) \Bigr] \tag{12} \\&\qquad \text{ where } \begin{pmatrix} G_i \\ \tilde{G}_i \end{pmatrix} \sim N\Bigl(0_2, \begin{pmatrix} 1 & \eta \\ \eta & 1 \end{pmatrix} \Bigr) \tag{13} \end{align}\] with \((G_i,\tilde{G}_i)\) independent of \((\ell_i,\varepsilon_i)\), and \((\ell_i,\varepsilon_i)\) are the same as in 3 4 with the noise \(\varepsilon_i\) independent of \(\boldsymbol{x}_i\). Using 9 , the above equation can be equivalently rewritten as \[\eta = F(\eta) \quad \text{where } F(\eta) \mathrel{\vcenter{:}}= q\frac{ \mathop{\mathrm{\mathbb{E}}}\bigl[ \bigl(\sigma G_i - \mathop{\mathrm{prox}}[\gamma \ell_i](\sigma G_i)\bigr) \bigl(\sigma G_i - \mathop{\mathrm{prox}}[\gamma \ell_i](\sigma \tilde{G}_i)\bigr) \bigr] }{ \mathop{\mathrm{\mathbb{E}}}[(\sigma G_i- \mathop{\mathrm{prox}}[\gamma\ell_i](\sigma G_i))^2 ] } \label{def-F}\tag{14}\] since \(\mathop{\mathrm{\mathbb{E}}}[(\sigma G_i- \mathop{\mathrm{prox}}[\gamma\ell_i](\sigma G_i))^2 ] = \sigma^2/(\delta q)\) in the denominator by 9 . This shows that any solution \(\eta\) must satisfy \(|\eta|\le q\) by the Cauchy-Schwarz inequality.

We will show in the next section that this equation in \(\eta\) has a unique solution. Our main results imply a close relationship between the solution \(\eta\) of 12 and the bagged estimates, in particular 30 satisfies \[\bigl(\hat{\boldsymbol{b}}(I_m)-\boldsymbol{\beta}^*\bigr)^T \bigl(\hat{\boldsymbol{b}}(I_{m'}) - \boldsymbol{\beta}^*\bigr) \to^P \eta \sigma^2.\] For two distinct and fixed \(m < m'\), the solution \(\eta\) further characterizes the joint distribution of two predicted values \(\boldsymbol{x}_i^T\hat{\boldsymbol{b}}(I_m)\) and \(\boldsymbol{x}_i^T\hat{\boldsymbol{b}}(I_m)\) with \(i\in I_m\cap I_{m'}\), by showing the existence of \((G_i,\tilde{G}_i)\) as in 13 , independent of \((\ell_i,U_i)\) and such that \[\begin{align} \boldsymbol{x}_i^T\hat{\boldsymbol{b}}(I_m) &=\mathop{\mathrm{prox}}[\gamma\ell_i](\sigma G_i) + o_P(1) , \\ \boldsymbol{x}_i^T\hat{\boldsymbol{b}}(I_{m'}) &= \mathop{\mathrm{prox}}[\gamma\ell_i](\sigma\tilde{G}_i) + o_P(1). \end{align}\]

2.3 Existence and uniqueness of solutions to the fixed-point equation↩︎

2.4 Main results in robust regression↩︎

For any \(I\subset[n]\) with \(|I|=qn = q\delta p\), the M-estimator \(\hat{\boldsymbol{b}}(I)=\mathop{\mathrm{arg\,min}}_{\boldsymbol{b}\in\mathbb{R}^p}\sum_{i\in I} \ell_i(\boldsymbol{x}_i^T\boldsymbol{b})\) satisfies the convergence in probability \[\|\hat{\boldsymbol{b}}(I)-\boldsymbol{\beta}^*\|^2 \to^P \sigma^2, \qquad \frac{1}{|I|}\sum_{i\in I} \Bigl(\ell_i'(\boldsymbol{x}_i^T\hat{\boldsymbol{b}}(I))\Bigr)^2 \to^P \frac{\sigma^2}{\gamma^2 q\delta}. \label{eq95convergence}\tag{17}\] The first convergence in probability was proved by many authors, e.g., [1]–[4]. The second can be obtained using the CGMT of [4], see for instance [16]. We will take the convergence in probability 17 for granted in our proof.

3 is proved in 4. It provides three messages. First, 18 states that the correlation \((\hat{\boldsymbol{b}}(I)-\boldsymbol{\beta}^*)^T(\hat{\boldsymbol{b}}(\tilde{I})-\boldsymbol{\beta}^*)\) between two estimators trained in independent subsets \(I,\tilde{I}\) both of cardinally \(qn\) converges to the unique solution \(\eta\) of 12 . A direct consequence is that the squared risk of the bagged estimate 11 satisfies \[\|\bar\boldsymbol{b}-\boldsymbol{\beta}^*\|^2 \to^P \sigma^2/M + (1-1/M)\sigma^2\eta. \label{risk95bagged95limit}\tag{21}\] Second, both terms in this risk decomposition of the bagged estimate \(\bar\boldsymbol{b}\) can be estimated using 19 averaged over all pairs \((I_m,I_{m'})_{m\ne m'}\), that is, \[\frac{1}{M^2} \sum_{m\ne m'} \frac{\hat{\gamma}(I_m) \hat{\gamma}(\tilde{I}_{m'})}{p} \sum_{i\in I_{m'} \cap \tilde{I}_m}\ell_i'\Bigl(\boldsymbol{x}_i^T\hat{\boldsymbol{b}}(I_{m'})\Bigr) \ell_i'\Bigl(\boldsymbol{x}_i^T\hat{\boldsymbol{b}}(\tilde{I})\Bigr) \to^P \Bigl(1-\frac{1}{M}\Bigr)\eta\sigma^2,\] and \(\frac{1}{M^2}\sum_{m=1}^M\frac{\hat{\gamma}(I_m)^2}{p} \sum_{i\in I_m} \ell_i'(\boldsymbol{x}_i^T\hat{\boldsymbol{b}}(I_m))^2 \to^P \sigma^2/M\). These estimators let us estimate the risk of the bagged estimate 21 , for instance to choose an optimal subsample size \(q\in(0,1)\), or to choose a large enough constant \(M>0\) so that 21 is close to the large-\(M\) limit given by \(\sigma^2\eta\). These estimators are reminiscent of the Corrected Generalized Cross-Validation developed for regularized least-squares estimators in [11].

As shown in 1, resampling and bagging is sometimes beneficial but not always. Whether the curve \(q\mapsto \sigma^2\eta\) is U-shaped and minimized at some \(q^*<1\) (i.e., bagging is beneficial) depends on the interplay between the oversampling ratio \(\delta=n/p\), the distribution of the noise \(\varepsilon_i\) and the robust loss function \(\rho\) used in 2 . In 1, we observe that if \(\varepsilon_i/\tau\) has t-distribution with 2 degrees of freedom and \(\delta=5\), subsampling is not beneficial for \(\tau=1\) but becomes beneficial for \(\tau\ge 1.5\). The generality of this phenomenon is unclear at this point.

The third message of 3 is the characterization of the limiting bivariate distribution of \((\boldsymbol{x}_i^T\hat{\boldsymbol{b}}(I),\boldsymbol{x}_i^T\hat{\boldsymbol{b}}(\tilde{I}))\) for an observation \(i\in I \cap \tilde{I}\) used to train both \(\hat{\boldsymbol{b}}(I)\) and \(\hat{\boldsymbol{b}}(\tilde{I})\). The convergence 20 implies that \((\boldsymbol{x}_i^T\hat{\boldsymbol{b}}(I),\boldsymbol{x}_i^T\hat{\boldsymbol{b}}(\tilde{I}))\) converges weakly to the distribution of \((\mathop{\mathrm{prox}}[\gamma\ell_i](\sigma G_i),\mathop{\mathrm{prox}}[\gamma\ell_i](\sigma \tilde{G}_i))\) where \((G_i,\tilde{G}_i)\) has distribution 13 , as in the fixed-point equation 12 satisfied by \(\eta\).

The setting of resampling without replacement in the proportional regime of the present paper is also studied in the recent paper [13]. There are some significant differences between our contributions and [13]. First, an additive Ridge penalty is imposed in [13] and multiple resampling schemes are studied, while our object of interest is the unregularized M-estimator 2 with a focus on resampling without replacement. The simple fixed-point equation 18 does not appear explicitly in [13], which instead focuses on self-consistent equations satisfied by bias and variance functionals [13] of the specific resampling scheme under study. Another distinctive contribution of the present paper is the proposed estimator 19 which can be used to optimally tune the subsample size, and the proof that the equation 12 admits a unique solution. The use of an additive Ridge penalty brings strong convexity to the optimization problem and simplifies the analysis, as observed in [1]; in this case this makes the analysis [12] based on [17] readily applicable.

2.5 Numerical simulations in robust regression↩︎

Let us verify 3 with numerical simulations. Throughout this section, we focus on the Huber loss \[\rho(t)= \begin{cases} t^2/2 &\text{ if }|t|<1,\\ |t|-1/2& \text{ if }|t|\ge 1. \end{cases} \label{huber}\tag{22}\] The oversampling ratio \(\delta=n/p\) is fixed to \(5\). First, we plot \(\eta\) and \(\sigma^2\eta\) as functions of \(\eta\in [1/\delta, 1]\) for different noise scales: we change the noise distribution as \(\{\text{scale}\}\times \text{t-dist (df=2)}\), \(\text{scale}\in \{1, 1.5, 2, 5, 10\}\). The left figures in 1 imply that the curve \(q\mapsto \eta\) is nonlinear. Note that the dashed line is the affine line \(q\mapsto (q-\delta^{-1})/(1-\delta^{-1})\). More interestingly, the larger the noise scale is, the larger the nonlinearity is. In the right figures in 1, we observe that the plot \(\eta\mapsto\eta\alpha^2\) takes a U-shape curve when the noise scale is sufficiently large. This simulation result suggests that as the scale of noise distribution increases, sub-sampling is eventually beneficial in the sense that the limit of 21 as \(M\to+\infty\) is smaller than the squared error of a single estimate trained on the full dataset. Next, we compare in simulations the correlation \[\frac{(\hat{\boldsymbol{b}}(I)-\boldsymbol{\beta}^*)^T(\hat{\boldsymbol{b}}(\tilde{I})-\boldsymbol{\beta}^*)}{\|(\hat{\boldsymbol{b}}(I)-\boldsymbol{\beta}^*)\|_2 \|(\hat{\boldsymbol{b}}(\tilde{I})-\boldsymbol{\beta}^*)\|_2}\] and the inner product \((\hat{\boldsymbol{b}}(I)-\boldsymbol{\beta}^*)^T(\hat{\boldsymbol{b}}(\tilde{I})-\boldsymbol{\beta}^*)\) with their theoretical limits \((\eta, \eta\sigma^2)\) as in 18 , as well as the estimator in 19 . Here, the noise distribution is fixed to \(3\cdot \text{t-dist(df=2)}\) with \((n, p)=(5000, 1000)\) and \({100}\) repetitions. 2 implies that the correlation and product are approximated well by the corresponding theoretical values and estimates.

3 Resampling without replacement in logistic regression↩︎

3.1 A review of existing results in logistic regression↩︎

Let \(\nu >0, q\in(0,1], \delta>1\) be fixed constants. If a single estimator \(\hat{\boldsymbol{b}}(I)\) is trained with 2 on a subset of observations \(I\subset [n]\) with \(|I|/n=q\) for some constant \(q\in(0,1]\) held fixed as \(n,p\to+\infty\), the behavior of \(\hat{\boldsymbol{b}}(I)\) is now well-understood when \((y_i,\boldsymbol{x}_i)_{i\in [n]}\) are iid with \(\boldsymbol{x}_i \sim N(\boldsymbol{0}_p,\boldsymbol{I}_p)\) normally distributed and the conditional distribution \(y_i \mid \boldsymbol{x}_i\) following a logistic model of the form \[\mathop{\mathrm{\mathbb{P}}}\Bigl(y_i= 1 \mid \boldsymbol{x}_i\Bigr) = \frac{1}{1+\exp(-\boldsymbol{x}_i^T\boldsymbol{\beta}^*)} = \frac{1}{1+\exp(- \nu \boldsymbol{x}_i^T\boldsymbol{w})}\] where \(\boldsymbol{\beta}^*\) is a ground truth with \(\|\boldsymbol{\beta}^*\|=\nu\), and \(\boldsymbol{w}= \boldsymbol{\beta}^*/\nu\) is the projection of \(\boldsymbol{\beta}^*\) on the unit sphere. In this logistic regression model, the limiting behavior of \(\hat{\boldsymbol{b}}(I)\) with the logistic loss 5 trained using \(|I|=(\delta q) p\) samples is characterized as follows: there exists a monotone continuous function \(h(\cdot)\) (with explicit expression given in [19]) such that:

\(\bullet\) If \(\delta q< h(\nu)\) then the logistic MLE 2 does not exist with high-probability.

\(\bullet\) If \(\delta q> h(\nu)\) then there exists a unique [20] solution \((\sigma_*,a_*,\gamma_*)\) to the following the low-dimensional system of equations \[\begin{align} \tag{23} \tfrac{\sigma^2}{\delta q}&= \mathop{\mathrm{\mathbb{E}}}[ (a U_i + \sigma G_i - \mathop{\mathrm{prox}}[\gamma\ell_i](a U_i + \sigma G_i))^2 ], \\ 0 &= \mathop{\mathrm{\mathbb{E}}}[ (a U_i + \sigma G_i - \mathop{\mathrm{prox}}[\gamma\ell_i](a U_i + \sigma G_i)) ], \tag{24} \\ 1 - \tfrac{1}{\delta q} &= \sigma^{-1} \mathop{\mathrm{\mathbb{E}}}[G_i \mathop{\mathrm{prox}}[\gamma\ell_i]'(a U_i + \sigma G_i) ] \tag{25} \end{align}\] where \(U_i = \boldsymbol{x}_i^T\boldsymbol{w}\) and \(G_i\sim N(0,1)\) is independent of \((\ell_i,U_i)\). Above, \(\mathop{\mathrm{prox}}[f](x_0)=\mathop{\mathrm{arg\,min}}_{x\in\mathbb{R}}(x_0-x)^2/2 + f(x)\) denotes the proximal operator of any convex function \(f\) for any \(x_0\in \mathbb{R}\). In this region \(\{\delta q> h(\nu)\}\) where the above system admits a unique solution \((a,\sigma,\gamma)\), the logistic MLE 2 exists with high-probability and the following convergence in probability holds, \[\begin{align} \boldsymbol{w}^T \hat{\boldsymbol{b}}(I) &\to^P a , \tag{26} \\ \qquad \|(\boldsymbol{I}_p-\boldsymbol{w}\boldsymbol{w}^T)\hat{\boldsymbol{b}}(I)\|^2 &\to^P \sigma^2, \tag{27} \\ \frac{1}{|I|}\sum_{i\in I}\ell_i'\Bigl(\boldsymbol{x}_i^T\hat{\boldsymbol{b}}(I)\Bigr)^2 &\to^P \frac{\sigma^2}{ \gamma^2q\delta}. \tag{28} \end{align}\] by [20], [21] for the first two lines and [16] for the third. Further results are obtained in [19], [20], [22], including asymptotic normality results for individual components \(\hat{b}_j\) of 2 . Note that the 3-unknowns system 23 25 is stated in these existing works after integration of the distribution of \(y_i\). We choose the equivalent formulation 23 25 without integrating the conditional distribution of \(y_i\) as the form 23 25 is closer to 9 10 from robust regression, and closer and to the quantities naturally appearing in our proofs. In 4, this common notation is useful to prove the main results simultaneously for robust linear regression and logistic regression.

While the limit in probability of the correlation \(\bar\boldsymbol{b}^T\boldsymbol{\beta}^*\) can be deduced directly from 26 , the case of Mean Squared Error (MSE) \(\|\bar\boldsymbol{b}- \boldsymbol{\beta}^*\|^2\) or the correlation \(\bar\boldsymbol{b}^T\boldsymbol{\beta}^*\) is more subtle. To see the crux of the problem, recall \(\boldsymbol{w}= \boldsymbol{\beta}^*/\|\boldsymbol{\beta}^*\|\), define \(\boldsymbol{P}= (\boldsymbol{I}_p-\boldsymbol{w}\boldsymbol{w}^T)\) for brevity, and consider the decomposition \[\label{MSE} \|\bar\boldsymbol{b}- \boldsymbol{\beta}^*\|^2 = \bigl(\boldsymbol{w}^T(\bar\boldsymbol{b}- \boldsymbol{\beta}^*)\bigr)^2 + \|\boldsymbol{P}\bar \boldsymbol{b}\|^2\tag{29}\] where the second term is \[\label{intro95inner95product95P} \bigl(\boldsymbol{w}^T(\bar\boldsymbol{b}- \boldsymbol{\beta}^*)\bigr)^2 + \frac{1}{M^2}\sum_{m,m'=1}^M \hat{\boldsymbol{b}}(I_m)^T\boldsymbol{P}\hat{\boldsymbol{b}}(I_{m'}).\tag{30}\] In order to characterize the limit of the MSE of \(\bar\boldsymbol{b}\), or to characterize the limit of the normalized correlation \(\|\bar\boldsymbol{b}\|^{-1} \bar\boldsymbol{b}^T\boldsymbol{\beta}^*\), we need to first understand the limit of the inner product \[\hat{\boldsymbol{b}}(I_m)^T\boldsymbol{P}\hat{\boldsymbol{b}}(I_{m'}) \label{corel95logi}\tag{31}\] where \(\hat{\boldsymbol{b}}(I_m)\) and \(\hat{\boldsymbol{b}}(I_{m'})\) are trained on two subsamples \(I_m\) and \(I_{m'}\) with non-empty intersection. This problem happens to be almost equivalent to the corresponding one in robust regression, and we will prove the following result and 3 simultaneously.

3.2 Main results for logistic regression↩︎

In other words, we assume a logistic model with parameters on the side of the phase transition where the MLE exists with high-probability. In this regime, the system 23 25 admits a unique solution \((a,\sigma,\gamma)\) and the convergence in probability 26 28 holds.

We omit the proof since it is exactly same as the proof of 2. Similarly to robust regression in 3 the solution \(\eta\) to ?? characterizes the limit in probability of the correlation 31 , the estimator 19 is still valid for estimating \(\eta\sigma^2\), and finally we can characterize the joint distribution of two predicted values \(\boldsymbol{x}_i^T\hat{\boldsymbol{b}}(I_m)\) and \(\boldsymbol{x}_i^T\hat{\boldsymbol{b}}(I_m)\) for an observation \(i\in I_m\cap I_{m'}\) appearing in both datasets.

3.3 Numerical simulations in logistic regression↩︎

Similarly to 2.5, we check the accuracy of 6 with numerical simulations. Here, \((n, p)\) is fixed to \((5000, 500)\) so that \(\delta=n/p=10\). For each signal strength \(\|\boldsymbol{\beta}_*\|\in\{1,2\}\), we compute the correlation \[\hat{\boldsymbol{b}}(I)^T \boldsymbol{P}\hat{\boldsymbol{b}}(\tilde{I})/\bigl(\|\boldsymbol{P}\hat{\boldsymbol{b}}(I)\|_2\|\boldsymbol{P}\hat{\boldsymbol{b}}(\tilde{I})\|_2\bigr)\] and the inner product \(\hat{\boldsymbol{b}}(I)^T \boldsymbol{P}\hat{\boldsymbol{b}}(\tilde{I})\) as we change the sub-sampling ratio \(q=k/n\in[0.4, 1]\) and the estimate constructed by 19 . We perform \(100\) repetitions. The theoretical limits \((\eta, \sigma^2\eta)\) are obtained by solving ?? numerically. 3 shows that the theoretical curves (\(q\mapsto \eta\) and \(q\mapsto \sigma^2\eta\)) match with the correlation and the inner product. The estimator 19 is accurate for medium to large subsample ratio \(q\), but appears slightly biased upwards for small values of \(q\). The source of this slight upward bias is unclear, although possibly due to the finite-sample nature of the simulations \((p=500)\).

In all simulations for logistic regression that we have performed, the curve \(q\mapsto \eta\) is affine, as in the left plot in 3. The reason for this is unclear to us at this point and this appears to be specific logistic regression; for instance the curve \(q\mapsto\eta\) in 1 for robust regression are clearly non-affine.

4 Proof of the main results↩︎

We prove here [thm:logi,thm:robust] simultaneously using the following notation:

• In Robust regression (3), set \(a=0\), let \((\sigma,\gamma)\) be the unique solution to 9 10 , let \(\boldsymbol{\beta}^*=0\) (without loss of generality thanks to translation invariance), and let \(\boldsymbol{P}= \boldsymbol{I}_p\). Furthermore, let \(U_i=0\).

• In logistic regression (6), let \((a,\sigma,\gamma)\) be the unique solution to 23 25 , let \(\boldsymbol{P}=\boldsymbol{I}_p - \boldsymbol{w}\boldsymbol{w}^\top\) for \(\boldsymbol{w}=\beta^*/\|\boldsymbol{\beta}^*\|\), and let \(U_i=\boldsymbol{x}_i^T\boldsymbol{w}\). Here, \(\boldsymbol{X}\boldsymbol{P}\) is independent of \((y_i,\ell_i,U_i)_{i\in[n]}\).

Thanks to \(\|\boldsymbol{X}/\sqrt n\|_{op}\to^P 1+\delta^{-1/2}\) and 17 or 27 26 , we have \(\|\boldsymbol{X}\hat{\boldsymbol{b}}(I)\|/\sqrt{I} \le K\) for \(K=2 q^{-1/2}(1+\delta^{-1/2})(a^2+\sigma^2)^{1/2}\) with probability approaching one. Thus \(\mathop{\mathrm{\mathbb{P}}}(\hat{\boldsymbol{b}}(I)=\boldsymbol{\hat{\beta}}(I))\to 1\) for \(\boldsymbol{\hat{\beta}}(I)\) in 45 , so we may argue with \(\boldsymbol{\hat{\beta}}= \boldsymbol{\hat{\beta}}(I)\). Similarly for \(\tilde{I}\) we have \(\mathop{\mathrm{\mathbb{P}}}(\hat{\boldsymbol{b}}(I)=\boldsymbol{\hat{\beta}}(I))\to 1\) for \(\boldsymbol{\hat{\beta}}(\tilde{I})\) in 45 , and we may argue with \(\boldsymbol{\tilde{\beta}}= \boldsymbol{\hat{\beta}}(\tilde{I})\). Let also \(\boldsymbol{\psi},\boldsymbol{\tilde{\psi}}\) be defined in 3 (in particular, we have \(\psi_i=0\) of \(i\notin I\) and \(\psi_i = -\ell_i(\boldsymbol{x}_i^T\hat{\boldsymbol{b}}(I))\) in the high-probability event \(\hat{\boldsymbol{b}}(I)=\boldsymbol{\hat{\beta}}\), and similarly for \(\boldsymbol{\tilde{\psi}},\hat{\boldsymbol{b}}(\tilde{I}),\boldsymbol{\tilde{\beta}}\).)

By 48 and 52 from the auxiliary lemmas, we have \[p \hat{\boldsymbol{\beta}}^\top\boldsymbol{P}\tilde{\boldsymbol{\beta}} = \gamma^2 \boldsymbol{\psi}^\top \tilde{\boldsymbol{\psi}} + O_P(\sqrt n)\] where \(\boldsymbol{\psi}^\top \tilde{\boldsymbol{\psi}}= \sum_{i\in I \cap \tilde{I}}\psi_i\tilde{\psi}_i\). With \(n/p=\delta\) and \(|I\cap \tilde{I}| = nq^2 + O_P(n^{1/2})\) thanks to the explicit formulae for the expectation and variance of the hyper-geometric distribution, we have \[\tilde{\boldsymbol{\beta}}^\top \boldsymbol{P}\hat{\boldsymbol{\beta}} = \delta q^2\gamma^2 \boldsymbol{\psi}^T\boldsymbol{\tilde{\psi}}/|I\cap \tilde{I}| + O_P(n^{-1/2}). \label{eq:key}\tag{34}\] By the Cauchy–Schwarz inequality and the concentration of sampling without replacement (11) (or see 7 for details), the absolute value of \(\boldsymbol{\psi}^T\boldsymbol{\tilde{\psi}}/|I\cap \tilde{I}| = \frac{1}{|I\cap \tilde{I}|}\sum_{i\in I\cap\tilde{I}}\psi_i\tilde{\psi}_i\) is smaller than \[\begin{gather} \Bigl(\frac{1}{|I\cap \tilde{I}|}\sum_{i\in I\cap\tilde{I}}\tilde{\psi}_i^2 \Bigr)^{1/2} \Bigl(\frac{1}{|I\cap \tilde{I}|}\sum_{i\in I\cap\tilde{I}}\psi_i^2 \Bigr)^{1/2} \\\le \Bigl(\frac{1}{|\tilde{I}|}\sum_{i\in \tilde{I}}\tilde{\psi}_i^2 \Bigr)^{1/2} \Bigl(\frac{1}{|\tilde{I}|}\sum_{i\in I}\psi_i^2 \Bigr)^{1/2} +o_P(1) = \frac{\sigma^2}{q\delta\gamma^2} + o_P(1) \end{gather}\] thanks to 17 (in robust regression) or 28 (in logistic regression) for the last equality. Combined with 34 , we have proved \(|\tilde{\boldsymbol{\beta}}^\top\boldsymbol{P}\hat{\boldsymbol{\beta}}| \le \delta q^2\gamma^2 \frac{\sigma^2}{q\delta\gamma^2} + o_P(1)= q\sigma^2 + o_P(1).\) Let \(\bar \mathop{\mathrm{\mathbb{E}}}\) be the conditional expectation given \(({I,\tilde{I}},\boldsymbol{X}\boldsymbol{\beta}^*,\boldsymbol{y})\) (In robust regression, \(\boldsymbol{\beta}^*=0\) so \(\bar\mathop{\mathrm{\mathbb{E}}}\) is the conditional expectation given \(\{{I,\tilde{I}}, (\varepsilon_i)_{i\in[n]}\}\)). Thanks to the the Gaussian Poincaré inequality in 8, \[\bar{\eta}\mathrel{\vcenter{:}}= \sigma^{-2} \bar{\mathop{\mathrm{\mathbb{E}}}}[\tilde{\boldsymbol{\beta}}^\top \boldsymbol{P}\hat{\boldsymbol{\beta}}] \quad\text{ satisfies }\quad \begin{cases} \bar \eta = \tilde{\boldsymbol{\beta}}^\top \boldsymbol{P}\hat{\boldsymbol{\beta}}/\sigma^2 + O_P(n^{-1/2}), \\ |\bar \eta |\le q+ o_P(1). \end{cases} \label{concentration-eta-bar}\tag{35}\] Similarly, by 8 we have \(\bar \mathop{\mathrm{\mathbb{E}}}[ \boldsymbol{\psi}^T \boldsymbol{\tilde{\psi}}] /{|I\cap \tilde{I}|} = \boldsymbol{\psi}^T \boldsymbol{\tilde{\psi}}/{|I\cap \tilde{I}|} + O_P(n^{-1/2})\). Combining the previous displays gives \[\bar \eta = \frac{\delta q^2 \gamma^2}{\sigma^2} \frac{1}{|I\cap \tilde{I}|}\bar{\mathop{\mathrm{\mathbb{E}}}}\Bigl[ \sum_{i\in I\cap\tilde{I}}\psi_i\tilde{\psi}_i \Bigr] + o_P(1). \label{ineq95eta95bar}\tag{36}\] For an overlapping observation \(i\in I\cap \tilde{I}\), using the derivative formula in 3 and the moment inequality in 7 conditionally on \((I,\tilde{I}, \boldsymbol{X}\boldsymbol{\beta}^*,\boldsymbol{y})\) and \((\boldsymbol{x}_l)_{l\ne i}\), applied to the standard normal \(\boldsymbol{P}\boldsymbol{x}_i + \boldsymbol{w}Z\) (for \(Z\sim N(0,1)\) independent of everything else) and \(\boldsymbol{W}=[\boldsymbol{P}\hat{\boldsymbol{\beta}}|\boldsymbol{P}\tilde{\boldsymbol{\beta}}]\in\mathbb{R}^{p\times 2}\), we find for the indicator function \(\mathbb{I}\{i\in I\cap \tilde{I}\}\) that \(\mathbb{I}\{i\in I\cap \tilde{I}\} \mathop{\mathrm{\mathbb{E}}}[\text{LHS}_i \mid I,\tilde{I}] \le C \mathop{\mathrm{\mathbb{E}}}[\sum_{j=1}^p \|\frac{\partial \hat{\boldsymbol{\beta}}}{\partial x_{ij}}\|^2 + \|\frac{\partial \tilde{\boldsymbol{\beta}}}{\partial x_{ij}}\|^2]\) where \[\text{LHS}_i \eqqcolon \Bigl\| \begin{pmatrix} \boldsymbol{x}_i^\top\boldsymbol{P}\hat{\boldsymbol{\beta}}- \mathop{\mathrm{tr}}[\boldsymbol{P}\boldsymbol{A}]\psi_i - (\hat{\boldsymbol{\beta}}^\top \boldsymbol{P}\boldsymbol{A}\boldsymbol{X}^\top \boldsymbol{D})\boldsymbol{e}_i \\ \boldsymbol{x}_i^\top\boldsymbol{P}\tilde{\boldsymbol{\beta}}- \mathop{\mathrm{tr}}[\boldsymbol{P}\tilde{\boldsymbol{A}}]\tilde{\psi}_i - (\tilde{\boldsymbol{\beta}}^\top\boldsymbol{P}\tilde{\boldsymbol{A}}\boldsymbol{X}^\top \tilde{\boldsymbol{D}})\boldsymbol{e}_i \end{pmatrix} - (\boldsymbol{W}^\top\boldsymbol{W})^{1/2} \boldsymbol{g}_i \Bigr\|^2\] for all \(i\in I\cap \tilde{I}\), where \(\boldsymbol{g}_i\sim N(\boldsymbol{0}_2, \boldsymbol{I}_2)\). After summing over \(i\in I\cap\tilde{I}\) and using 50 , we get \(\sum_{i\in I\cap\tilde{I}}\mathop{\mathrm{\mathbb{E}}}[\text{LHS}_i \mid I,\tilde{I}] \le C\) and \[\frac{1}{n} \sum_{i\in I\cap\tilde{I}}\mathop{\mathrm{\mathbb{E}}}\Bigl[\text{LHS}_i \mid I,\tilde{I}\Bigr] \le \frac{C'}{n}, \qquad \frac{1}{n} \sum_{i\in I\cap\tilde{I}}\text{LHS}_i = O_{{p}}\Bigl(\frac{1}{n}\Bigr)\] for some deterministic constants \(C,C'\) independent of \(n,p\).

Using 27 in logistic regression or 17 in robust regression, we know \(\|\boldsymbol{P}\hat{\boldsymbol{\beta}}\|^2\to^P \sigma^2\) and similarly \(\|\boldsymbol{P}\tilde{\boldsymbol{\beta}}\|^2\to^P\sigma^2\), as well as \(\mathop{\mathrm{tr}}[\boldsymbol{P}\boldsymbol{A}]\to^P\gamma\), and \(\mathop{\mathrm{tr}}[\boldsymbol{P}\tilde{\boldsymbol{A}}]\to^P\gamma\) by 52 . Using the Lipschitz inequality for the matrix square root \(\|\sqrt{\boldsymbol{M}}-\sqrt{\boldsymbol{N}}\|_{op} \le \|(\sqrt{\boldsymbol{M}}+\sqrt{\boldsymbol{N}})^{-1}\|_{op} \|\boldsymbol{M}- \boldsymbol{N}\|_{op}\) for positive definite matrices \(\boldsymbol{N},\boldsymbol{M}\) [23] [24] which follows from \(\boldsymbol{x}^T(\sqrt{\boldsymbol{M}}+\sqrt{\boldsymbol{N}})\boldsymbol{x}\lambda = \boldsymbol{x}^T(\boldsymbol{M}-\boldsymbol{N})\boldsymbol{x}\) for any unit eigenvector \(\boldsymbol{x}\) of \(\sqrt{\boldsymbol{M}}-\sqrt{\boldsymbol{N}}\) with eigenvalue \(\lambda\), here we get \[\label{39} \Big\| \begin{pmatrix} 1 & \bar\eta \\ \bar\eta & 1 \\ \end{pmatrix}^{-1} \Big\|_{op} = \frac{1}{1-\bar\eta} \le \frac{2}{1-q} \text{ and } \Big\| \sigma \begin{pmatrix} 1 & \bar\eta \\ \bar\eta & 1 \\ \end{pmatrix}^{1/2} - (\boldsymbol{W}^\top\boldsymbol{W})^{1/2} \Big\|_{op} = o_P(1)\tag{37}\] on the event \(|\bar\eta|\le (1+q)/2<1\) which has probability approaching one thanks to 35 . Using the moment bounds 47 to bound from above \(\sum_{i=1}^n((\hat{\boldsymbol{\beta}}^\top \boldsymbol{P}\boldsymbol{A}\boldsymbol{X}^\top \boldsymbol{D})\boldsymbol{e}_i)^2 = \|\boldsymbol{D}\boldsymbol{X}\boldsymbol{A}\boldsymbol{P}\boldsymbol{\hat{\beta}}\|^2\), we find \[\frac{1}{n} \sum_{i\in I\cap \tilde{I}} \| \begin{pmatrix} \boldsymbol{x}_i^\top\boldsymbol{P}\hat{\boldsymbol{\beta}}- \gamma\psi_i \\ \boldsymbol{x}_i^\top\boldsymbol{P}\tilde{\boldsymbol{\beta}}- \gamma\tilde{\psi}_i \end{pmatrix} - \sigma \begin{pmatrix} 1 & \bar\eta \\ \bar\eta & 1 \\ \end{pmatrix}^{1/2} \boldsymbol{g}_i \|^2 = o_P(1) + o_P(1)\frac{1}{n} \sum_{i=1}^n\bigl(\|\boldsymbol{g}_i\|^2+\psi_i^2+\tilde{\psi}_i^2\bigr),\] and thanks to \(\frac{1}{n} \sum_{i=1}^n(\|\boldsymbol{g}_i\|^2+\psi_i^2+\tilde{\psi}_i^2)= O_P(1)\), the previous display converges to 0 in probability. Since \(\boldsymbol{x}_i^T\boldsymbol{P}\boldsymbol{\hat{\beta}}= \boldsymbol{x}_i^T\boldsymbol{\hat{\beta}}- U_i \boldsymbol{w}^T\boldsymbol{\hat{\beta}}\) for \(U_i = \boldsymbol{x}_i^T \boldsymbol{w} =^d N(0,1)\) and given \(\hat{\boldsymbol{\beta}}^\top\boldsymbol{w}\to^P a\) by 26 , together with \(\frac{1}{n} \sum_{i=1} U_i^2 = O_P(1)\) since \(\sum_{i=1} U_i^2\sim \chi^2_n\) we find \[\frac{1}{n} \sum_{i\in I\cap \tilde{I}} \| \begin{pmatrix} \boldsymbol{x}_i^\top\hat{\boldsymbol{\beta}} - a U_i - \gamma\psi_i - \sigma G_i \\ \boldsymbol{x}_i^\top \tilde{\boldsymbol{\beta}} - aU_i - \gamma\tilde{\psi}_i - \sigma\tilde{G}_i \end{pmatrix} \|^2 = o_P(1) ~~ \text{ where} ~~ \begin{pmatrix} G_i \\ \tilde{G}_i \end{pmatrix} = \begin{pmatrix} 1 & \bar \eta \\ \bar\eta & 1 \\ \end{pmatrix}^{1/2} \boldsymbol{g}_i.\] With probability approaching one, the second term in 45 is 0 for the large enough \(K\) that we took at the beginning, and in this event the modified M-estimator \(\hat{\boldsymbol{\beta}}\) equals to the original M-estimator \(\hat{\boldsymbol{b}}(I)\) so that \(\psi_i = - \ell_i(\boldsymbol{x}_i^\top\hat{\boldsymbol{b}})\) (cf. 2), and similarly for \(\tilde{\boldsymbol{\psi}}\). We have established \[\begin{align} \frac{1}{n}\sum_{i\in I\cap\tilde{I}} \| \underbrace{\boldsymbol{x}_i^\top\hat{\boldsymbol{b}}+ \gamma \ell_i'(\boldsymbol{x}_i^\top\hat{\boldsymbol{b}}) - a U_i - \sigma G_i }_{-\text{Rem}_i} \|^2 = o_P(1). \end{align}\] Define \(\text{Rem}_i\) by \(\boldsymbol{x}_i^\top\hat{\boldsymbol{b}}+ \gamma \ell_i'(\boldsymbol{x}_i^\top\hat{\boldsymbol{b}}) = a U_i + \sigma G_i + \text{Rem}_i\) so that \(\boldsymbol{x}_i^\top\hat{\boldsymbol{b}}= \mathop{\mathrm{prox}}[\gamma\ell_i](a U_i + \sigma G_i + \text{Rem}_i)\) by definition of the proximal operator. Now set \(\hat{p}_i = \mathop{\mathrm{prox}}[\gamma\ell_i](a U_i + \sigma G_i)\). Because \(\mathop{\mathrm{prox}}[\gamma\ell_i](\cdot)\) is 1-Lipschitz, \[\Bigl( \sum_{i\in I\cap \tilde{I}} \|\hat{p}_i - \boldsymbol{x}_i^\top \hat{\boldsymbol{b}}\|^2 \Bigr)^{1/2} \le \Bigl( \sum_{i\in I\cap \tilde{I}} \|\text{Rem}_i\|^2 \Bigr)^{1/2} = o_P({\sqrt n}).\] Similarly, a proximal approximation holds for \(\boldsymbol{x}_i^\top \tilde{\boldsymbol{\beta}}\) using \((U_i,\tilde{G}_i)\) instead. We have to be a little careful here because \(\bar \eta\) is independent of the \((G_i,\tilde{G}_i)\) but not of the \((U_i,y_i)\). Using that \(|\ell_i'|\le 1\), and that \(\bar E[|A-B|] = o_P(1)\) if \(A,B\) are bounded random variables such that \(|A-B|=o_P(1)\), 36 gives \[\begin{align} \bar \eta &= \frac{\delta q^2\gamma^2}{\sigma^2} \sum_{i\in I\cap \tilde{I}} \bar{\mathop{\mathrm{\mathbb{E}}}}\Bigl[ \frac{ \ell_i'(\mathop{\mathrm{prox}}[\gamma\ell_i](a U_i+\sigma G_i)) \ell_i'(\mathop{\mathrm{prox}}[\gamma\ell_i](a U_i+\sigma \tilde{G}_i) ) }{|I\cap \tilde{I}|} \Bigr] + o_P(1) \end{align}\] where inside the conditional expectation \(\bar\mathop{\mathrm{\mathbb{E}}}[\cdot]\), \((\bar\eta,U_i,\ell_i)\) are fixed are integration is performed with respect to the distribution of \((G_i,\tilde{G}_i)\), so that \(\bar\eta = \frac{\delta q^2\gamma^2}{\sigma^2} \varphi(\bar\eta) + o_P(1)\) where \[\begin{align} \varphi(t) &= \frac{1}{|I\cap\tilde{I}|} \sum_{i\in I\cap \tilde{I}} \bar\mathop{\mathrm{\mathbb{E}}}\Bigl[ \ell_i'\Bigl(\mathop{\mathrm{prox}}[\gamma\ell_i](a U_i+\sigma G_i^t)\Bigr) \ell_i'\Bigl(\mathop{\mathrm{prox}}[\gamma\ell_i](a U_i+\sigma \tilde{G}_i^t)\Bigr) \Bigr] \\&= \frac{1}{|I\cap\tilde{I}|} \sum_{i\in I\cap \tilde{I}} \iint \ell_i'\Bigl(\mathop{\mathrm{prox}}[\gamma\ell_i](a U_i+\sigma g)\Bigr) \ell_i'\Bigl(\mathop{\mathrm{prox}}[\gamma\ell_i](a U_i+\sigma \tilde{g})\Bigr) \phi_t(g,\tilde{g})\mathrm{d}g\mathrm{d}\tilde{g} \end{align}\] where \(\phi_t\) is the joint density of two jointly normal \((G^t,\tilde{G}^t)\) with \(\mathop{\mathrm{\mathbb{E}}}[(G^t)^2]=\mathop{\mathrm{\mathbb{E}}}[(\tilde{G}^t)^2]=1\) and \(\mathop{\mathrm{\mathbb{E}}}[G^t\tilde{G}^t]=t\), and in the first line \((G_i^t,\tilde{G}_i^t)\sim \phi_t\) is independent of \((\bar\eta,U_i,y_i,\ell_i)_{i\in[n]}\). Because of the law of large numbers for the deterministic solution \(\eta\) of 12 (with \(a=0\) in robust regression) or ?? (in logistic regression), we have \(\eta = \frac{\delta^2q^2\gamma^2}{\sigma^2}\varphi(\eta) {+o_P(1)}\). Taking the difference between the previous display and this fixed-point equation satisfied by \(\eta\), using the mean-value theorem, \[\bar \eta-\eta = \frac{\delta q^2\gamma^2}{\sigma^2} \Bigl(\varphi(\bar\eta) - \varphi(\eta)\Bigr) + o_P(1) = \frac{\delta q^2\gamma^2}{\sigma^2} \Bigl(\bar\eta-\eta\Bigr) \varphi'(t)+ o_P(1)\] for some \(t\in [\bar\eta,\eta]\). By calculation similar to 15 16 thanks to 1, if \((G_i^t,\tilde{G}_i^t)\) has density \(\phi_t\), \[\begin{align} \varphi'(t)&= \frac{ 1 }{|I\cap \tilde{I}|} \frac{\sigma^2}{\gamma^2} \sum_{i\in I\cap \tilde{I}} \bar \mathop{\mathrm{\mathbb{E}}}\Bigl[ \frac{ \gamma \ell_i''(\mathop{\mathrm{prox}}[\gamma\ell_i](aU_i+\sigma G_i^t)) }{ 1 + \gamma \ell_i''(\mathop{\mathrm{prox}}[\gamma\ell_i](aU_i+\sigma G_i^t)) } \frac{\gamma \ell_i''(\mathop{\mathrm{prox}}[\gamma\ell_i](aU_i+\sigma \tilde{G}_i^t))}{1 + \gamma \ell_i''(\mathop{\mathrm{prox}}[\gamma\ell_i](aU_i+\sigma \tilde{G}_i^t))} \Bigr]\\ &\le \frac{ 1 }{|I\cap \tilde{I}|} \frac{\sigma^2}{\gamma^2} \sum_{i\in I\cap \tilde{I}} \bar \mathop{\mathrm{\mathbb{E}}}\Bigl[ \frac{ \gamma \ell_i''(\mathop{\mathrm{prox}}[\gamma\ell_i](aU_i+\sigma G_i^t)) }{ 1 + \gamma \ell_i''(\mathop{\mathrm{prox}}[\gamma\ell_i](aU_i+\sigma G_i^t)) } \Bigr] \\ &= \frac{ 1 }{|I\cap \tilde{I}|} \frac{\sigma^2}{\gamma^2} \sum_{i\in I\cap \tilde{I}} \int\Bigl[ \frac{ \gamma \ell_i''(\mathop{\mathrm{prox}}[\gamma\ell_i](aU_i+\sigma g)) }{ 1 + \gamma \ell_i''(\mathop{\mathrm{prox}}[\gamma\ell_i](aU_i+\sigma g)) } \Bigr]\frac{e^{-g^2/2}}{\sqrt{2\pi}}\mathrm{d}g \quad\text{since } G_i^t\sim N(0,1) \\ &= \frac{\sigma^2}{\gamma^2 q\delta} + o_P(1), \end{align}\] where we have used the law of large numbers and the nonlinear system with equation 10 in robust regression and equation 25 in logistic regression. Combining the above displays, we are left with \[|\bar\eta - \eta| = |\bar \eta - \eta| \frac{\delta q^2\gamma^2}{\sigma^2} |\varphi'(t)| + o_P(1) \le |\bar \eta - \eta| \frac{\delta q^2\gamma^2}{\sigma^2} \frac{\sigma^2}{\gamma^2 q\delta} + o_P(1) = q|\bar{\eta}-\eta| + o_P(1)\] and \(\bar\eta - \eta = o_P(1)\) thanks to \(q\in (0,1)\). Since \(\bar\eta = \hat{\boldsymbol{\beta}}^\top\boldsymbol{P}\tilde{\boldsymbol{\beta}}/\sigma^2 + o_P(1)\) by 35 , the proof of 18 and 32 is complete. Next, 19 follows from 48 and 52 .

Finally for 20 and 33 , by symmetry \(\mathop{\mathrm{\mathbb{E}}}[\text{LHS}_i \mid I, \tilde{I}]\) is the same for all \(i\in I\cap \tilde{I}\). In particular, the maximum of the conditional expectation is the same as the average over \(I\cap \tilde{I}\), so that \(\sum_{i\in I\cap \tilde{I}}\mathop{\mathrm{\mathbb{E}}}[\text{LHS}_i \mid I, \tilde{I}] \le C\) proved above gives \(\max_{i\in I\cap \tilde{I}} \mathop{\mathrm{\mathbb{E}}}[\text{LHS}_i \mid I, \tilde{I}] = O_P(1/n)\) since \(I\cap \tilde{I}\) has cardinality of order \(n\). Finally, we have \[\boldsymbol{W}^\top\boldsymbol{W}\to^P \sigma^2 \begin{pmatrix} 1 & \eta \\ \eta & 1 \end{pmatrix}, \qquad \Bigl(\boldsymbol{W}^\top\boldsymbol{W}\Bigr)^{1/2}\to^P \sigma \begin{pmatrix} 1 & \eta \\ \eta & 1 \end{pmatrix}^{1/2}, \label{cont95mapping95matrix95sq}\tag{38}\] by continuity of the matrix square root and the continuous mapping theorem (or, alternatively, by reusing the argument in 37 ). Using again \(\mathop{\mathrm{tr}}[\boldsymbol{P}\boldsymbol{A}]\to^P\gamma\), \(\boldsymbol{\hat{\beta}}^T\boldsymbol{w}\to^Pa\), \((\hat{\boldsymbol{\beta}}^\top \boldsymbol{P}\boldsymbol{A}\boldsymbol{X}^\top \boldsymbol{D})\boldsymbol{e}_i\to^P 0\), and similarly for \(\boldsymbol{\tilde{\beta}}\), combined with 38 , we obtain 20 and 33 .

5 Auxiliary lemmas↩︎

5.1 Approximate multivariate normality↩︎

This moment inequality is a matrix-generalization of [25] and [26]. It is particularly useful to show that \[\boldsymbol{W}(\boldsymbol{z})^\top \boldsymbol{z}- \sum_{j=1}^p \frac{\partial \boldsymbol{W}(\boldsymbol{z})^\top \boldsymbol{e}_j}{\partial z_j}\] is approximately multivariate normal with covariance approximated by \(\boldsymbol{W}(\boldsymbol{z})^\top\boldsymbol{W}(\boldsymbol{z})\).

5.2 Derivative of \(F(\eta)\)↩︎

5.3 Modified loss and moment inequalities↩︎

This subsection provides useful approximations to study two estimators \(\hat{\boldsymbol{b}}(I),\hat{\boldsymbol{b}}(\tilde{I})\) trained on two subsampled datasets indexed in \(I\) and \(\tilde{I}\). These approximations are used in the proof of the main result in 4, with the key ingredient being 4. The approximations in this subsection are obtained as a consequence of the moment inequalities given in [lm:chi_square,lm:second_order_stein] developed in [31] for estimating the out-of-sample error of a single estimator. Because the moment inequalities in [lm:chi_square,lm:second_order_stein] requires us to bound from above expectations involving \(\hat{\boldsymbol{b}}(I),\hat{\boldsymbol{b}}(\tilde{I})\) and their derivatives, we resort to the following modification of the M-estimators (introduced in [32]) to guarantee that any finite moment \(\hat{\boldsymbol{b}}(I),\hat{\boldsymbol{b}}(\tilde{I})\) and their derivatives are suitably bounded.

Moment inequalities↩︎

5.4 Concentration of sampling without replacement↩︎

References↩︎