.tocmtchapter

1 Introduction↩︎

In machine learning, statistics, and the natural sciences, cross-validation (CV) [1], [2] is routinely used to compare the performance of learning algorithms [3], [4]. In practice, it is not uncommon to pair CV’s point estimates with uncertainty quantification in the form of estimated standard errors or putative confidence intervals. Yet the validity of such uncertainty quantification has been poorly understood until recently, and it is now understood to be closely related to notions of algorithmic stability [5], [6]. Stability of algorithms has long been studied in the learning theory literature [7], allowing existing stability results to be applied to CV uncertainty quantification for assessing the performance of a single algorithm. However, when comparing two algorithms’ performances, their individual stabilities do not directly translate to the type of stability needed for valid CV uncertainty quantification, raising the question of when such uncertainty quantification is valid. As a concrete empirical illustration, 1 shows CV-based confidence intervals accurately covering the performance of a single stable algorithm but badly failing to cover the comparison between two individually stable instances of that same model with slightly different tuning parameters.

1.0.0.1 Our contributions

This work lies at the interface of algorithmic stability and cross-validation. We demonstrate the importance of considering relative stability by studying the Lasso [8] and its close cousin, soft-thresholding (ST) [9]. In the canonical fixed-dimensional linear regression setting of 3, we tightly characterize the components of relative stability and show that the comparison of two ST or Lasso fits with slightly different tuning parameters does not satisfy relative stability (1), even though the assessment of a single ST or Lasso fit does (2), calling into question the validity of CV confidence intervals for such comparisons. 4 provides more details about the relationship between stability and CV. Simulations in 5 support these conclusions, showing that CV confidence intervals provide accurate coverage of the test error of a single ST or Lasso fit even for moderate sample sizes, while they fail to cover the difference in test errors between two ST or Lasso fits even for very large sample sizes.

1.0.0.2 Related work

The importance of the stability of an algorithm with respect to its generalization error [10] has prompted numerous studies of the stability of popular classes of algorithms [10]–[14]. Across the years, different notions of stability have been introduced [15]–[20] and building upon the domain of algorithmic stability, multiple papers [5], [6], [13], [19], [20] have established interesting relationships between the theoretical properties of cross-validation and the stability properties of the algorithms involved. [5] and [6] derive central limit theorems and consistent variance estimators for the CV estimator under sufficient conditions on algorithmic stability. The former paper does so with the mean-square stability [19] and the latter with loss stability [20], both of which are known to decay to zero for a variety of algorithms. However, to our knowledge, no prior works have assessed the sufficient conditions for asymptotic normality in the case when the asymptotic variance in these central limit theorems goes to zero, as would be expected in the common scenario of comparing the performance of two algorithms that converge to the same prediction rule (e.g., if they are both consistent for the optimal prediction rule). This is the focus of this paper, leading to novel negative results about stability and validity of CV confidence intervals even in very regular settings.

[21] studies the role of algorithmic stability in model comparison, but while our work focuses on drawing inferences about the test error \(R_n\), defined in [eq:rcond-def-rscv], theirs focuses on the expected test error \(\mathbb{E}[R_n]\) and shows that inference concerning \(\mathbb{E}[R_n]\) is often difficult even when inference concerning \(R_n\) is easy. We note that some recent works [22]–[24] have studied various other aspects of asymptotic distributional properties of CV, but none present negative results comparable to ours.

1.0.0.3 Notation

For each \(n\in\mathbb{N}\), we define the set \([n] \triangleq\{1,\dots, n\}\). We denote by \(\lambda_{\mathrm{min}}\mathopen{}\mathclose{\left({\mathbf{A}}}\right)\) the minimum eigenvalue of a matrix \(\mathbf{A}\). For deterministic sequences \((f_n)_n\) and \((g_n)_n\), \(f_n = o(g_n)\) if \(\frac{f_n}{g_n} \rightarrow 0\) as \(n \rightarrow \infty\), and \(f_n = O(g_n)\) if \(\frac{f_n}{g_n}\) is asymptotically bounded. Following canonical notations, we write \(f_n = \omega(g_n)\) to mean that \(g_n = o(f_n)\) as \(n \rightarrow \infty\), we write \(f_n = \Omega(g_n)\) to mean that \(g_n = O(f_n)\) as \(n \rightarrow \infty\), and we write \(f_n = \Theta(g_n)\) to mean that \(f_n = O(g_n)\) and \(f_n = \Omega(g_n)\) as \(n \rightarrow \infty\). Finally, we write \(f_n \sim g_n\) to mean that \(\frac{f_n}{g_n} \rightarrow 1\) as \(n \rightarrow \infty\).

2 Preliminaries↩︎

Before presenting our results, we establish some necessary definitions, largely following the notation and nomenclature of [6]. We will consider a sequence \((Z_i)_{i\geq 0}\) of random data points taking values in a set \(\mathcal{Z}\) and a scalar loss function \(h_n(Z_0,\mathbf{Z})\) where \(\mathbf{Z}\) is a training set of size \(n\). A typical choice for \(h_n\) in the regression setting is squared error loss, \[h_n(Z_0, \mathbf{Z}) = (Y_0 - \hat{f}(X_0; \mathbf{Z}))^2,\] applied to the predicted response value of a test point \(Z_0 = (X_0, Y_0)\), obtained from an algorithm fitting a prediction rule \(\hat{f}(\cdot; \mathbf{Z})\) to training data \(\mathbf{Z}\). When comparing the performance of two algorithms, we will choose \(h_n\) to be the difference between the losses of two prediction rules. In order to ensure a smooth read when we switch between the settings of single algorithm assessment and comparison of algorithms, we will make the distinction clear by adding a superscript to \(h_n\): \(h_n^{\mathrm{sing}}\) and \(h_n^{\mathrm{diff}}\), respectively. In addition, our asymptotic statements should all be interpreted as taking \(n\to\infty\).

For the purpose of illustrating the importance of considering stability in a relative sense rather than an absolute sense, we will now define a notion of relative stability based on loss stability. We define loss stability for algorithms whose learned prediction rules do not depend on the order of the training points, which will be our focus here.

We introduced these quantities for a function \(h_n\), but we will generically refer to the loss stability and the relative loss stability of an algorithm or a comparison of algorithms when \(h_n\) is clear from context. Note that we include the factor of \(n\) in the numerator of [eq:rel95stab] because it sets the baseline rate of \(r(h_n)\) to constant order: we will say that an algorithm or a comparison of algorithms satisfies the relative loss stability condition if \(r(h_n) = o(1)\), which is equivalent to a key sufficient condition for the central limit theorem and consistent variance estimation for CV proved in [6]; see 4 for more on the connection between relative stability and CV. We will illustrate the importance of relative stability for CV by studying the soft-thresholded least squares regression and Lasso algorithms in the linear regression setting.

Throughout, we will consider \(\textrm{i.i.d.}\) data points \(Z_i = (X_i, Y_i)\) from the linear model align Y_i = X_i^^+ _i,
X_i~(0,), _i ~(0, ^2), _i X_i, parametrized by the unknown vector \(\beta^\star\in \mathbb{R}^p\) and noise level \(\tau > 0\). Here, \(\mathbf{Y} = (Y_1, \ldots, Y_n) \in \mathbb{R}^n\) is the vector of response variables or targets, \(\mathbf{X} = (X_1, \ldots, X_n)^\top \in \mathbb{R}^{n \times p}\) is the matrix of regressors or features, and \(\varepsilon = (\varepsilon_1, \ldots, \varepsilon_n) \in \mathbb{R}^n\) is the noise vector.

The loss function considered for a linear prediction rule will be the squared error loss \[h_n^{\mathrm{sing}}(Z_0, \mathbf{Z}) \triangleq(Y_0 - X_0^\top \hat{\beta})^2,\] where the estimated parameter vector \(\hat{\beta}\) is learned from the training set \(\mathbf{Z} = (Z_1, \ldots, Z_n)\). When the focus is on the comparison of two prediction rules, the loss function will be defined as the difference of two such individual losses in the form \[\label{eq:hdiff} h_n^{\mathrm{diff}}(Z_0, \mathbf{Z}) \triangleq(Y_0 - X_0^\top \hat{\beta}^{(1)})^2 - (Y_0 - X_0^\top \hat{\beta}^{(2)})^2\tag{1}\] for \(\hat{\beta}^{(1)}\) and \(\hat{\beta}^{(2)}\) both learned on the training set \(\mathbf{Z}\).

A classical way to estimate \(\beta^\star\) is the ordinary least squares (OLS) estimator defined as align _ols(^)^-1 ^. To simplify notation, we will leave the dependence on the sample size \(n\) implicit even though our asymptotic results will involve \(\hat{\beta}_{\mathrm{\small ols}}\). When we expect the parameter vector \(\beta^\star\) to exhibit some level of sparsity, that is to say it has some number of zero entries, a popular estimator used is the Lasso estimator [8], for some choice of penalization parameter \(\lambda\) which determines the level of sparsity in the learned parameter vector. As a convenient stepping stone for our analysis, we first study soft-thresholding, a close cousin of the Lasso.

The close relationship between ST and the Lasso when \(\mathbf{X}^\top \mathbf{X}/n\) is close to the identity, formalized in the following lemma, will allow us to derive (in)stability results for the Lasso from our ST (in)stability results.

The proof of 1 in 9 follows from viewing \(\hat{\beta}^{\mathrm{\small lasso}}_{\lambda_n}\) and \(\hat{\beta}_{\lambda_n}\) as the optimizers of closely related objective functions and using the optimizer comparison lemma of [25] to deduce their proximity.

3 Main Results↩︎

We now state the main results of this work. Our primary theoretical result is that even a simple learning algorithm (ST) in a simple, well-behaved learning setting can fail to generate relatively stable comparisons.

The proof of 1 can be found in 10. Notably, since this main result is a lower bound, the stringent assumptions like linearity, sparsity, covariate independence, and Gaussianity serve to strengthen the result, indicating that failure occurs even in this best-case scenario.

Perhaps surprisingly, under the same conditions, we have relative stability for single algorithm assessment.

The proof of 2 can be found in 13. This secondary result is also significant as it shows that a CV user can be easily duped into thinking that confidence intervals that rely on stability will yield valid inference in the algorithm comparison setting simply because an algorithm is stable when considered in isolation.

We can think of 1 as a stylized version of a setting where one wants to compare two similar machine learning algorithms, such as when the two only differ by a tuning parameter. Note that \(\lambda_n = O(\sqrt{n})\) implies \(\lambda_n = o(n)\), and that \(\Theta(1)\) is also \(o(n)\) as it is asymptotically lower- and upper-bounded by a constant, which means that \(\lambda_n + \delta_n = o(n)\) under the conditions of 1 and thus both ST(\(\lambda_n\)) and ST\((\lambda_n + \delta_n)\) individually satisfy the relative loss stability condition thanks to 2. So, taken together, 1 2 show that even if two learning algorithms are individually well-behaved, their comparison may not be, even when the data comes from a very regular distribution.

We now discuss the penalty parameter regimes appearing in 1 2. For simulations with features and targets sampled in the same conditions as the theorems, we observed that the values selected for \(\lambda_n\) via CV are concentrated around a constant times \(\sqrt{n}\). It therefore makes sense to compare two versions of ST with penalization of order \(\sqrt{n}\) in 1, and we do so by setting the base level of penalization to \(\lambda_n\) of order \(\sqrt{n}\) and parameterizing the difference in penalization of the ST algorithms by \(\delta_n\) of order \(1\). Note that both \(\lambda_n\) and \(\delta_n\) are assumed deterministic in the theorems, but we will present simulations with stochastic \(\lambda_n\) selected via inner CV in 5. Under some regularity conditions on the features, [26] proved that choosing \(\lambda_n = o(n)\) ensures weak consistency of the Lasso estimator for \(\beta^\star\), \(\textrm{i.e.}\) it converges in probability to \(\beta^\star\), and it is therefore natural that the regimes we study are always within this weak consistency regime. As for the \(\sqrt{n}\) order of the penalization specific to our primary result, it has been shown to be a regime of interest for variable selection consistency [27], [28].

The powerful Lasso-ST proximity bound of 1 allows us to translate 1 2 into identical results for the popular Lasso algorithm, showing that our conclusions about ST are by no means specific to that method.

Our proof in 16 combines the ST instability bounds of 1 with the powerful Lasso-ST proximity bound of 1.

As with ST, Lasso comparison instability occurs even though the Lasso algorithm itself is relatively stable.

Our Lasso stability proof in 17 may be of independent interest as the Lasso is known to be unstable under the more stringent notion of uniform stability [29].

4 Importance of Relative Stability for Cross-validation↩︎

To connect our results on relative stability back to CV and prepare for our numerical experiments, we need to introduce some further notation. We have been using \(n\) for the size of the training sets used in the iterations of CV, while [6] use it for the sample size of the larger set of all the data on which CV is run. For the sake of simplicity, we will write \(k\) instead of \(k_n\) to denote the number of folds even though it can depend on \(n\) (leave-one-out CV corresponds to \(k=n+1\)), and we will assume that \(k-1\) evenly divides \(n\). The full data sample size is then simply equal to \(\frac{nk}{k-1}\).

Consider \(k\) vectors of integers, \(\{B_j'\}_{j=1}^k\), each of length \(\frac{n}{k-1}\), whose elements partition \([\frac{n k}{k-1}]\). For each \(B'_j\), define \(B_j\) as a vector of the \(n\) indices in \([\frac{n k}{k-1}]\) that are not in \(B'_j\), so that we can consider each \((B_j,B'_j)\) as a train-validation split. For \(B\) a vector of indices in \([\frac{n k}{k-1}]\), we denote by \(Z_B\) the subvector of \((Z_1, \ldots, Z_{\frac{n k}{k-1}})\) corresponding to the entries of \(B\). Then for a scalar loss function \(h_n(Z_i,Z_B)\), we define the \(k\)-fold cross-validation error align _n _j=1^k _iB_j’ h_n(Z_i, Z_B_j) and the inferential target, the \(k\)-fold test error align R_n _j=1^k _iB_j’ [h_n(Z_i, Z_B_j)Z_B_j], where the conditioning in \(\mathbb{E}[h_n(Z_i, Z_{B_j})\mid Z_{B_j}]\) is on the data points from the \(j\)-th validation set \(Z_{B_j}\), and thus the expectation is taken over only the test point \(Z_i\), as the function \(h_n\) is treated as non-random.

In our notation, [6] use the stability condition \(\gamma(h_n) = o(\frac{\sigma^2(h_n)}{n})\), equivalent to \(r(h_n) = o(1)\), to prove the central limit theorem align (_n- R_n) (0,1). Along with an estimator \(\hat{\sigma}_n^2(h_n)\) provided in [6] and proved to be consistent for \(\sigma^2(h_n)\) therein, this central limit theorem enables the construction of asymptotically valid confidence intervals for \(R_n\). Note that it is indeed possible to use the training sample size in the denominator of the stability condition rather than the full sample size of the dataset on which CV is run, as \(\frac{n k}{k-1} = \Theta(n)\) for any choice of \(k>1\).

When assessing a single algorithm, unless we are in a fully noiseless setting, we might expect \(\sigma^2(h_n^{\mathrm{sing}})\) to be of constant order in general. This means the loss stability condition simplifies to a condition on stability in the absolute sense: \(\gamma(h_n^{\mathrm{sing}}) = o(1/n)\). For instance, we show in 4 that in the linear model with noise, for any linear predictor satisfying some consistency condition, \(\sigma^2(h_n^{\mathrm{sing}})\) converges to a positive constant. However, when comparing two consistent algorithms, we then expect \(\sigma^2(h_n^{\mathrm{diff}})\) to go to \(0\) for algorithms whose performances become increasingly similar when the sample size grows, and this is when reasoning about stability in an absolute sense, as has been the focus in past literature, becomes insufficient. In fact, in 1 it turns out that \(\gamma(h_n^{\mathrm{diff}}) = O(1/n^2)\) (see 12), so the ST comparison is loss stable in the absolute sense. However, the relative loss stability condition does not hold because it properly accounts for the fact that \(\sigma^2(h_n^{\mathrm{diff}})\) goes to zero at a \(1/n^2\) rate.

5 Numerical Experiments↩︎

We performed numerical experiments to empirically confirm the theoretical results of 3. We sampled the features, the independent noise terms and the target variables from the linear model [eq:linear-model] with parameter vector \(\beta^\star= (3, 1, -5, 3, 0, 0, 0, 0, 0, 0)\) of dimension \(10\), and with noise level \(\tau = 10\). We fix \(k = 10\). To satisfy the assumptions of 2 1, we choose \(\lambda_n = \sqrt{n}\) for the base level of penalization, and when comparing algorithms, we set \(\delta_n = 1\) for the difference in the penalization parameters. To explore the asymptotic regime in our simulations, we work with \(n\) ranging from \(90\) to \(90{,}000\) (so the total sample size \(\frac{nk}{k-1}\) ranged from \(10^2\) to \(10^5\)). We used Monte Carlo estimation to compute both \(\sigma^2(h_n)\) and \(\gamma(h_n)\), leveraging 7 8 proved in 18. We provide Python code replicating all experiments at https://github.com/alexandre-bayle/ricv and additional details about the experiments in 18.

We present two types of plot. The first type displays the rates for \(\sigma^2(h_n)\), \(\gamma(h_n)\) and \(r(h_n)\) on the log–log scale by plotting their empirical values for increasing \(n\) with dots and plotting lines for the corresponding rates predicted by our theory. We display the values with a \(\pm \, 2\) standard error confidence band, with details on how to obtain it for \(r(h_n)\) in 18. Thanks to the large number of Monte Carlo replications used, the error bars are very small and thus are not visible. Note that we use the \(x\)-axis labels \(n / 900\) to provide a better scale for visualization. \(n / 900\) goes up to \(10^2\), which is consistent with \(n\) going up to \(90{,}000\). For the second type of plot, using kernel density estimation (KDE), we plot the probability density function across sample sizes of both \(\frac{\sqrt{\frac{n k}{k-1}}}{\sigma(h_n)} (\hat{R}_n- R_n)\) and \(\frac{\sqrt{\frac{n k}{k-1}}}{\hat{\sigma}_n(h_n)} (\hat{R}_n- R_n)\), where \(\hat{\sigma}_n^2(h_n)\) is the within-fold variance estimator introduced in [5] and proved to be consistent for \(\sigma^2(h_n)\) under the relative loss stability condition in [6]. We expect convergence in distribution to \(\mathcal{N}(0, 1)\) under the relative loss stability condition thanks to the combination of results of [6], we thus shade the area below the curve of the probability density function of \(\mathcal{N}(0, 1)\) to make it clearer when the probability density function curves match or not. From its definition [eq:rcond-def-rscv], note that \(R_n\) is straightforward to compute in the simulations thanks to 8.

The simulation results for ST are presented in 2. For the single algorithm assessment of ST, the rates of \(\sigma^2(h_n^{\mathrm{sing}})\), \(\gamma(h_n^{\mathrm{sing}})\), and \(r(h_n^{\mathrm{sing}})\) are constant order, \(1/n^2\) order, and \(1/n\) order, respectively, as stated in 2, and for the algorithm comparison of ST, when \(\delta_n = 1\), we have the expected \(1/n^2\) rate for \(\sigma^2(h_n^{\mathrm{diff}})\) and we actually observe that \(\gamma(h_n^{\mathrm{diff}})\) and \(r(h_n^{\mathrm{diff}})\) seem to be scaling as \(1/(n^2 \sqrt{n})\) and \(\sqrt{n}\), respectively, even though 1 only established them being \(\Omega\) of these rates. As we can see for both choices of the dividing standard deviation in the KDE plots of 2, the asymptotic distribution seems to be Gaussian, but the asymptotic variance does not go to \(1\) when the relative loss stability condition does not hold, that is to say in the comparison setting. The variance estimator \(\hat{\sigma}^2_n(h_n^{\text{diff}})\) has been proved to be a consistent estimator of the targeted variance of \(\sqrt{\frac{n k}{k-1}} (\hat{R}_n- R_n)\) under the loss stability condition in [6]. In the setting we explored, the condition does not hold for the comparison, and we observed empirically that \(\hat{\sigma}^2_n(h_n^{\text{diff}})\) underestimates the targeted variance of \(\sqrt{\frac{n k}{k-1}} (\hat{R}_n - R_n)\), and overestimates \(\sigma^2(h_n^{\text{diff}})\). While the intervals proposed in [6] are valid when the loss stability condition holds, they will not be wide enough when \(\hat{\sigma}^2_n(h_n^{\text{diff}})\) underestimates the targeted variance of \(\sqrt{\frac{n k}{k-1}} (\hat{R}_n - R_n)\), leading to undercoverage and hence asymptotic invalidity.

Next, in light of our analogous theoretical results for the Lasso, we provide analogous simulations for the Lasso as well, though we make them even more realistic by choosing \(\lambda_n\) via inner cross-validation. In particular, we ran simulations for the Lasso with \(\lambda_n\) selected via an inner CV (see 18) for each of the \(k\) iterations of the CV run, still with constant order \(\delta_n = 1\) for the comparison. As mentioned in 3, we actually observed in simulations that the values selected for \(\lambda_n\) are concentrated around a constant times \(\sqrt{n}\). The results for this setting are displayed in 3 and confirm that the same conclusions hold empirically for the cross-validated Lasso as for ST.

We note that the dichotomy exhibited by ST and Lasso is not universal: there are instances when an algorithm satisfies the relative loss stability condition both in its individual form and in the comparison setting. One example of this is ridge regression and we present the corresponding simulations in 4 in 19. [10] proved that ridge regression, with bounded targets, has \(O(1/n)\) uniform stability. This means it has \(O(1/n^2)\) loss stability by [19] and [20]. In the simulations, we see that for individual ridge, with no boundedness assumption, with isotropic features, loss stability scales as \(1/n^2\) and the relative loss stability condition then holds since \(\sigma^2(h_n^{\mathrm{sing}})\) is of constant order. And loss stability scales as \(1/n^4\) in the comparison setting, which, when compared to the observed \(1/n^2\) rate of \(\sigma^2(h_n^{\mathrm{diff}})\), means the relative loss stability condition also holds for comparison.

As a matter of fact, when \(\beta^\star\) has no zero coefficients, the ST estimator can also be an example of an algorithm which satisfies the relative loss stability condition in both its individual form and in the comparison setting. The theory sheds light on the importance of the zero coefficients in the true parameter vector. When \(\beta^\star\) has no zero coefficients, \(\textrm{i.e.}\) \(\|\beta^\star\|_0 = p\), ST actually becomes stable for the algorithm comparison setting. The results of the simulations for this setting, with the choice \(\beta^\star= (3, 1, -5, 3, 4, -3, 10, 8, 5, 2)\), are presented in 5 in 19 and show how the convergence rate of \(\gamma(h_n^{\mathrm{diff}})\) changes compared to the \(\|\beta^\star\|_0 < p\) setting. It now scales as \(1/n^4\), which means that ST satisfies the relative loss stability condition \(r(h_n^{\mathrm{diff}}) = o(1)\) in the comparison setting, since \(\frac{n^2}{\delta_n^2} \, \sigma^2(h_n^{\mathrm{diff}})\) still goes to \(4 \tau^2 \|\beta^\star\|_0\) when \(\|\beta^\star\|_0 = p\). Nonetheless, we reiterate that even a single zero coefficient in \(\beta^\star\) leads to instability for ST, and more generally Lasso, in the comparison setting.

6 Conclusions and Future Work↩︎

Cross-validation is a powerful tool, but given its widespread use for comparing and selecting models, scrutiny of its statistical properties is critical for safe model deployment. This work highlights the importance of relative stability for CV and the challenges posed by relative instability for model comparison. In particular, we proved that even simple, absolutely-stable learning algorithms can generate relatively unstable comparisons. In practice, this led to invalid and highly misleading confidence intervals for the test error difference with \(\sigma^2(h_n^{\mathrm{diff}})\) being well below the targeted variance of \(\sqrt{\frac{n k}{k-1}} (\hat{R}_n- R_n)\). Since CV is often used to conduct formal hypothesis tests for an improvement in test error between two learning algorithms [6], [30]–[34], our work shows that such tests can be misleading even for simple, absolutely stable algorithms and that method developers and consumers should first verify the relative stability of a comparison before applying them.

This paper uses ST and the Lasso to illustrate the dichotomy between algorithm evaluation and comparison when using CV for uncertainty quantification. While it is true that we expect this dichotomy to extend to other ML algorithms as well, we do not attempt to make any claims concerning other ML algorithms in this work. Importantly, we did not aim to show that the CV central limit theorem [eq:clt] is always a poor choice for algorithm comparison. Indeed, 5 presented examples (ST with fully dense \(\beta^*\) and ridge regression) in which relative comparisons are stable. That said, what we have shown is that even a simple ML algorithm, in the linear model setting, applied to very well-behaved data, can fail to satisfy relative stability in the comparison setting, which we hope is enough to at least convince users of CV that they should not expect by default that relative stability holds when comparing two algorithms (even if they are individually stable), which we feel is an important and practical realization that was previously unknown. The main goal is to increase awareness of the pitfalls of CV, highlighting how simple it is for it to be misleading, especially if not studied through the proper lens of relative stability.

However, this work is not without its limitations. First, our analyses are fairly specific to our particular data distribution and ST/Lasso models. Establishing broad, easily verified conditions under which an algorithm comparison is or is not relatively stable is an important direction for future work. Second, while we prove the relative instability of ST/Lasso comparisons and demonstrate the invalidity of their CV confidence intervals, we leave open the question of whether relative instability always implies CV invalidity.

The focus of this work is on exposing a surprising failure mode of the commonly used CV procedure, and not on identifying the best inference procedure for test error. This is why our experiments, designed principally to corroborate our theory, focused on CV intervals alone. While we have shown that the CV central limit theorem [eq:clt] and hence the CV confidence interval construction of [6] can break down in the presence of relative unstable comparisons, we do however introduce below 5 which provides an alternative (possibly conservative) CI construction that yields asymptotic validity even if the comparison is not relatively stable. Specifically, it yields validity whenever each algorithm is individually relatively stable or, more generally, whenever one can construct a valid interval separately for each algorithm’s test loss.

Proof align &_n(R_n^(1)-R_n^(2) )
&-_n(R_n^(1)
& R_n^(2) )
&-_n(R_n^(1) )
& -_n(R_n^(2) )
&-/2-/2 = 1-. 0◻
This approach would ensure valid asymptotic coverage under individual algorithm stability without requiring any additional stability assumption on the comparison. However, the interval could also be significantly wider than the interval derived from [6], due to strong positive correlations between \(\hat{R}_n^{(1)}\) and \(\hat{R}_n^{(2)}\) ignored in the construction of 5. An open question for the reader is whether one can derive tighter confidence intervals for algorithm comparisons when it is only known that each algorithm is individually stable.

Impact Statement↩︎

By highlighting a surprising failure mode of a commonly used procedure for quantifying confidence in the difference between learning algorithms, this paper’s potential broader impact is to reduce the overinterpretation of small empirical CV differences between two learning algorithms, helping to more rigorously distinguish legitimate improvements from inconsequential changes.

.tocmtappendix

7 Additional Notation↩︎

Let \(\stackrel{\mathrm{a.s.}}{\longrightarrow}\) denote almost sure convergence. Let \(\mathbb{1}\mathopen{}\mathclose{\left[{A}}\right]\) denote the indicator function of a subset \(A\). We will denote by \(\Phi\) the cumulative distribution function of the standard Normal and by \(\varphi\) its probability density function. We define the sign function as \(\text{sign}(x) = \frac{x}{|x|} \mathbb{1}\mathopen{}\mathclose{\left[{x \neq 0}}\right]\) and the positive part as \(x_+ = \max(x, 0)\). We write \(\mathbf{M} \sim W_p^{-1}(\Sigma, n)\) to indicate \(\mathbf{M}\) follows the inverse-Wishart distribution with \(n\) degrees of freedom and scale matrix \(\Sigma \in \mathbb{R}^{p \times p}\).

8 Experimental Details for Figure 1↩︎

The experimental setup for 1 is very similar to the one described in the first paragraph of 5. We consider the Lasso estimator here, with \(\lambda_n = \sqrt{n}\) for the base level of penalization, and when comparing algorithms, \(\delta_n = 1\) for the difference in the penalization parameters, where \(\beta^\star= (3, 1, -5, 3, 0, 0, 0, 0, 0, 0)\). For the largest sample size under consideration, we are plotting the actual coverage probability of the confidence interval \(\hat{R}_n\pm q_{1-\alpha/2} \, \sigma(h_n) / \sqrt{\frac{n k}{k-1}}\), over the full range \(\alpha \in [0, 1]\), where \(q_{1-\alpha/2}\) is the \((1 - \alpha/2)\)-quantile of the standard normal distribution, built from the CV central limit theorem of [6] using the true variance \(\sigma^2(h_n)\), against the target coverage, in the single algorithm setting and the comparison setting.

9 Proof of 1:↩︎

Define the objective functions align f_1() (̍^)^-1^- ̍_2^2 + ̍̍_1 f_2() -̍ ̍_2^2 + ̍̍_1 so that \(\hat{\beta}_{\lambda_n}\in \mathop{\mathrm{arg min}}_{\beta\in\mathbb{R}^p}f_1(\beta)\) and \(\hat{\beta}^{\mathrm{\small lasso}}_{\lambda_n}\in \mathop{\mathrm{arg min}}_{\beta\in\mathbb{R}^p}f_2(\beta)\). For any \(\beta\in\mathbb{R}^p\) we have align f_2() - f_1() &= ^- ^- (- (^)^-1^)
&= (- (^/n)^-1) (^- ^)/n
&= (^/n - ) (- (^)^-1^) = (^/n - ) (- _ols). Moreover, by the definitions of the operator norm and \(\hat{\beta}_{\lambda_n}\) (2), align f̍_2(__n) - f_1(__n)̍_2 &= (̍^/n - ) (__n- _ols)̍_2
& ^̍/n - ̍_ _̍_n-_ols̍_2 ^̍/n - I̍_ _n/n. Finally, since \(f_2\) is \(\mu_n\) strongly convex, the optimizer comparison lemma of [25] implies that \(\mu_n \|{\hat{\beta}_{\lambda_n}- \hat{\beta}^{\mathrm{\small lasso}}_{\lambda_n}}\|_2^2 \leq \|{\hat{\beta}_{\lambda_n}- \hat{\beta}^{\mathrm{\small lasso}}_{\lambda_n}}\|_2 \|{\nabla f_2(\hat{\beta}_{\lambda_n}) - \nabla f_1(\hat{\beta}_{\lambda_n})}\|_2\), yielding the first result.

Now fix any \(q\in \mathbb{N}\), and suppose \(X_i\overset{\textrm{\tiny\textrm{i.i.d.}}}{\sim}\mathcal{N}(0,\mathbf{I})\) and \((n-p+1)/2 > 2q\). Then \(V = (\mathbf{X}^\top \mathbf{X})^{-1}\) has an inverse-Wishart distribution with \(n\) degrees of freedom, and each diagonal entry \(V_{jj}\) has an inverse-gamma distribution with shape \(=\frac{n-p+1}{2}\) and scale \(=\frac{1}{2}\) [35]. Therefore, by Jensen’s inequality and the moment formula for an inverse-gamma, align [1/_n^2q] &= [_(nV)^2q] [\(\operatorname{tr}\)(nV)^2q] [p^2q-1_j=1^p (nV_jj)^2q] = (np)^2q[V_11^2q]
&= ()^2q ()^2q p^2q = O(1).

Next let \(W = \mathbf{X}^\top\mathbf{X}\) so that each entry \(W_{jk} = \sum_{i=1}^n X_{ij}X_{ik}\). Then, we may apply Jensen’s inequality, the Marcinkiewicz-Zygmund [36] inequality, Jensen’s inequality again, and finally the moment formula for a Gaussian random variable to find that align [^̍/n - ̍_^2q] &= [^̍/n - ̍_F^2q] [_j=1^p_k=1^p |W_jk-[W_jk]|^2q]
& (p[|X_11^2-1|^2q] + (p^2-p) [|X_11X_12|^2q])
& (p2^2q-1(1+[|X_11|^4q]) + (p^2-p) [|X_11X_12|^2q])
&= (p2^2q-1 (1+(4q-1)!!) + (p^2-p)((2q-1)!!)^2) = O(1/n^q). The second advertised result now follows by Cauchy–Schwarz.

10 Proof of 1:↩︎

1 follows immediately from the following two propositions, proved in [sec:rate-sigma-comp-stlasso] [sec:rate-lstab-comp-stlasso], respectively. Note that the first proposition holds for \(\lambda_n = o(n)\) and \(\delta_n = o(n)\), and does not require the assumption \(\|\beta^\star\|_0 < p\), which makes this proposition a stronger result than what is needed for the proof of 1 assuming \(\lambda_n = O(\sqrt{n})\), \(\lambda_n = \omega(1)\), \(\delta_n = \Theta(1)\) and \(\|\beta^\star\|_0 < p\).

11 Proof of 6:↩︎

We start by introducing a lemma which provides key equations in the comparison setting.

Proof The first equation follows from the first equation of 3. The remaining equations are then derived from there using the same arguments as those mentioned in 3. 0◻

We will show that \[\frac{n}{\delta_n} \mathbb{E}[\hat{\beta}_{\lambda_n + \delta_n} - \hat{\beta}_{\lambda_n}] \rightarrow - \text{sign}(\beta^\star)\] and \[\frac{n}{\delta_n} \mathbb{E}[\hat{\beta}_{\lambda_n + \delta_n} \hat{\beta}_{\lambda_n + \delta_n}^\top - \hat{\beta}_{\lambda_n} \hat{\beta}_{\lambda_n}^\top] \rightarrow - (\text{sign}(\beta^\star) \beta^{\star\top}+ \beta^\star\text{sign}(\beta^\star)^\top)\] where \(\text{sign}(\beta^\star) = (\text{sign}(\beta^\star_i))_{i \in [p]}\), in order to conclude that \(\frac{n^2}{\delta_n^2} \, \sigma^2(h_n^{\mathrm{diff}}) \rightarrow 4\tau^2 \|\beta^\star\|_0\).

Indeed, if the convergences of these two expectations hold, starting from the expression of \(\sigma^2(h_n)\) in 2, since \(\frac{n}{\delta_n} \mathbb{E}[\hat{\beta}_{\lambda_n + \delta_n} - \hat{\beta}_{\lambda_n}]\) and \(\frac{n}{\delta_n} \mathbb{E}[\hat{\beta}_{\lambda_n + \delta_n} \hat{\beta}_{\lambda_n + \delta_n}^\top - \hat{\beta}_{\lambda_n} \hat{\beta}_{\lambda_n}^\top]\) are non-random, we can expand the square, use linearity of expectation, take the limits and factorize back to obtain the following convergence \[\frac{n^2}{\delta_n^2} \sigma^2(h_n^{\mathrm{diff}}) \rightarrow \mathbb{E}[(2 (Y_0 X_0^\top - \beta^{\star\top}\mathbb{E}[X_0 X_0^\top]) (-\text{sign}(\beta^\star)) + \operatorname{tr}((X_0 X_0^\top - \mathbb{E}[X_0 X_0^\top]) (\text{sign}(\beta^\star) \beta^{\star\top}+ \beta^\star\text{sign}(\beta^\star)^\top)))^2]\] where, for \(Y_0 = X_0^\top \beta^\star+ \varepsilon_0\) with \(\mathbb{E}[X_0] = 0\) and \(\mathrm{Var}(X_0) = \mathbf{I}\), align &[(2 (Y_0 X_0^- ^[X_0 X_0^]) (-(^)) + \(\operatorname{tr}\)((X_0 X_0^- [X_0 X_0^]) ((^) ^+ ^(^)^)))^2]
&= [(- 2 Y_0 X_0^(^) + 2 ^(^) + 2 X_0^^X_0^(^) - 2^(^)))^2]
&= [(- 2 _0 X_0^(^))^2]
&= 4[_0^2] [(X_0^(^))^2] _0, X_0
&= 4 ^2 ̍^̍_0 since align [(X_0^(^))^2] = ((^)^X_0) = (^)^(X_0) (^) = (^)^(^) = ̍^̍_0.

We have for \(i=1, \ldots, p,\) align __n + _n, i - __n, i &= (_ols,i) (|_ols,i| - )_+ - (_ols,i) (|_ols,i| - )_+
&= -(_ols,i)

.

Since \(\hat{\beta}_{\mathrm{\small ols}}\mid \mathbf{X} \sim \mathcal{N}(\beta^\star, \tau^2 (\mathbf{X}^\top \mathbf{X})^{-1})\), we can write \(\hat{\beta}_{\mathrm{\small ols},i}= \beta^\star_i + \tilde{\tau}_n Z\) where \(\tilde{\tau}_n = \frac{\tau}{\sqrt{n}} \sqrt{(\frac{\mathbf{X}^\top \mathbf{X}}{n})_{i,i}^{-1}}\) and \(Z \mid \mathbf{X} \sim \mathcal{N}(0, 1)\). Note that we could have \(i\) as a subscript of \(\tilde{\tau}_n\) and \(Z\), but we will only consider one \(i\) at a time in our computations and we can thus omit this subscript for both of them for the sake of notational simplicity, and we will also omit it for some additional notation we define in the rest of the proof.

We now show that \(\frac{n}{\delta_n} \mathbb{E}[\hat{\beta}_{\lambda_n + \delta_n, i} - \hat{\beta}_{\lambda_n, i}] \rightarrow - \text{sign}(\beta^\star_i)\).

Using the law of total expectation, align &[__n + _n, i - __n, i ]
&= - (_ols,i> ) + (_ols,i< - )
& - [_ols,i- _ols,i, ]   (_ols,i)
& - [_ols,i+ _ols,i, ]   (_ols,i)

Define \(\alpha_n^{(1)} = \frac{1}{\tilde{\tau}_n} (\frac{\lambda_n}{n} - \beta^\star_i)\), \(\alpha_n^{(2)} = \frac{1}{\tilde{\tau}_n} (\frac{\lambda_n}{n} + \beta^\star_i)\), \(\theta_n^{(1)} = \frac{1}{\tilde{\tau}_n}(\frac{\lambda_n + \delta_n}{n} - \beta^\star_i)\) and \(\theta_n^{(2)} = \frac{1}{\tilde{\tau}_n}(\frac{\lambda_n + \delta_n}{n} + \beta^\star_i)\).

In the order they appear, the four probabilities above are equal to \[\mathbb{P}(Z > \theta_n^{(1)} \mid \mathbf{X}) = 1 - \Phi(\theta_n^{(1)}),\] \[\mathbb{P}(Z < -\theta_n^{(2)} \mid \mathbf{X}) = \Phi(-\theta_n^{(2)}),\] \[\mathbb{P}(Z \in [\alpha_n^{(1)}, \theta_n^{(1)}] \mid \mathbf{X}) = \Phi(\theta_n^{(1)}) - \Phi(\alpha_n^{(1)}),\] \[\mathbb{P}(Z \in [-\theta_n^{(2)}, -\alpha_n^{(2)}] \mid \mathbf{X}) = \Phi(-\alpha_n^{(2)}) - \Phi(-\theta_n^{(2)}).\]

Using the first moment of the truncated normal [37], we have align [_ols,i- _ols,i, ] &= ^_i - + _n   [Z Z , ]
&= ^_i - - _n and align [_ols,i+ _ols,i, ] &= ^_i + + _n   [Z Z , ]
&= ^_i + - _n .

Therefore align &[__n + _n, i - __n, i ]
&= - (_ols,i> ) + (_ols,i< - )
& - [_ols,i- _ols,i, ]   (_ols,i)
& - [_ols,i+ _ols,i, ]   (_ols,i)
&= - (1 - (_n^(1))) + (-_n^(2))
& - (^_i - ) ((_n^(1)) - (_n^(1))) + _n ((_n^(1)) - (_n^(1)))
& - (^_i + ) ((-_n^(2)) - (-_n^(2))) + _n ((-_n^(2)) - (-_n^(2)))
&= - (1 - (_n^(1))) + (-_n^(2))
& - (^_i - ) (_n^(1) - _n^(1)) ‘(c_n^(1)) + _n (_n^(1) - _n^(1))’(d_n^(1))
& - (^_i + ) (_n^(2) - _n^(2)) ‘(-c_n^(2)) + _n (_n^(2) - _n^(2))’(-d_n^(2)) where \(c_n^{(1)}, d_n^{(1)} \in [\alpha_n^{(1)}, \theta_n^{(1)}]\) and \(c_n^{(2)}, d_n^{(2)} \in [\alpha_n^{(2)}, \theta_n^{(2)}]\) using first-order Taylor expansions.

We have \(\theta_n^{(1)} - \alpha_n^{(1)} = \theta_n^{(2)} - \alpha_n^{(2)} = \frac{1}{\tilde{\tau}_n} \frac{\delta_n}{n}\), \(\Phi' = \varphi\) and \(\varphi'(x) = - x \varphi(x)\), thus align [__n + _n, i - __n, i ] &= - (1 - (_n^(1))) + (-_n^(2))
& - (^_i - ) (c_n^(1)) - _n d_n^(1) (d_n^(1))
& - (^_i + ) (-c_n^(2)) - _n (-d_n^(2) (-d_n^(2)))
&= - (1 - (_n^(1))) + (-_n^(2))
& - (^_i - ) (c_n^(1)) - d_n^(1) (d_n^(1))
& - (^_i + ) (-c_n^(2)) - (-d_n^(2) (-d_n^(2)))
&= - (1 - (_n^(1))) + (-_n^(2))
& + _n^(1) (c_n^(1)) - d_n^(1) (d_n^(1))
& - _n^(2) (-c_n^(2)) - (-d_n^(2) (-d_n^(2))).

We first consider \(\beta^\star_i > 0\).

Since \(\lambda_n = o(n)\) and \(\delta_n = o(n)\), for \(n\) large enough, \(\frac{\lambda_n + \delta_n}{n} < \beta^\star_i\), so \(\alpha_n^{(1)} \leq \theta_n^{(1)} < 0\), thus for \(c_n^{(1)} \in [\alpha_n^{(1)}, \theta_n^{(1)}]\), we have \(|\alpha_n^{(1)} \varphi(c_n^{(1)})| \leq |\alpha_n^{(1)}| \varphi(\theta_n^{(1)}) = |\frac{\alpha_n^{(1)}}{\theta_n^{(1)}}| |\theta_n^{(1)}| \varphi(\theta_n^{(1)})\), where the ratio \(\frac{\alpha_n^{(1)}}{\theta_n^{(1)}} = \frac{\frac{\lambda_n}{n} - \beta^\star_i}{\frac{\lambda_n + \delta_n}{n} - \beta^\star_i}\) is deterministic and goes to \(1\).

As \(-\theta_n^{(2)} \leq -\alpha_n^{(2)} < 0\), for \(c_n^{(1)} \in [-\theta_n^{(2)}, -\alpha_n^{(2)}]\), we have \(|-\alpha_n^{(2)} \varphi(-c_n^{(2)})| \leq |-\alpha_n^{(2)} \varphi(-\alpha_n^{(2)})|\).

Since \(\frac{\mathbf{X}^\top \mathbf{X}}{n} \stackrel{\mathrm{a.s.}}{\longrightarrow}\mathbb{E}[X_0 X_0^\top]\) (strong law of large numbers), \(\lambda_n = o(n)\) and \(\delta_n = o(n)\), we have \(\tilde{\tau}_n \stackrel{\mathrm{a.s.}}{\longrightarrow}0^+\), and using the continuous mapping theorem, \(\alpha_n^{(1)} \stackrel{\mathrm{a.s.}}{\longrightarrow}-\infty\), \(\theta_n^{(1)} \stackrel{\mathrm{a.s.}}{\longrightarrow}-\infty\), \(\alpha_n^{(2)} \stackrel{\mathrm{a.s.}}{\longrightarrow}+\infty\) and \(\theta_n^{(2)} \stackrel{\mathrm{a.s.}}{\longrightarrow}+\infty\). We then also have \(d_n^{(1)} \stackrel{\mathrm{a.s.}}{\longrightarrow}-\infty\) and \(d_n^{(2)} \stackrel{\mathrm{a.s.}}{\longrightarrow}+\infty\).

\(\Phi\) and \(x \mapsto x \varphi(x)\) are continuous bounded functions so we get \(L^1\) convergence of \(\Phi(\theta_n^{(1)})\), \(\Phi(-\theta_n^{(2)})\), \(\theta_n^{(1)} \varphi(\theta_n^{(1)})\), \(-\alpha_n^{(2)} \varphi(-\alpha_n^{(2)})\), \(d_n^{(1)} \varphi(d_n^{(1)})\) and \(-d_n^{(2)} \varphi(-d_n^{(2)})\) to 0. By putting everything together, we obtain \[\frac{n}{\delta_n} \mathbb{E}[\hat{\beta}_{\lambda_n + \delta_n, i} - \hat{\beta}_{\lambda_n, i}] = \frac{n}{\delta_n} \mathbb{E}[\mathbb{E}[\hat{\beta}_{\lambda_n + \delta_n, i} - \hat{\beta}_{\lambda_n, i} \mid \mathbf{X}]] \rightarrow -1 = -\text{sign}(\beta^\star_i).\]

When \(\beta^\star_i < 0\), we show in a similar manner that \[\frac{n}{\delta_n} \mathbb{E}[\hat{\beta}_{\lambda_n + \delta_n, i} - \hat{\beta}_{\lambda_n, i}] \rightarrow 1 = -\text{sign}(\beta^\star_i).\]

If \(\beta^\star_i = 0\), \(\alpha_n^{(1)} = \alpha_n^{(2)}\) and \(\theta_n^{(1)} = \theta_n^{(2)}\) so \(1 - \Phi(\alpha_n^{(1)}) = \Phi(-\alpha_n^{(2)})\), \(\varphi(\alpha_n^{(1)}) = \varphi(-\alpha_n^{(2)})\), \(1 - \Phi(\theta_n^{(1)}) = \Phi(-\theta_n^{(2)})\) and \(\varphi(\theta_n^{(1)}) = \varphi(-\theta_n^{(2)})\) which leads to align &[__n + _n, i - __n, i ]
&= - (1 - (_n^(1))) + (-_n^(2))
& - (^_i - ) ((_n^(1)) - (_n^(1))) + _n ((_n^(1)) - (_n^(1)))
& - (^_i + ) ((-_n^(2)) - (-_n^(2))) + _n ((-_n^(2)) - (-_n^(2)))
&= 0 and thus \(\mathbb{E}[\hat{\beta}_{\lambda_n + \delta_n, i} - \hat{\beta}_{\lambda_n, i}] = 0 = \text{sign}(\beta^\star_i)\).

Thus, we have convergence component-wise and can conclude \(\frac{n}{\delta_n} \mathbb{E}[\hat{\beta}_{\lambda_n + \delta_n} - \hat{\beta}_{\lambda_n}] \rightarrow -\text{sign}(\beta^\star)\).

We now show that \(\frac{n}{\delta_n} \mathbb{E}[\hat{\beta}_{\lambda_n + \delta_n, i} \hat{\beta}_{\lambda_n + \delta_n, j} - \hat{\beta}_{\lambda_n, i} \hat{\beta}_{\lambda_n, j}] \rightarrow - (\text{sign}(\beta^\star_i) \beta^\star_j + \beta^\star_i \text{sign}(\beta^\star_j))\).

Note that align & [ (__n + _n, i __n + _n, j - __n, i __n, j) + (^_i) ^_j + ^_i (^_j)]
&= [ (__n + _n, i - __n, i) __n + _n, j + (^_i) ^_j] + [__n, i   (__n + _n, j - __n, j) + ^_i (^_j)] with align &[ (__n + _n, i - __n, i) __n + _n, j + (^_i) ^_j]
&= [( (__n + _n, i - __n, i) + (^_i)) (__n + _n, j - ^_j)]
&+ ^_j   [ (__n + _n, i - __n, i) + (^_i)] - (^_i) [__n + _n, j - ^_j] where, using Cauchy–Schwarz, align &[( (__n + _n, i - __n, i) + (^_i)) (__n + _n, j - ^_j)]
&.

We can do the same with \(\mathbb{E}[\hat{\beta}_{\lambda_n, i} \, \frac{n}{\delta_n} (\hat{\beta}_{\lambda_n + \delta_n, j} - \hat{\beta}_{\lambda_n, j}) + \beta^\star_i \text{sign}(\beta^\star_j)]\).

Therefore, proving \(\mathbb{E}[\frac{n}{\delta_n} (\hat{\beta}_{\lambda_n + \delta_n, i} \hat{\beta}_{\lambda_n + \delta_n, j} - \hat{\beta}_{\lambda_n, i} \hat{\beta}_{\lambda_n, j}) + \text{sign}(\beta^\star_i) \beta^\star_j + \beta^\star_i \text{sign}(\beta^\star_j)] \rightarrow 0\) for all \(i, j\) comes down to proving \(\mathbb{E}[(\frac{n}{\delta_n} (\hat{\beta}_{\lambda_n + \delta_n, i} - \hat{\beta}_{\lambda_n, i}) + \text{sign}(\beta^\star_i))^2] = O(1)\) for all \(i\) given that we have already shown for all \(i\), accounting for the fact that both \(\lambda_n\) and \(\delta_n\) are \(o(n)\),

• \(\mathbb{E}[\hat{\beta}_{\lambda_n, i}] \rightarrow \beta^\star_i\) and \(\mathbb{E}[\hat{\beta}_{\lambda_n + \delta_n, i}] \rightarrow \beta^\star_i\),

• \(\mathbb{E}[(\hat{\beta}_{\lambda_n, i} - \beta^\star_i)^2] \rightarrow 0\) and \(\mathbb{E}[(\hat{\beta}_{\lambda_n + \delta_n, i} - \beta^\star_i)^2] \rightarrow 0\),

• \(\frac{n}{\delta_n} \mathbb{E}[\hat{\beta}_{\lambda_n + \delta_n, i} - \hat{\beta}_{\lambda_n, i}] \rightarrow - \text{sign}(\beta^\star_i)\).

The first two bullet points were proved in 14 and the third one earlier in this proof.

As a reminder, we have align __n + _n, i - __n, i &= -(_ols,i)

thus align (__n + _n, i - __n, i)^2 and align ( (__n + _n, i - __n, i) + (^_i))^2 ( (__n + _n, i - __n, i)^2 + (^_i)^2) . Hence, \(\mathbb{E}[(\frac{n}{\delta_n} (\hat{\beta}_{\lambda_n + \delta_n, i} - \hat{\beta}_{\lambda_n, i}) + \text{sign}(\beta^\star_i))^2] = O(1)\).

Therefore, we get \[\frac{n}{\delta_n} \mathbb{E}[\hat{\beta}_{\lambda_n + \delta_n} \hat{\beta}_{\lambda_n + \delta_n}^\top - \hat{\beta}_{\lambda_n} \hat{\beta}_{\lambda_n}^\top] \rightarrow -(\text{sign}(\beta^\star) \beta^{\star\top}+ \beta^\star\text{sign}(\beta^\star)^\top).\]

We can then conclude that \(\frac{n^2}{\delta_n^2} \, \sigma^2(h_n^{\mathrm{diff}}) \rightarrow 4\tau^2 \|\beta^\star\|_0\) as mentioned earlier in the proof.

12 Proof of 7:↩︎

Starting from the expression for \(\gamma(h_n)\) stated in 2, we have align (h_n^) &= [(2 (Y_0 X_0^- ^[X_0 X_0^]) _n + \(\operatorname{tr}\)((X_0 X_0^- [X_0 X_0^]) _n))^2]. where

• \(\nu_n \triangleq\hat{\beta}_{\lambda_n + \delta_n} - \hat{\beta}_{\lambda_n} - (\hat{\beta}'_{\lambda_n + \delta_n} - \hat{\beta}'_{\lambda_n})\),

• \(\Psi_n \triangleq\hat{\beta}_{\lambda_n} \hat{\beta}_{\lambda_n}^\top - \hat{\beta}_{\lambda_n + \delta_n} \hat{\beta}_{\lambda_n + \delta_n}^\top - (\hat{\beta}'_{\lambda_n} \hat{\beta}_{\lambda_n}^{'\top} - \hat{\beta}'_{\lambda_n + \delta_n} \hat{\beta}_{\lambda_n + \delta_n}^{'\top})\).

\(\mathbb{E}[X_0 X_0^\top] = \mathbf{I}\) since the features are drawn from \(\mathcal{N}(0, \mathbf{I})\), and using independence of \(Z_0\) from the training points, we have align (h_n^) &= [(2 _i (Y_0 X_0, i - ^_i) _n, i + _i, j (X_0, i X_0, j - ) _n, i, j))^2]
&= 4 _i [(Y_0 X_0, i - ^_i)^2] [_n, i^2]
& + 4 _i j [(Y_0 X_0, i - ^_i) (Y_0 X_0, j - ^_j)] [_n, i _n, j]
& + 4 _i, j, k [(Y_0 X_0, i - ^_i) (X_0, j X_0, k - )] [_n, i _n, j, k]
& + _i, j, k, l [(X_0, i X_0, j - ) (X_0, k X_0, l - )] [_n, i, j _n, k, l].

Since \(Y_0 = X_0^\top \beta^\star+ \varepsilon_0 = \sum_k X_{0, k} \beta^\star_k + \varepsilon_0\) with \(X_0 \sim \mathcal{N}(0, \mathbf{I})\) and \(\varepsilon_0 \mathpalette{\independenT}{\perp}X_0\), we have align [Y_0 X_0, i] = ^_i [X_0, i^2] + _k i ^_k [X_0, i X_0, k] + [_0 X_0, i] = ^_i and \(Y_0^2 = \sum_{k, l} X_{0, k} X_{0, l} \beta^\star_k \beta^\star_l + 2 \varepsilon_0 \sum_k X_{0, k} \beta^\star_k + \varepsilon_0^2\), so for \(i \neq j\), align [Y_0^2 X_0, i X_0, j] = _k, l [X_0, i X_0, j X_0, k X_0, l] ^_k ^_l + 2 _k [_0 X_0, i X_0, j X_0, k] ^_k + [_0^2 X_0, i X_0, j] = 2 ^_i ^_j since the expectation in the first sum is equal to \(1\) when \(k = i, l = j\) or \(k = j, l = i\), and equal to \(0\) otherwise, and thus, for \(i \neq j\), align [(Y_0 X_0, i - ^_i) (Y_0 X_0, j - ^_j)] = [Y_0^2 X_0, i X_0, j] - ^_i [Y_0 X_0, j] - ^_j [Y_0 X_0, i] + ^_i ^_j = ^_i ^_j. For the case \(i = j\), align [Y_0^2 X_0, i^2] &= _k, l [X_0, i^2 X_0, k X_0, l] ^_k ^_l + 2 _k [_0 X_0, i^2 X_0, k] ^_k + [_0^2 X_0, i^2]
&= [X_0, i^4] ^_i + _k i [X_0, i^2 X_0, k^2] ^_k + ^2
&= [X_0, i^4] ^_i + _k i ^_k + ^2 and then, for \(\beta^\star_i = 0\), align [(Y_0 X_0, i - ^_i)^2] = [Y_0^2 X_0, i^2] = _k i ^_k + ^2 ^2 > 0.

Therefore align (h_n^) &= 4 _i, ^_i = 0 [Y_0^2 X_0, i^2] [_n, i^2]
& + 4 _i, ^_i [(Y_0 X_0, i - ^_i)^2] [_n, i^2]
& + 4 _i j, ^_i , ^_j ^_i ^_j [_n, i _n, j]
& + 4 _i, j, k [(Y_0 X_0, i - ^_i) (X_0, j X_0, k - )] [_n, i _n, j, k]
& + _i, j, k, l [(X_0, i X_0, j - ) (X_0, k X_0, l - )] [_n, i, j _n, k, l]. where importantly we were able to remove the \(i, j\) terms in the third sum when \(\beta^\star_i = 0\) or \(\beta^\star_j = 0\).

We will now prove the following results:

• \(\mathbb{E}[\nu_{n, i}^2] = O(\frac{\delta_n^2}{n^2})\) for all \(i\),

• \(\mathbb{E}[\nu_{n, i}^2] = \Omega(\frac{\delta_n^2}{n^2 \sqrt{n}})\) for all \(i\) such that \(\beta^\star_i = 0\),

• \(\mathbb{E}[\nu_{n, i}^2] = o(\frac{\delta_n^2}{n^2 \sqrt{n}})\) for all \(i\) such that \(\beta^\star_i \neq 0\),

• \(\mathbb{E}[\Psi_{n, i, j}^2] = O(\frac{\delta_n^2}{n^4})\) for all \(i, j\).

Once we prove these, Cauchy–Schwarz will yield the following upper-bounding rates for terms appearing in the expression of \(\gamma(h_n^{\mathrm{diff}})\):

• for \(i, j\) such that \(\beta^\star_i \neq 0\) and \(\beta^\star_j \neq 0\), \(|\mathbb{E}[\nu_{n, i} \nu_{n, j}]| \leq \sqrt{\mathbb{E}[\nu_{n, i}^2] \mathbb{E}[\nu_{n, j}^2]} = o(\frac{\delta_n^2}{n^2 \sqrt{n}})\),

• \(|\mathbb{E}[\nu_{n, i} \Psi_{n, j, k}]| \leq \sqrt{\mathbb{E}[\nu_{n, i}^2] \mathbb{E}[\Psi_{n, j, k}^2]} = O(\sqrt{\frac{\delta_n^2}{n^2} \frac{\delta_n^2}{n^4}}) = O(\frac{\delta_n^2}{n^3}) = o(\frac{\delta_n^2}{n^2 \sqrt{n}})\),

• \(|\mathbb{E}[\Psi_{n, i, j} \Psi_{n, k, l}]| \leq \sqrt{\mathbb{E}[\Psi_{n, i, j}^2] \mathbb{E}[\Psi_{n, k, l}^2]} = O(\sqrt{\frac{\delta_n^2}{n^4} \frac{\delta_n^2}{n^4}}) = O(\frac{\delta_n^2}{n^4}) = o(\frac{\delta_n^2}{n^2 \sqrt{n}})\),

and it will therefore be clear that \(\gamma(h_n^{\mathrm{diff}}) = \Omega(\frac{\delta_n^2}{n^2 \sqrt{n}})\) as the terms of leading order in \(\gamma(h_n^{\mathrm{diff}})\) will be the \(\mathbb{E}[\nu_{n, i}^2]\) terms for \(i\) such that \(\beta^\star_i = 0\).

We will now prove the first result \(\mathbb{E}[\nu_{n, i}^2] = O(\frac{\delta_n^2}{n^2})\) for all \(i\).

We have align _n, i &= __n + _n, i - __n, i - (‘__n + _n, i -’__n, i)
&= (_ols,i) (|_ols,i| - )_+ - (_ols,i) (|_ols,i| - )_+
&- ((_ols,i’) (|_ols,i’| - )_+ - (_ols,i’) (|_ols,i’| - )_+)
&= (_ols,i)

&- (_ols,i’)

. We can observe that both \(|\hat{\beta}_{\lambda_n + \delta_n, i} - \hat{\beta}_{\lambda_n, i}|\) and \(|\hat{\beta}'_{\lambda_n + \delta_n, i} - \hat{\beta}'_{\lambda_n, i}|\) are upper-bounded by \(\frac{\delta_n}{n}\) and thus \(\nu_{n, i}^2 \leq 4 \frac{\delta_n^2}{n^2}\), which implies \(\mathbb{E}[\nu_{n, i}^2] = O(\frac{\delta_n^2}{n^2})\) for all \(i\).

We will now prove the second result \(\mathbb{E}[\nu_{n, i}^2] = \Omega(\frac{\delta_n^2}{n^2 \sqrt{n}})\) for all \(i\) such that \(\beta^\star_i = 0\).

Based on the previous expression, we can further detail \(\nu_{n, i}\) as follows align _n, i &=

-

which means there are \(25\) possible cases that form a partition and we can write \(\nu_{n, i}\) as the sum of \(25\) terms that are of the form: an indicator of one of the \(25\) events multiplied by the value of \(\nu_{n, i}\) for this event. We can then similarly write \(\nu_{n, i}^2\) as the sum of \(25\) terms that are of the form: an indicator of one of the \(25\) events multiplied by the value of \(\nu_{n, i}^2\) for this event.

We can then lower-bound \(\mathbb{E}[\nu_{n, i}^2]\) by the expectation of any one of the \(25\) terms since they are all non-negative. In particular, we can do it using the term coming from the combination of the first case on the left side and the last case on the right side align [_n, i^2] &[ ]
&= (_ols,i> , |_ols,i’| < ).

Since \(\lambda_n = \omega(1)\) and \(\delta_n = \Theta(1)\), \(\frac{\lambda_n}{n} - \frac{\delta_n}{n} > 0\) for \(n\) large enough, and we then have \(\{\hat{\beta}_{\mathrm{\small ols},i}> \hat{\beta}_{\mathrm{\small ols},i}' + 2 \frac{\delta_n}{n}, \hat{\beta}_{\mathrm{\small ols},i}' \in [\frac{\lambda_n}{n} - \frac{\delta_n}{n}, \frac{\lambda_n}{n}]\} \subseteq \{\hat{\beta}_{\mathrm{\small ols},i}> \frac{\lambda_n+\delta_n}{n}, |\hat{\beta}_{\mathrm{\small ols},i}'| < \frac{\lambda_n}{n}\}\), therefore align &(_ols,i> , |_ols,i’| < )
&(_ols,i> _ols,i’ + 2 , _ols,i’ )
&= (n (_ols,i’ - _ols,i) < - 2 _n, _ols,i’ ). We have align (_ols,i, _ols,i’ , ) &= (^+ (^)^-1 ^, ^+ (^)^-1 ^’ , )
&= (^)^-1 ^(, ‘) (^)^-1
&= ^2 (^)^-1 ^ (^)^-1 where \(\mathbf{\tilde{X}} \triangleq(X_2, \ldots, X_n)^\top\) is the matrix of regressors for the training points except for the first one that is being changed, since \(\mathrm{Cov}(\varepsilon_i, \varepsilon'_j)\) is equal to \(\tau^2\) if \(i = j \geq 2\) and \(0\) otherwise. Then align (_ols,i’ - _ols,i, _ols,i’ , ) &= ^2 (^)^-1 - ^2 (^)^-1 ^ (^)^-1
&= ^2 ( - (^)^-1 ^) (^)^-1. Hence, the bivariate normal vector \((\hat{\beta}_{\mathrm{\small ols},i}' - \hat{\beta}_{\mathrm{\small ols},i}, \hat{\beta}_{\mathrm{\small ols},i}')\) has uncorrelated components in the limit, with zero correlation being equivalent to independence for multivariate normal vectors. Since \(n (\hat{\beta}'_{\text{OLS}} - \hat{\beta}_{\mathrm{\small ols}}) \stackrel{\mathrm{a.s.}}{\longrightarrow}V \triangleq(Y'_1 - X_1^{\prime\top} \beta^\star) X'_1 - (Y_1 - X_1^\top \beta^\star) X_1\), proved in 15, and \(\delta_n = \Theta(1)\), we have align (n (_ols,i’ - _ols,i) < - 2 _n, _ols,i’ ) = ((_ols,i’ )). We can then focus on the rate of \(\mathbb{P}(\hat{\beta}_{\mathrm{\small ols},i}' \in [\frac{\lambda_n}{n} - \frac{\delta_n}{n}, \frac{\lambda_n}{n}])\). align (_ols,i’ ) = [(_ols,i’ )] where, using \(\beta^\star_i = 0\) and \(\hat{\beta}'_{\text{OLS}} \mid \mathbf{X} \sim \mathcal{N}(\beta^\star, \tau^2 (\mathbf{X'}^\top \mathbf{X'})^{-1})\) and defining \(\tilde{\tau}'_n = \frac{\tau}{\sqrt{n}} \sqrt{(\frac{\mathbf{X'}^\top \mathbf{X'}}{n})_{i,i}^{-1}}\), align (_ols,i’ ) &= ( ) - ( ( - ))
&= ( ) - ’(c_n) for \(c_n \in [\frac{1}{\tilde{\tau}'_n} (\frac{\lambda_n}{n} - \frac{\delta_n}{n}), \frac{1}{\tilde{\tau}'_n} \frac{\lambda_n}{n}]\) by a second-order Taylor expansion. We have align ( ) = ( ) whose expectation is \(\Theta(\frac{1}{\sqrt{n}})\) since \(\lambda_n = O(\sqrt{n})\) and \(\delta_n = \Theta(1)\) yield \(\frac{\delta_n}{\sqrt{n}} = \Theta(\frac{1}{\sqrt{n}})\) and \(\mathbb{E}[\frac{1}{\tilde{\tau}'_n} \frac{\sqrt{n}}{n} \varphi(\frac{1}{\tilde{\tau}'_n} \frac{\lambda_n}{n})] = \Theta(1)\).

As for the second part of the Taylor expansion, its expectation is a \(o(\frac{1}{\sqrt{n}})\) since \(\varphi'\) is bounded, \(\delta_n = \Theta(1)\) and we have align [] = [((^)_i,i^-1)^-1] = n - p + 1 using the fact that for \(X_i \overset{\textrm{\tiny\textrm{i.i.d.}}}{\sim}\mathcal{N}(0, \mathbf{I})\), we know \((\mathbf{X}^\top \mathbf{X})^{-1} \sim W_p^{-1}(\mathbf{I}, n)\) and then the diagonal element \((\mathbf{X}^\top \mathbf{X})_{i, i}^{-1}\) follows an inverse gamma distribution with shape parameter \(\frac{n-p+1}{2}\) and scale parameter \(\frac{1}{2}\), and the expectation of the reciprocal of an inverse gamma distributed variable is the ratio of the shape and the scale.

We can then conclude that align (_ols,i’ ) = [(_ols,i’ )] = () and thus \(\mathbb{E}[\nu_{n, i}^2] = \Omega(\frac{\delta_n^2}{n^2 \sqrt{n}})\).

We will now prove the third result \(\mathbb{E}[\nu_{n, i}^2] = o(\frac{\delta_n^2}{n^2 \sqrt{n}})\) for all \(i\) such that \(\beta^\star_i \neq 0\).

Consider \(i\) such that \(\beta^\star_i > 0\), since the combination of the first case on the left side and the first case on the right side in the expression of \(\nu_{n, i}\) corresponds to a value of \(0\) for \(\nu_{n, i}\), we can write \(\nu_{n, i}^2\) as the sum of \(24\) terms that are of the form: an indicator of one of the \(24\) other events multiplied by the value of \(\nu_{n, i}^2\) for this event. Since \(\nu_{n, i}^2 \leq 4 \frac{\delta_n^2}{n^2}\), we can upper-bound \(\nu_{n, i}^2\) by \(4 \frac{\delta_n^2}{n^2}\) multiplied by the sum of the \(24\) indicators and we then need to show that all \(24\) indicators have an expectation which is \(o(\frac{1}{\sqrt{n}})\). Using \(\mathbb{E}[\mathbb{1}\mathopen{}\mathclose{\left[{A}}\right]] = \mathbb{P}(A)\), \(\mathbb{P}(A \cap B) \leq \min(\mathbb{P}(A), \mathbb{P}(B))\) and the fact that \(\hat{\beta}_{\mathrm{\small ols}}\) and \(\hat{\beta}'_{\text{OLS}}\) have the same unconditional distribution, we can upper-bound all \(24\) indicator expectations by one of the following four probabilities

• \(\mathbb{P}(\hat{\beta}_{\mathrm{\small ols},i}< -\frac{\lambda_n+\delta_n}{n}) = \mathbb{E}[\mathbb{P}(\hat{\beta}_{\mathrm{\small ols},i}< -\frac{\lambda_n+\delta_n}{n} \mid \mathbf{X})]\),

• \(\mathbb{P}(\hat{\beta}_{\mathrm{\small ols},i}\in [\frac{\lambda_n}{n}, \frac{\lambda_n+\delta_n}{n}]) = \mathbb{E}[\mathbb{P}(\hat{\beta}_{\mathrm{\small ols},i}\in [\frac{\lambda_n}{n}, \frac{\lambda_n+\delta_n}{n}] \mid \mathbf{X})]\),

• \(\mathbb{P}(\hat{\beta}_{\mathrm{\small ols},i}\in [-\frac{\lambda_n+\delta_n}{n}, -\frac{\lambda_n}{n}]) = \mathbb{E}[\mathbb{P}(\hat{\beta}_{\mathrm{\small ols},i}\in [-\frac{\lambda_n+\delta_n}{n}, -\frac{\lambda_n}{n}] \mid \mathbf{X})]\),

• \(\mathbb{P}(|\hat{\beta}_{\mathrm{\small ols},i}| < \frac{\lambda_n}{n}) = \mathbb{E}[\mathbb{P}(|\hat{\beta}_{\mathrm{\small ols},i}| < \frac{\lambda_n}{n} \mid \mathbf{X})]\).

Since \(\hat{\beta}_{\mathrm{\small ols}}\mid \mathbf{X} \sim \mathcal{N}(\beta^\star, \tau^2 (\mathbf{X}^\top \mathbf{X})^{-1})\), we can write \(\hat{\beta}_{\mathrm{\small ols},i}= \beta^\star_i + \tilde{\tau}_n Z\) where \(\tilde{\tau}_n = \frac{\tau}{\sqrt{n}} \sqrt{(\frac{\mathbf{X}^\top \mathbf{X}}{n})_{i,i}^{-1}}\) and \(Z \mid \mathbf{X} \sim \mathcal{N}(0, 1)\). Note that we could have \(i\) as a subscript of \(\tilde{\tau}_n\) and \(Z\), but we will only consider one \(i\) at a time in our computations and we can thus omit this subscript for both of them for the sake of notational simplicity, and we will also omit it for some additional notation we define in the rest of the proof.

Define \(\alpha_n^{(1)} = \frac{1}{\tilde{\tau}_n} (\frac{\lambda_n}{n} - \beta^\star_i)\), \(\alpha_n^{(2)} = \frac{1}{\tilde{\tau}_n} (\frac{\lambda_n}{n} + \beta^\star_i)\), \(\theta_n^{(1)} = \frac{1}{\tilde{\tau}_n}(\frac{\lambda_n + \delta_n}{n} - \beta^\star_i)\) and \(\theta_n^{(2)} = \frac{1}{\tilde{\tau}_n}(\frac{\lambda_n + \delta_n}{n} + \beta^\star_i)\).

In the order they appear, the four conditional probabilities above are equal to \[\mathbb{P}(Z < -\theta_n^{(2)} \mid \mathbf{X}) = \Phi(-\theta_n^{(2)}),\] \[\mathbb{P}(Z \in [\alpha_n^{(1)}, \theta_n^{(1)}] \mid \mathbf{X}) = \Phi(\theta_n^{(1)}) - \Phi(\alpha_n^{(1)}),\] \[\mathbb{P}(Z \in [-\theta_n^{(2)}, -\alpha_n^{(2)}] \mid \mathbf{X}) = \Phi(-\alpha_n^{(2)}) - \Phi(-\theta_n^{(2)}),\] \[\mathbb{P}(Z \in [-\alpha_n^{(2)}, \alpha_n^{(1)}] \mid \mathbf{X}) = \Phi(\alpha_n^{(1)}) - \Phi(-\alpha_n^{(2)}).\]

Since \(\frac{\mathbf{X}^\top \mathbf{X}}{n} \stackrel{\mathrm{a.s.}}{\longrightarrow}\mathbb{E}[X_0 X_0^\top]\) (strong law of large numbers), \(\lambda_n = o(n)\) and \(\delta_n = o(n)\), we have \(\tilde{\tau}_n \stackrel{\mathrm{a.s.}}{\longrightarrow}0^+\), and using the continuous mapping theorem, \(\alpha_n^{(1)} \stackrel{\mathrm{a.s.}}{\longrightarrow}-\infty\), \(\theta_n^{(1)} \stackrel{\mathrm{a.s.}}{\longrightarrow}-\infty\), \(\alpha_n^{(2)} \stackrel{\mathrm{a.s.}}{\longrightarrow}+\infty\) and \(\theta_n^{(2)} \stackrel{\mathrm{a.s.}}{\longrightarrow}+\infty\) as \(\beta^\star_i > 0\).

If we show that \(\sqrt{n} \, \Phi(\alpha_n^{(1)})\) goes to \(0\) in \(L^1\), then all other similar convergences will follow and we will get that all four unconditional probabilities listed above are \(o(\frac{1}{\sqrt{n}})\) and thus \(\mathbb{E}[\nu_{n, i}^2] = o(\frac{\delta_n^2}{n^2 \sqrt{n}})\) for all \(i\) such that \(\beta^\star_i > 0\).

We have align   (_n^(1)) &= _n^(1) (_n^(1)) thus, by Cauchy–Schwarz, \[\mathbb{E}[|\sqrt{n} \, \Phi(\alpha_n^{(1)})|] \leq \sqrt{\mathbb{E}\mathopen{}\mathclose{\left[\frac{n}{(\alpha_n^{(1)})^2}}\right] \mathbb{E}\mathopen{}\mathclose{\left[(\alpha_n^{(1)} \Phi(\alpha_n^{(1)}))^2}\right]}.\]

\(\alpha_n^{(1)} \stackrel{\mathrm{a.s.}}{\longrightarrow}-\infty\) so \(\alpha_n^{(1)} \Phi(\alpha_n^{(1)}) \stackrel{\mathrm{a.s.}}{\longrightarrow}0\). This comes from the fact that \(x \, \Phi(x) \rightarrow 0\) for \(x \rightarrow -\infty\), as we notice that for \(x < 0\), we have \(0 < -x \, \Phi(x) = -x \, (1 - \Phi(-x)) = -x \int_{-x}^{+\infty} \varphi(t) dt \leq \int_{-x}^{+\infty} t \varphi(t) dt\) where this last expression goes to \(0\) when \(x \rightarrow -\infty\).

Since \(\lambda_n = o(n)\) and \(\delta_n = o(n)\), for \(n\) large enough, \(\frac{\lambda_n + \delta_n}{n} < \beta^\star_i\), so \(\alpha_n^{(1)} \leq \theta_n^{(1)} < 0\). Since the function \(x \mapsto x \, \Phi(x)\) is continuous bounded for \(x < 0\), we get \(L^1\) convergence of \((\alpha_n^{(1)} \Phi(\alpha_n^{(1)}))^2\) to 0.

Moreover, \(\frac{n}{(\alpha_n^{(1)})^2} = \frac{n \tilde{\tau}_n^2}{(\frac{\lambda_n}{n} - \beta^\star_i)^2} = \frac{\tau^2}{(\frac{\lambda_n}{n} - \beta^\star_i)^2} (\frac{\mathbf{X}^\top \mathbf{X}}{n})_{i,i}^{-1}\) and it is thus sufficient to have \(\mathbb{E}[(\frac{\mathbf{X}^\top \mathbf{X}}{n})_{i,i}^{-1}] = O(1)\), which is the case for features drawn \(\textrm{i.i.d.}\) from \(\mathcal{N}(0, \mathbf{I})\) as \(\mathbb{E}[(\frac{\mathbf{X}^\top \mathbf{X}}{n})_{i,i}^{-1}] = \frac{n}{n-p-1}\).

Hence, \(\sqrt{n} \, \Phi(\alpha_n^{(1)})\) goes to \(0\) in \(L^1\).

The proof is similar for \(i\) such that \(\beta^\star_i < 0\).

Finally, we show the fourth result \(\mathbb{E}[\Psi_{n, i, j}^2] = O(\frac{\delta_n^2}{n^4})\) for all \(i, j\) or equivalently \(\mathbb{E}[\Psi_{n, i, j}^2] = O(\frac{\delta_n^4}{n^4})\) since \(\delta_n = \Theta(1)\).

Similarly to previous computations and upper-bounding with Cauchy–Schwarz, we can upper-bound \(\mathbb{E}[\Psi_{n, i, j}^2]\) using products of \(\mathbb{E}[\nu_{n, i}^4]\) and the fourth moment of \(\hat{\beta}_{\lambda_n, i}\), \(\hat{\beta}_{\lambda_n + \delta_n, i}\), \(\hat{\beta}'_{\lambda_n, i}\) or \(\hat{\beta}'_{\lambda_n + \delta_n, i}\) and their counterparts for \(j\). Since \(\nu_{n, i}^2 \leq 4 \frac{\delta_n^2}{n^2}\), we have \(\mathbb{E}[\nu_{n, i}^4] = O(\frac{\delta_n^4}{n^4})\). Additionally, the fourth moments are bounded as we showed the \(L^4\) consistency of soft-thresholding for \(\beta^\star\). This gives us \(\mathbb{E}[\Psi_{n, i, j}^2] = O(\frac{\delta_n^4}{n^4})\).

With the four results proved, we can conclude that \(\gamma(h_n^{\mathrm{diff}}) = \Omega(\frac{\delta_n^2}{n^2 \sqrt{n}})\).

13 Proof of 2:↩︎

2 follows immediately from the following two propositions, proved in [sec:rate-sigma-single-algo-stlasso] [sec:rate-lstab-single-algo-stlasso], respectively.

14 Proof of 8:↩︎

We start by introducing a lemma which provides equations that will prove useful in the single algorithm setting.

Proof align h_n(Z_0, ) &= (Y_0 - X_0^)^2
&= Y_0^2 - 2 Y_0 X_0^+ (X_0^)^2
&= Y_0^2 - 2 Y_0 X_0^+ X_0^^X_0
&= Y_0^2 - 2 Y_0 X_0^+ \(\operatorname{tr}\)(X_0 X_0^^)

Note that \(\mathbb{E}[Y_0 X_0^\top] = \mathbb{E}[\mathbb{E}[Y_0 \mid X_0] X_0^\top] = \mathbb{E}[X_0^\top \beta^\star X_0^\top] = \beta^{\star\top}\mathbb{E}[X_0 X_0^\top]\).

Since \(\hat{\beta}\) is only a function of \(\mathbf{Z}\), the independence of \(Z_0\) and \(\mathbf{Z}\) yields the next three equations.

The fifth equation comes from noticing align ^2(h_n) = ([h_n(Z_0, ) Z_0]) = [([h_n(Z_0, ) Z_0] - [h_n(Z_0, )])^2]. And the last one comes from the definition of \(\gamma(h_n)\) as align (h_n) = [(h(Z_0, ) - h(Z_0, ) - ([h(Z_0, ) ] - [h(Z_0, ) ]))^2]. 0◻

In addition to giving a first glimpse into the differences between the single algorithm and comparison settings, the following lemma plays an important role in our proof via its result for a single linear predictor.

Proof Let \(Y_0 = X_0^\top \beta^\star+ \varepsilon_0\) be the response variable with \(\mathrm{Var}(\varepsilon_0) = \tau^2\). Using the information on the distribution of \(X_0\) and the independence of \(X_0\) and \(\varepsilon_0\), note that align [Y_0^2] = (Y_0) + [Y_0]^2 = (X_0^^+ _0) + 0 = ^(X_0) ^+ (_0) = ^^+ ^2. For the single linear predictor, starting from the expression of \(\sigma^2(h_n)\) in 3, since \(\mathbb{E}[\hat{\beta}_n]\) and \(\mathbb{E}[\hat{\beta}_n \hat{\beta}_n^\top]\) are non-random, we can expand the square, use linearity of expectation, take the limits and factorize back to obtain the convergence align ^2(h_n^) &= [(Y_0^2 - [Y_0^2] - 2 (Y_0 X_0^- ^[X_0 X_0^]) [_n] + \(\operatorname{tr}\)((X_0 X_0^- [X_0 X_0^]) [_n _n^]))^2]
&[(Y_0^2 - [Y_0^2] - 2 (Y_0 X_0^- ^[X_0 X_0^]) ^+ \(\operatorname{tr}\)((X_0 X_0^- [X_0 X_0^]) ^^))^2]
& = [(Y_0^2 - ^^- ^2 - 2 Y_0 X_0^^+ 2^^+ \(\operatorname{tr}\)((X_0 X_0^^^- ^^))^2]
& = [((X_0^^+ _0)^2 - ^^- ^2 - 2 (X_0^^+ _0) X_0^^+ 2 ^^+ (X_0^^)^2 - ^^)^2]
& = [(_0^2 - ^2)^2] = (_0^2) = [_0^4] - [_0^2]^2 = 3^4 - ^4 = 2^4. Similarly, we derive the second result with two linear predictors by starting from the expression of \(\sigma^2(h_n)\) in 2. 0◻

We will show that \(\mathbb{E}[\hat{\beta}_{\lambda_n}] \rightarrow \beta^\star\) and \(\mathbb{E}[\hat{\beta}_{\lambda_n} \hat{\beta}_{\lambda_n}^\top] \rightarrow \beta^\star\beta^{\star\top}\) in order to obtain the convergence of \(\sigma^2(h_n^{\mathrm{sing}})\) as an application of 4.

We have for \(i=1, \ldots, p,\) align __n, i &= (_ols,i) (|_ols,i| - )_+
&= (_ols,i)

.

A classic result for the OLS estimator is \(\hat{\beta}_{\mathrm{\small ols}}\mid \mathbf{X} \sim \mathcal{N}(\beta^\star, \tau^2 (\mathbf{X}^\top \mathbf{X})^{-1})\). We can write \(\hat{\beta}_{\mathrm{\small ols},i}= \beta^\star_i + \tilde{\tau}_n Z\) where \(\tilde{\tau}_n = \frac{\tau}{\sqrt{n}} \sqrt{(\frac{\mathbf{X}^\top \mathbf{X}}{n})_{i,i}^{-1}}\) and \(Z \mid \mathbf{X} \sim \mathcal{N}(0, 1)\). Note that we could have \(i\) as a subscript of \(\tilde{\tau}_n\) and \(Z\), but we will only consider one \(i\) at a time in our computations and we can thus omit this subscript for both of them for the sake of notational simplicity, and we will also omit it for some additional notation we define in the rest of the proof.

We now show that \(\mathbb{E}[\hat{\beta}_{\lambda_n, i}] \rightarrow \beta^\star_i\).

Using the law of total expectation, align [__n, i ] &= [_ols,i- _ols,i, ]   (_ols,i )
& + [_ols,i+ _ols,i-, ]   (_ols,i- ).

Define \(\alpha_n^{(1)} = \frac{1}{\tilde{\tau}_n} (\frac{\lambda_n}{n} - \beta^\star_i)\) and \(\alpha_n^{(2)} = \frac{1}{\tilde{\tau}_n} (\frac{\lambda_n}{n} + \beta^\star_i)\).

The first probability is equal to \[\mathbb{P}(Z \geq \alpha_n^{(1)} \mid \mathbf{X}) = 1 - \Phi(\alpha_n^{(1)})\]

and the second probability to \[\mathbb{P}(Z \leq - \alpha_n^{(2)} \mid \mathbf{X}) = \Phi(- \alpha_n^{(2)}) = 1 - \Phi(\alpha_n^{(2)}).\]

Using the first moment of the truncated normal [37], we have align [_ols,i- _ols,i, ] &= ^_i - + _n   [Z Z _n^(1), ]
&= ^_i - + _n and align [_ols,i+ _ols,i-, ] &= ^_i + + _n [Z Z -_n^(2)]
&= ^_i + - _n .

Therefore align [__n, i ] &= [_ols,i- _ols,i, ]   (_ols,i )
& + [_ols,i+ _ols,i-, ]   (_ols,i- )
&= (^_i - ) (1 - (_n^(1))) + _n (_n^(1)) + (^_i + ) (-_n^(2)) - _n (-_n^(2))
&= (^_i - ) (1 - (_n^(1))) + (^_i + ) (-_n^(2)) + _n ((_n^(1)) - (-_n^(2))).

Note that \(\varphi'(x) = -x \varphi(x)\). A straightforward study of the behavior of the function \(x \mapsto x \varphi(x)\) shows it is bounded. We denote the maximum of its absolute value by \(M\).

Using the mean value inequality for \(\varphi\), we have align |_n ((_n^(1)) - (-_n^(2)))| &_n |_n^(1) - (-_n^(2))| |’|
&M _n |_n^(1) - (-_n^(2))|
&= M _n ( - ^_i + + ^_i)
&= 2 M .

Therefore, since \(\lambda_n = o(n)\), \(\tilde{\tau}_n (\varphi(\alpha_n^{(1)}) - \varphi(-\alpha_n^{(2)}))\) goes to \(0\) in \(L^1\).

We first consider \(\beta^\star_i > 0\).

Since \(\frac{\mathbf{X}^\top \mathbf{X}}{n} \stackrel{\mathrm{a.s.}}{\longrightarrow}\mathbb{E}[X_0 X_0^\top]\) (strong law of large numbers), and \(\lambda_n = o(n)\), we have \(\tilde{\tau}_n \stackrel{\mathrm{a.s.}}{\longrightarrow}0^+\), and using the continuous mapping theorem, \(\alpha_n^{(1)} \stackrel{\mathrm{a.s.}}{\longrightarrow}-\infty\) and \(\alpha_n^{(2)} \stackrel{\mathrm{a.s.}}{\longrightarrow}+\infty\). \(\Phi\) is continuous bounded so we get \(L^1\) convergence of \(\Phi(\alpha_n^{(1)})\) and \(\Phi(-\alpha_n^{(2)})\) to 0. By putting everything together, we obtain \[\mathbb{E}[\hat{\beta}_{\lambda_n, i}] = \mathbb{E}[\mathbb{E}[\hat{\beta}_{\lambda_n, i} \mid \mathbf{X}]] \rightarrow \beta^\star_i.\]

When \(\beta^\star_i < 0\), we show in a similar manner that \(\mathbb{E}[\hat{\beta}_{\lambda_n, i}] \rightarrow \beta^\star_i\).

If \(\beta^\star_i = 0\), \(\alpha_n^{(1)} = \alpha_n^{(2)}\) so \(1 - \Phi(\alpha_n^{(1)}) = \Phi(-\alpha_n^{(2)})\) and \(\varphi(\alpha_n^{(1)}) = \varphi(-\alpha_n^{(2)})\) which leads to \(\mathbb{E}[\hat{\beta}_{\lambda_n, i} \mid \mathbf{X}] = 0\) and thus \(\mathbb{E}[\hat{\beta}_{\lambda_n, i}] = 0\).

Thus, we have convergence component-wise and can conclude \(\mathbb{E}[\hat{\beta}_{\lambda_n}] \rightarrow \beta^\star\).

We now show that \(\mathbb{E}[\hat{\beta}_{\lambda_n, i} \hat{\beta}_{\lambda_n, j}] \rightarrow \beta^\star_i \beta^\star_j\).

Note that align [__n, i __n, j - ^_i ^_j] = [(__n, i - ^_i) __n, j] + ^_i [__n, j - ^_j] where, using Cauchy–Schwarz and the fact that \((a + b)^2 \leq 2 (a^2 + b^2)\), align |[(__n, i - ^_i) __n, j]| .

Therefore, proving \(\mathbb{E}[\hat{\beta}_{\lambda_n, i} \hat{\beta}_{\lambda_n, j}] \rightarrow \beta^\star_i \beta^\star_j\) for all \(i,j\) comes down to proving \(\mathbb{E}[(\hat{\beta}_{\lambda_n, i} - \beta^\star_i)^2] \rightarrow 0\) for all \(i\) given that we have already shown \(\mathbb{E}[\hat{\beta}_{\lambda_n, i}] \rightarrow \beta^\star_i\) for all \(i\).

As a reminder, we have \[\hat{\beta}_{\lambda_n, i} = \text{sign}(\hat{\beta}_{\mathrm{\small ols},i}) \begin{cases} |\hat{\beta}_{\mathrm{\small ols},i}| - \frac{\lambda_n}{n} & \text{if } |\hat{\beta}_{\mathrm{\small ols},i}| \geq \frac{\lambda_n}{n}\\ 0 & \text{if } |\hat{\beta}_{\mathrm{\small ols},i}| < \frac{\lambda_n}{n}\\ \end{cases}\] thus align &[(__n, i - ^_i)^2 ]
&= [(_ols,i- ^_i - )^2 _ols,i, ]   (_ols,i )
& + [(_ols,i- ^_i + )^2 _ols,i-, ]   (_ols,i- ).

We introduce an intermediate lemma where we derive the second and fourth moments of the truncated normal. The second moment we will use right now and the fourth moment will be used later.

Proof First, we can derive a recursive formula on the \(m_k\)’s using integration by parts with the fact that \(\varphi'(x) = - x \varphi(x)\). For \(k \in \mathbb{N}\), we have align m_k+2 &= _a^b dx = []_a^b + (k+1) _a^b dx
&= + (k+1) m_k. Since \(m_0 = \int_a^b \frac{\varphi(x)}{\Phi(b) - \Phi(a)} dx = 1\), we immediately obtain \(m_2 = 1 + \frac{a \varphi(a) - b \varphi(b)}{\Phi(b) - \Phi(a)}\). And consequently, we have align m_4 &= + 3 m_2 = + 3 (1 + )
&= 3 + . 0◻

Using the second moment derived in 5, we have align &[(_ols,i- ^_i - )^2 _ols,i, ]
&= [(_n Z - )^2 Z _n^(1), ]
&= _n^2   [Z^2 Z _n^(1), ] - 2 _n [Z Z _n^(1), ] +
&= _n^2 (1 + ) - 2 _n + and align &[(_ols,i- ^_i + )^2 _ols,i-, ]
&= _n^2   [Z^2 Z -_n^(2), ] + 2 _n [Z Z -_n^(2), ] +
&= _n^2 (1 + ) - 2 _n + .

Thus align &[(__n, i - ^_i)^2 ]
&= [(_ols,i- ^_i - )^2 _ols,i, ]   (_ols,i )
& + [(_ols,i- ^_i + )^2 _ols,i-, ]   (_ols,i- )
&= _n^2 (1 - (_n^(1)) + _n^(1) (_n^(1))) - 2 _n (_n^(1)) + (1 - (_n^(1)))
& + _n^2 ((-_n^(2)) + _n^(2) (-_n^(2))) - 2 _n (-_n^(2)) + (-_n^(2)).

For \(X_i \overset{\textrm{\tiny\textrm{i.i.d.}}}{\sim}\mathcal{N}(0, \mathbf{I})\), we know \((\mathbf{X}^\top \mathbf{X})^{-1} \sim W_p^{-1}(\mathbf{I}, n)\), therefore \(\mathbb{E}[(\mathbf{X}^\top \mathbf{X})^{-1}] = \frac{\mathbf{I}}{n-p-1}\) and \(\mathbb{E}[(\frac{\mathbf{X}^\top \mathbf{X}}{n})_{i, i}^{-1}] = \frac{n}{n-p-1} = o(n)\).

Thus, using Jensen’s inequality, \(\mathbb{E}[\sqrt{(\frac{\mathbf{X}^\top \mathbf{X}}{n})_{i, i}^{-1}}] \leq \sqrt{\mathbb{E}[(\frac{\mathbf{X}^\top \mathbf{X}}{n})_{i, i}^{-1}]} = \sqrt{\frac{n}{n-p-1}} = o(\sqrt{n})\).

As a reminder, \(\tilde{\tau}_n = \frac{\tau}{\sqrt{n}} \sqrt{(\frac{\mathbf{X}^\top \mathbf{X}}{n})_{i,i}^{-1}}\). We then have \(L^1\) convergence of both \(\tilde{\tau}_n\) and \(\tilde{\tau}_n^2\) to 0. As previously mentioned, the function \(x \mapsto x \varphi(x)\) is bounded. Since \(\Phi\) and \(\varphi\) are also bounded, and \(\lambda_n = o(n)\), then \[\mathbb{E}[(\hat{\beta}_{\lambda_n, i} - \beta^\star_i)^2] = \mathbb{E}[\mathbb{E}[(\hat{\beta}_{\lambda_n, i} - \beta^\star_i)^2 \mid \mathbf{X}]] \rightarrow 0.\]

Therefore, we get \[\mathbb{E}[\hat{\beta}_{\lambda_n} \hat{\beta}_{\lambda_n}^\top] \rightarrow \beta^\star\beta^{\star\top}.\]

We can then conclude that \(\sigma^2(h_n^{\mathrm{sing}}) \rightarrow 2 \tau^4\) by 4.

15 Proof of 9:↩︎

As a reminder, to study the loss stability, we consider \(Z'_1 = (X'_1, Y'_1)\) an \(\textrm{i.i.d.}\) copy of \(Z_1 = (X_1, Y_1)\) used as replacement for the first point of the training set.

Define the vector \(V \triangleq(Y'_1 - X_1^{\prime\top} \beta^\star) X'_1 - (Y_1 - X_1^\top \beta^\star) X_1\) and the symmetric matrix \(M \triangleq- (V \beta^{\star\top}+ \beta^\star V^\top)\).

Starting from the expression of \(\gamma(h_n)\) in 3 and using the fact that \(X_0 \sim \mathcal{N}(0, \mathbf{I})\), we have align (h_n^) = [(2 (Y_0 X_0^- ^) (‘__n - __n) + \(\operatorname{tr}\)((X_0 X_0^- ) (__n __n^-’__n __n^’)))^2].

We will show that align (h_n^) ~ [(2 (Y_0 X_0^- ^) V + \(\operatorname{tr}\)((X_0 X_0^- ) M))^2]. by proving that the difference align W_n &(2 (Y_0 X_0^- ^) (‘__n - __n) + \(\operatorname{tr}\)((X_0 X_0^- ) (__n __n^-’__n __n^’)))^2
& - (2 (Y_0 X_0^- ^) + \(\operatorname{tr}\)((X_0 X_0^- ) ))^2 goes to \(0\) in \(L^1\).

Since \(a^2 - b^2 = (a - b)(a + b)\), we have align W_n &= (D_n,1 + D_n,2) (S_n,1 + S_n,2). where align D_n,1 & (Y_0 X_0^- ^) (‘__n - __n - ),
D_n,2 &\(\operatorname{tr}\)((X_0 X_0^- ) (__n __n^-’__n __n^’ - )),
S_n,1 & (Y_0 X_0^- ^) (‘__n - __n + ),
S_n,2 &\(\operatorname{tr}\)((X_0 X_0^- ) (__n __n^-’__n __n^’ + )).

Using Cauchy–Schwarz and the fact that \((a+b)^2 \leq 2 (a^2 + b^2)\), align [|W_n|] &
& .

To obtain convergence of \(W_n\) to 0 in \(L^1\), we will thus prove that \(\mathbb{E}[D_{n,1}^2] \rightarrow 0\), \(\mathbb{E}[S_{n,1}^2] = O(1)\), \(\mathbb{E}[D_{n,2}^2] \rightarrow 0\) and \(\mathbb{E}[S_{n,2}^2] = O(1)\).

We have align [D_n,1^2] &= [4 (Y_0 X_0^- ^) (‘__n - __n - ) (’__n - __n - )^(Y_0 X_0 - ^)]
&= [4 \(\operatorname{tr}\)((Y_0 X_0^- ^) (‘__n - __n - ) (’__n - __n - )^(Y_0 X_0 - ^))]
&= [4 \(\operatorname{tr}\)((Y_0 X_0 - ^) (Y_0 X_0^- ^) (‘__n - __n - ) (’__n - __n - )^)]
&= 4 \(\operatorname{tr}\)([(Y_0 X_0 - ^) (Y_0 X_0^- ^) (‘__n - __n - ) (’__n - __n - )^])
&= 4 \(\operatorname{tr}\)([(Y_0 X_0 - ^) (Y_0 X_0^- ^)] [(‘__n - __n - ) (’__n - __n - )^]) as \(\hat{\beta}'_{\lambda_n} - \hat{\beta}_{\lambda_n} - \frac{V}{n}\) is a function of the training points and using independence of \(Z_0\) from the training points.

By Cauchy–Schwarz, for all \(i, j\), align [|(‘__n, i - __n, i - ) (’__n, j - __n, j - )|] & thus, if we show \(\mathbb{E}[(\hat{\beta}'_{\lambda_n, i} - \hat{\beta}_{\lambda_n, i} - \frac{V_i}{n})^2] \rightarrow 0\) for all \(i\), then we obtain align [(‘__n - __n - ) (’__n - __n - )^] and therefore \(\mathbb{E}[D_{n,1}^2] \rightarrow 0\). We are going to hold off on proving \(\mathbb{E}[(\hat{\beta}'_{\lambda_n, i} - \hat{\beta}_{\lambda_n, i} - \frac{V_i}{n})^2] \rightarrow 0\) as we will actually show the stronger convergence \(\mathbb{E}[(\hat{\beta}'_{\lambda_n, i} - \hat{\beta}_{\lambda_n, i} - \frac{V_i}{n})^4] \rightarrow 0\) in the context of proving \(\mathbb{E}[D_{n,2}^2] \rightarrow 0\).

With similar computations and upper-bounding, we can show that \(\mathbb{E}[S_{n,1}^2] = O(1)\) if we prove that for all \(i\), \(\mathbb{E}[(\hat{\beta}'_{\lambda_n, i} - \hat{\beta}_{\lambda_n, i} + \frac{V_i}{n})^2] = O(1)\).

As we have shown in 14 that the soft-thresholding Lasso estimator is consistent for \(\beta^\star\) in \(L^2\) when \(\lambda_n = o(n)\), both \(\mathbb{E}[\hat{\beta}_{\lambda_n, i}^2]\) and \(\mathbb{E}[\hat{\beta}_{\lambda_n, i}^{'2}]\) are bounded and thus \(\mathbb{E}[(\hat{\beta}'_{\lambda_n, i} - \hat{\beta}_{\lambda_n, i} + \frac{V_i}{n})^2] = O(1)\) since \((\hat{\beta}'_{\lambda_n, i} - \hat{\beta}_{\lambda_n, i} + \frac{V_i}{n})^2 \leq 3 (\hat{\beta}_{\lambda_n, i}^{'2} + \hat{\beta}_{\lambda_n, i}^2 + \frac{V_i^2}{n^2})\) by Cauchy–Schwarz.

We now focus on proving \(\mathbb{E}[D_{n,2}^2] \rightarrow 0\).

We have align D_n, 2 &= \(\operatorname{tr}\)((X_0 X_0^- ) (__n __n^- ‘__n __n^’ - ))
&= X_0^(__n __n^- ‘__n __n^’ - ) X_0 - \(\operatorname{tr}\)(__n __n^- ‘__n __n^’ - )
&= _i, j (X_0,i X_0,j - ) (__n, i __n, j - ‘__n, i’__n, j - )
&= _i, j U_i, j (__n, i __n, j - ‘__n, i’__n, j - ) where \(U_{i, j} \triangleq X_{0,i} X_{0,j} - \mathbb{1}\mathopen{}\mathclose{\left[{i = j}}\right]\), and thus align D_n, 2^2 &= _i, j, k, l U_i, j U_k, l (__n, i __n, j - ‘__n, i’__n, j - ) (__n, k __n, l - ‘__n, k’__n, l - ).

Using independence of \(Z_0\) and the training points, we have align [D_n, 2^2] &= _i, j, k, l [U_i, j U_k, l] [(__n, i __n, j - ‘__n, i’__n, j - ) (__n, k __n, l - ‘__n, k’__n, l - )] where, using Cauchy–Schwarz, align &[|(__n, i __n, j - ‘__n, i’__n, j - ) (__n, k __n, l - ‘__n, k’__n, l - )|]
&.

We thus want to show \(\mathbb{E}[(\hat{\beta}_{\lambda_n, i} \hat{\beta}_{\lambda_n, j} - \hat{\beta}'_{\lambda_n, i} \hat{\beta}'_{\lambda_n, j} - \frac{M_{i,j}}{n})^2] \rightarrow 0\) for all \(i, j\).

Since \(M = - (V \beta^{\star\top}+ \beta^\star V^\top)\), we have \(M_{i,j} = - V_i \beta^\star_j - \beta^\star_i V_j\) and then align &__n, i __n, j - ‘__n, i’__n, j -
&= __n, i __n, j - ‘__n, i’__n, j + ^_j + ^_i
&= -(‘__n, i - __n, i - ) __n, j -’__n, i (‘__n, j - __n, j - ) - (__n, j - ^_j) - (’__n, i - ^_i) .

By Cauchy–Schwarz, align &(__n, i __n, j - ‘__n, i’__n, j - )^2
&= ((‘__n, i - __n, i - ) __n, j +’__n, i (‘__n, j - __n, j - ) + (__n, j - ^_j) + (’__n, i - ^_i) )^2
& ((‘__n, i - __n, i - )^2 __n, j^2 + __n, i^’2 (’__n, j - __n, j - )^2 + (__n, j - ^_j)^2 + (‘__n, i - ^_i)^2 ) and the probability version of Cauchy–Schwarz yields align &[(__n, i __n, j -’__n, i ’__n, j - )^2]
& ( +
& + + ). Hence, we will get \(\mathbb{E}[D_{n,2}^2] \rightarrow 0\) if we prove that for all \(i\)

• \(\mathbb{E}[(\hat{\beta}_{\lambda_n, i} - \beta^\star_i)^4] \rightarrow 0\), the proof will be the same for \(\mathbb{E}[(\hat{\beta}'_{\lambda_n, i} - \beta^\star_i)^4] \rightarrow 0\),

• \(\mathbb{E}[(\hat{\beta}'_{\lambda_n, i} - \hat{\beta}_{\lambda_n, i} - \frac{V_i}{n})^4] \rightarrow 0\).

Note that we will automatically get \(L^2\) convergence of \(\hat{\beta}'_{\lambda_n, i} - \hat{\beta}_{\lambda_n, i} - \frac{V_i}{n}\) to 0 for all \(i\), which implies \(\mathbb{E}[D_{n,1}^2] \rightarrow 0\) as mentioned earlier.

We now introduce a lemma that will allow us to upper-bound quantities of interest.

Proof For \((x_1, \ldots, x_k), (y_1, \ldots, y_k) \in \mathbb{R}^k\) and \(p, q \in (1, +\infty)\) such that \(\frac{1}{p} + \frac{1}{q} = 1\), Hölder’s inequality gives us align _i=1^k |x_i y_i| (_i=1^k |x_i|^p)^ (_i=1^k |y_i|^q)^ and therefore the lemma is an application of it with \(x_i = a_i, y_i = 1, p = \ell\). 0◻

Combining 6 for \(\ell = 4\) with similar computations and upper-bounding as above, we can show that \(\mathbb{E}[S_{n,2}^2]\) is bounded if for all \(i\), \(\mathbb{E}[\hat{\beta}_{\lambda_n, i}^4]\) and \(\mathbb{E}[\hat{\beta}_{\lambda_n, i}^{'4}]\) are bounded, which automatically comes from the \(L^4\) convergence of the soft-thresholding Lasso estimator to \(\beta^\star\) needed for \(\mathbb{E}[D_{n,2}^2] \rightarrow 0\).

We start by showing \(\mathbb{E}[(\hat{\beta}_{\lambda_n, i} - \beta^\star_j)^4] \rightarrow 0\).

As a reminder, we have \[\hat{\beta}_{\lambda_n, i} = \text{sign}(\hat{\beta}_{\mathrm{\small ols},i}) \begin{cases} |\hat{\beta}_{\mathrm{\small ols},i}| - \frac{\lambda_n}{n} & \text{if } |\hat{\beta}_{\mathrm{\small ols},i}| \geq \frac{\lambda_n}{n}\\ 0 & \text{if } |\hat{\beta}_{\mathrm{\small ols},i}| < \frac{\lambda_n}{n}\\ \end{cases}\] thus, using \((a + b)^4 \leq 8 (a^4 + b^4)\), which is an application of 6 for \(\ell = 4\), align &[(__n, i - ^_i)^4 ]
&= [(_ols,i- ^_i - )^4 _ols,i, ] (_ols,i )
& + [(_ols,i- ^_i + )^4 _ols,i-, ] (_ols,i- )
& ([(_ols,i- ^_i)^4 _ols,i, ] + ) (_ols,i )
& + 8 ([(_ols,i- ^_i)^4 _ols,i-, ] + ) (_ols,i- ).

Since \(\hat{\beta}_{\mathrm{\small ols}}\mid \mathbf{X} \sim \mathcal{N}(\beta^\star, \tau^2 (\mathbf{X}^\top \mathbf{X})^{-1})\), we can write \(\hat{\beta}_{\mathrm{\small ols},i}= \beta^\star_i + \tilde{\tau}_n Z\) where \(\tilde{\tau}_n = \frac{\tau}{\sqrt{n}} \sqrt{(\frac{\mathbf{X}^\top \mathbf{X}}{n})_{i,i}^{-1}}\) and \(Z \mid \mathbf{X} \sim \mathcal{N}(0, 1)\). Note that we could have \(i\) as a subscript of \(\tilde{\tau}_n\) and \(Z\), but we will only consider one \(i\) at a time in our computations and we can thus omit this subscript for both of them for the sake of notational simplicity, and we will also omit it for some additional notation we define in the rest of the proof.

Define \(\alpha_n^{(1)} = \frac{1}{\tilde{\tau}_n} (\frac{\lambda_n}{n} - \beta^\star_i)\) and \(\alpha_n^{(2)} = \frac{1}{\tilde{\tau}_n} (\frac{\lambda_n}{n} + \beta^\star_i)\).

Using the fourth moment derived in 5, we have align & [(_ols,i- ^_i)^4 _ols,i, ]
&= [(_n Z)^4 Z _n^(1), ]
&= _n^4   [Z^4 Z _n^(1), ]
&= _n^4 (3 + ) and align &[(_ols,i- ^_i)^4 _ols,i-, ]
&= _n^4   [Z^4 Z -_n^(2), ]
&= _n^4 (3 + ).

Since \(\mathbb{P}(\hat{\beta}_{\mathrm{\small ols},i}\geq \frac{\lambda_n}{n} \mid \mathbf{X}) = 1 - \Phi(\alpha_n^{(1)})\) and \(\mathbb{P}(\hat{\beta}_{\mathrm{\small ols},i}\leq -\frac{\lambda_n}{n} \mid \mathbf{X}) = \Phi(-\alpha_n^{(2)})\), align &[(__n, i - ^_i)^4 ]
& ([(_ols,i- ^_i)^4 _ols,i, ] + ) (_ols,i )
& + 8 ([(_ols,i- ^_i)^4 _ols,i-, ] + ) (_ols,i- )
&= 8 (3 _n^4 (1 - (_n^(1))) + _n^4 ((_n^(1))^3 + 3 _n^(1)) (_n^(1)) + (1 - (_n^(1))))
& + 8 (3 _n^4 (-_n^(2)) + _n^4 ((_n^(2))^3 + _n^(2)) (-_n^(2)) + (-_n^(2))).

For \(X_i \overset{\textrm{\tiny\textrm{i.i.d.}}}{\sim}\mathcal{N}(0, \mathbf{I})\), we know \((\mathbf{X}^\top \mathbf{X})^{-1} \sim W_p^{-1}(\mathbf{I}, n)\) and then the diagonal element \((\mathbf{X}^\top \mathbf{X})_{i, i}^{-1}\) follows an inverse gamma distribution with shape parameter \(\frac{n-p+1}{2}\) and scale parameter \(\frac{1}{2}\). Therefore, \(\mathbb{E}[((\mathbf{X}^\top \mathbf{X})_{i, i}^{-1})^2] = \frac{1}{(n-p-1)(n-p-3)}\) and \(\mathbb{E}[((\frac{\mathbf{X}^\top \mathbf{X}}{n})_{i, i}^{-1})^2] = \frac{n^2}{(n-p-1)(n-p-3)} = o(n^2)\).

As a reminder, \(\tilde{\tau}_n = \frac{\tau}{\sqrt{n}} \sqrt{(\frac{\mathbf{X}^\top \mathbf{X}}{n})_{i,i}^{-1}}\). We then have \(L^1\) convergence of \(\tilde{\tau}_n^4\) to 0. As previously mentioned, the function \(x \mapsto x \varphi(x)\) is bounded. Similarly, a straightforward study of the behavior of the function \(x \mapsto x^3 \varphi(x)\) shows it is bounded. Since \(\Phi\) is also bounded, and \(\lambda_n = o(n)\), then \[\mathbb{E}[(\hat{\beta}_{\lambda_n, i} - \beta^\star_i)^4] = \mathbb{E}[\mathbb{E}[(\hat{\beta}_{\lambda_n, i} - \beta^\star_i)^4 \mid \mathbf{X}]] \rightarrow 0.\]

We now show that \(\mathbb{E}[(\hat{\beta}'_{\lambda_n, i} - \hat{\beta}_{\lambda_n, i} - \frac{V_i}{n})^4] \rightarrow 0\).

We have align &‘__n, i - __n, i
&= (_ols,i’) (|_ols,i’| - )_+ - (_ols,i) (|_ols,i| - )_+
&= (_ols,i’)

-(_ols,i)

&=

-

.

As an intermediate step, we need to show \(\hat{\beta}'_{\text{OLS}} - \hat{\beta}_{\mathrm{\small ols}}- \frac{V}{n} \stackrel{\mathrm{a.s.}}{\longrightarrow}0\).

Let \(\mathbf{\tilde{X}} \triangleq(X_2, \ldots, X_n)^\top\) be the matrix of regressors for the training points except for the first one that is being changed.

We have align &‘_ - _ols
&= ( )^-1 ^ - (^)^-1 ^
&= (^ + X’_1 X_1^)^-1 (^ + Y’_1 X’_1) - (^ + X_1 X_1^)^-1 (^ + Y_1 X_1)
&= [(^ + X’_1 X_1^)^-1 - (^ + X_1 X_1^)^-1] ^
& + (^ + X’_1 X_1^)^-1 Y’_1 X’_1 - (^ + X_1 X_1^)^-1 Y_1 X_1.

Using the Sherman–Morrison–Woodbury formula, align (^ + X_1 X_1^)^-1 &= (^)^-1 - (^)^-1 X_1 ( + X_1^(^)^-1 X_1)^-1 X_1^(^)^-1
&= ()^-1 - ()^-1 X_1 ( + X_1^()^-1 X_1)^-1 X_1^()^-1
&= A_n - B_n where, by the strong law of large numbers,

• \(A_n \triangleq(\frac{\mathbf{\tilde{X}}^\top \mathbf{\tilde{X}}}{n})^{-1} \stackrel{\mathrm{a.s.}}{\longrightarrow}\mathbb{E}[X_0 X_0^\top]^{-1} = \mathbf{I}\),

• \(B_n \triangleq(\frac{\mathbf{\tilde{X}}^\top \mathbf{\tilde{X}}}{n})^{-1} X_1 (\mathbf{I} + \frac{1}{n} X_1^\top (\frac{\mathbf{\tilde{X}}^\top \mathbf{\tilde{X}}}{n})^{-1} X_1)^{-1} X_1^\top (\frac{\mathbf{\tilde{X}}^\top \mathbf{\tilde{X}}}{n})^{-1} \stackrel{\mathrm{a.s.}}{\longrightarrow}X_1 X_1^\top\).

Similarly, align (^ + X’_1 X_1^)^-1 &= A_n - B’_n with align B’_n ()^-1 X’_1 ( + X_1^ ()^-1 X’_1)^-1 X_1^ ()^-1 X’_1 X_1^.

Then align ‘_ - _ols &= (B_n - B’_n) ^ + ( A_n - B’_n) Y’_1 X’_1 - ( A_n - B_n) Y_1 X_1
&= (B_n - B’_n) + ( A_n - B’_n) Y’_1 X’_1 - ( A_n - B_n) Y_1 X_1 where \(\frac{\mathbf{\tilde{X}}^\top \mathbf{\tilde{Y}}}{n} \stackrel{\mathrm{a.s.}}{\longrightarrow}\mathbb{E}[Y_0 X_0] = \beta^\star\), by the strong law of large numbers.

Therefore, align n (‘_ - _ols) &(X_1 X_1^- X’_1 X_1^) ^+ Y’_1 X’_1 - Y_1 X_1
& = (Y’_1 - X_1^ ^) X’_1 - (Y_1 - X_1^^) X_1
& = V.

We can write align (‘__n, i - __n, i - )^4 &= (_ols,i’ - _ols,i- )^4  
& + (_ols,i’ - _ols,i- )^4  
& + (_ols,i’ - _ols,i- 2 - )^4  
& + (_ols,i’ - _ols,i+ 2 - )^4  
& + (_ols,i’ - - )^4  
& + (_ols,i’ + - )^4  
& + (_ols,i- + )^4  
& + (_ols,i+ + )^4  
& + ()^4   and we have a similar expression for \((\hat{\beta}'_{\lambda_n, i} - \hat{\beta}_{\lambda_n, i} - \frac{V_i}{n})^6\) with terms taken to the sixth power.

Since \(\hat{\beta}_{\mathrm{\small ols}}\mid \mathbf{X} \sim \mathcal{N}(\beta^\star, \tau^2 (\mathbf{X}^\top \mathbf{X})^{-1})\) and we can bound the central moments of a Normal with the powers of its variance, there exists \(C > 0\) such that \(\mathbb{E}[(\hat{\beta}_{\mathrm{\small ols},i}- \beta^\star_i)^6 \mid \mathbf{X}] \leq C (\tau^2 (\mathbf{X}^\top \mathbf{X})_{i, i}^{-1})^3 = C \tau^6 ((\mathbf{X}^\top \mathbf{X})_{i, i}^{-1})^3\).

For \(X_i \overset{\textrm{\tiny\textrm{i.i.d.}}}{\sim}\mathcal{N}(0, \mathbf{I})\), we know \((\mathbf{X}^\top \mathbf{X})^{-1} \sim W_p^{-1}(\mathbf{I}, n)\) and then the diagonal element \((\mathbf{X}^\top \mathbf{X})_{i, i}^{-1}\) follows an inverse gamma distribution with shape parameter \(\frac{n-p+1}{2}\) and scale parameter \(\frac{1}{2}\). Therefore, \(\mathbb{E}[((\mathbf{X}^\top \mathbf{X})_{i, i}^{-1})^3] = \frac{1}{(n-p-1)(n-p-3)(n-p-5)}\), which means \(\mathbb{E}[(\hat{\beta}_{\mathrm{\small ols},i}- \beta^\star_i)^6]\) and thus \(\mathbb{E}[\hat{\beta}_{\mathrm{\small ols},i}^6]\), by an application of 6 for \(\ell = 6\), are bounded. Similarly, \(\mathbb{E}[\hat{\beta}_{\mathrm{\small ols},i}^{'6}]\) is bounded.

Consequently, since \(\lambda_n = o(n)\) and \(\mathbb{E}[\hat{\beta}_{\mathrm{\small ols},i}^6]\) and \(\mathbb{E}[\hat{\beta}_{\mathrm{\small ols},i}^{'6}]\) are bounded, the almost sure convergence of the fourth moment turns into \(L^1\) convergence to \(0\).

Therefore, align (h_n^) &~ [(2 (Y_0 X_0^- ^) V + \(\operatorname{tr}\)((X_0 X_0^- ) M))^2]
&= [(2 Y_0 X_0^V - 2 ^V - \(\operatorname{tr}\)((X_0 X_0^- ) (V ^+ ^V^)))^2]
&= [(2 Y_0 X_0^V - 2 X_0^^X_0^V)^2]
&= [(2 (Y_0 - X_0^^) X_0^V)^2] where \(V = (Y'_1 - X_1^{\prime\top} \beta^\star) X'_1 - (Y_1 - X_1^\top \beta^\star) X_1\).

16 Proof of 3:↩︎

Instantiate the ST notation of 1, and define the shorthand align V &(Y_0 X_0^- ^[X_0 X_0^])^, M (X_0 X_0^- [X_0 X_0^]),
A & V^[^lasso__n+_n - ^lasso__n] + \(\operatorname{tr}\)(M, [^lasso__n^lasso__n- ^lasso__n+_n^lasso__n+_n]),
B & V^[__n+_n - __n] + \(\operatorname{tr}\)(M, [__n__n^- __n+_n__n+_n^]). We will establish the \(\sigma^2(\tilde{h}_n^{\mathrm{diff}})\) upper bound in 16.1 and the \(\gamma(\tilde{h}_n^{\mathrm{diff}})\) lower bound in 16.2.

16.1 \(\sigma^2(\tilde{h}_n^{\mathrm{diff}})\) upper bound↩︎

By 2 and Cauchy-Schwarz, we have align &^2(_n^)- ^2(h_n^) = [A^2 - B^2] = [(A-B)^2]+[(A-B)2B]
& [(A-B)^2]+2 = [(A-B)^2] + 2. Meanwhile, Cauchy-Schwarz, the triangle inequality, and the definition of the operator norm imply align |A-B| &= |V^[^lasso__n+_n-__n+_n + __n- ^lasso__n]
&+ \(\operatorname{tr}\)(M, [^lasso__n^lasso__n- __n__n^+^lasso__n+_n^lasso__n+_n-__n+_n__n+_n^]])|
& V̍̍_2[^̍lasso__n+_n-__n+_n̍_2 + _̍_n- ^lasso__n̍_2]
&+ M̍̍_[ ^̍lasso__n-__n̍_2(^̍lasso__n-__n̍_2+ 2_̍_n̍_2)
&+ ^̍lasso__n+_n-__n+_n̍_2(^̍lasso__n+_n-__n+_n̍_2+2_̍_n+_n̍_2)] Now, since \(\mathbb{E}[ \|{\hat{\beta}^{\mathrm{\small lasso}}_{\lambda_n}-\hat{\beta}_{\lambda_n}}\|_2^2] = O(\frac{\lambda_n^2}{n^3})\) by 1, \(\mathbb{E}[\|{\hat{\beta}_{\lambda_n}}\|_2^2]\leq \mathbb{E}[\|{\hat{\beta}_{\mathrm{\small ols}}}\|_2^2] = O(1)\) by [38], and \(\delta_n = O(\lambda_n)\), we have \(\mathbb{E}[(A-B)^2] = O(\frac{\lambda_n^2+(\lambda_n+\delta_n)^2}{n^{3}}) = O(\frac{\lambda_n^2}{n^{3}})\). Therefore, since \(\sigma^2(h_n^{\mathrm{diff}})= O(\frac{1}{n^2})\) by 1 and \(\lambda_n = O(\sqrt{n})\), we can conclude from inequality [lsdiff-sdiff] that \(\sigma^2(\tilde{h}_n^{\mathrm{diff}})= O(\frac{1}{n^2})\) as well.

16.2 \(\gamma(\tilde{h}_n^{\mathrm{diff}})\) lower bound↩︎

Let \((A',B')\) be an exchangeable copy of \((A,B)\) in which the first datapoint \(Z_1\) has been replaced by an i.i.d. copy \(Z'_1\). Then, by the triangle inequality and exchangeability, align - &= -
& - - = -2. Since \(\sqrt{\gamma(h_n^{\mathrm{diff}})} =\Omega(\frac{1}{n^{5/4}})\) by 1, \(\lambda_n = O(\sqrt{n})\) by assumption, and \(\sqrt{\mathbb{E}[(A-B)^2]} = O(\frac{\lambda_n}{n^{3/2}}) = O(\frac{1}{n})\) by 16.1, we also have \(\sqrt{\gamma(\tilde{h}_n^{\mathrm{diff}})} = \Omega(\frac{1}{n^{5/4}})\) by [gdiff-lgdiff].

17 Proof of 4:↩︎

Instantiate the ST notation of 2, and define the shorthand align V &(Y_0 X_0^- ^[X_0 X_0^])^, M (X_0 X_0^- [X_0 X_0^]),
B & -V^[^lasso__n] + \(\operatorname{tr}\)(M, [^lasso__n^lasso__n]),
A & -V^[__n] + \(\operatorname{tr}\)(M, [__n__n^]). We will establish the \(\sigma^2(\tilde{h}_n^{\mathrm{sing}})\) lower bound in 17.1 and the \(\gamma(\tilde{h}_n^{\mathrm{sing}})\) upper bound in 17.2.

17.1 \(\sigma^2(\tilde{h}_n^{\mathrm{sing}})\) lower bound↩︎

By 3 and Cauchy-Schwarz, we have align &^2(h_n^)- ^2(_n^) = [A^2 - B^2] = [(A-B)2A] - [(A-B)^2]
& = 2. Meanwhile, Cauchy-Schwarz, the triangle inequality, and the definition of the operator norm imply align |A-B| &= |V^[__n- ^lasso__n] + \(\operatorname{tr}\)(M, [^lasso__n^lasso__n- __n__n^])|
& V̍̍_2 [_̍_n- ^lasso__n̍_2] + M̍̍_[ ^̍lasso__n-__n̍_2(^̍lasso__n-__n̍_2+ 2_̍_n̍_2). Now, since \(\mathbb{E}[ \|{\hat{\beta}^{\mathrm{\small lasso}}_{\lambda_n}-\hat{\beta}_{\lambda_n}}\|_2^2] = O(\frac{\lambda_n^2}{n^3})\) by 1 and \(\mathbb{E}[\|{\hat{\beta}_{\lambda_n}}\|_2^2]\leq \mathbb{E}[\|{\hat{\beta}_{\mathrm{\small ols}}}\|_2^2] = O(1)\) by [38], we have \(\mathbb{E}[(A-B)^2] = O(\frac{\lambda_n^2}{n^{3}})\). Therefore, since \(\sigma^2(h_n^{\mathrm{sing}})= \Theta(1)\) by 2 and \(\lambda_n = o(n^{3/2})\), we can conclude from inequality [ssing-lssing] that \(\sigma^2(\tilde{h}_n^{\mathrm{sing}})= \Omega(1)\) as well.

17.2 \(\gamma(\tilde{h}_n^{\mathrm{sing}})\) upper bound↩︎

Let \((A',B')\) be an exchangeable copy of \((A,B)\) in which the first datapoint \(Z_1\) has been replaced by an i.i.d. copy \(Z'_1\). Then, by 3, Cauchy-Schwarz, and exchangeability, align &(_n^)- (h_n^) = [(B-B’)^2 - (A-A’)^2]
&= [(B-A+A’-B’)^2] + [(B-A+A’-B’)2(A-A’)]
&[(A-B)^2] + 4 = 4[(A-B)^2] + 4. Since \(\gamma(h_n^{\mathrm{sing}})= O(\frac{1}{n^{2}})\) by 1, \(\sqrt{\mathbb{E}[(A-B)^2]} = O(\frac{\lambda_n}{n^{3/2}})\) by 17.1, and \(\lambda_n = o(n)\), we conclude from inequality [lgsing-gsing] that \(\gamma(\tilde{h}_n^{\mathrm{sing}})= O(\frac{\lambda_n^2}{n^3}+\frac{\lambda_n}{n^{7/2}}) = o(\frac{1}{n})\).

18 Experimental Setup Details↩︎

We provide additional details about the numerical experiments presented in 5.

In our simulations, we work with the following sample sizes for the full set size \(\frac{n k}{k-1}\): 100, 1,000, 10,000, 100,000, which means \(n\) takes the following values: 90, 900, 9,000, 90,000.

For simulations with the Lasso estimator, we used the implementation from \(\texttt{scikit-learn}\). For the KDE plots, we called \(\texttt{kdeplot}\) from the \(\texttt{seaborn}\) library.

We perform 50,000 replications to sample from \(\frac{\sqrt{\frac{n k}{k-1}}}{\sigma(h_n)} (\hat{R}_n- R_n)\) and \(\frac{\sqrt{\frac{n k}{k-1}}}{\hat{\sigma}_n(h_n)} (\hat{R}_n- R_n)\). We ensured reproducibility by setting random seeds at the start of all replications.

Regarding the inner cross-validation used to determine \(\lambda_n\) in each iteration of the outer cross-validation, we performed an adaptive grid search via \((k-1)\)-fold cross-validation on the training set of size \(n\), based on the initial split of the cross-validation on the full set of size \(\frac{n k}{k-1}\). For the adaptive grid search scheme, we started with powers of 10, identified the best choice of penalization, subdivided around this choice with 10 values with an exponential scaling, and did so 3 additional times to identify the optimal penalization with precision.

We now introduce two lemmas that allow us to properly estimate \(\sigma^2(h_n)\), \(\gamma(h_n)\) and \(R_n\).

Proof align ^2(h_n) &= ([h_n(Z_0, ) Z_0])
&= [[h_n(Z_0, ) Z_0]^2] - [h_n(Z_0, )]^2
&= [[h_n(Z_0, ) h_n(Z_0, ) Z_0]] - [h_n(Z_0, ) h_n(Z_0, )]
&= [h_n(Z_0, ) h_n(Z_0, )] - [h_n(Z_0, ) h_n(Z_0, )]
&= [h_n(Z_0, ) (h_n(Z_0, ) - h_n(Z_0, ))] 0◻

Proof Starting from the expression of \(\mathbb{E}[h_n(Z_0, \mathbf{Z}) \mid \mathbf{Z}]\) in 3, we have align [h_n^(Z_0, ) ] &= [Y_0^2] - 2 ^[X_0 X_0^] + \(\operatorname{tr}\)([X_0 X_0^] ^)
&= (Y_0) + [Y_0]^2 - 2 ^[X_0 X_0^] + \(\operatorname{tr}\)([X_0 X_0^] ^)
&= ^(X_0) ^+ ^2 + ([X_0]^^)^2 - 2 ^[X_0 X_0^] + \(\operatorname{tr}\)([X_0 X_0^] ^)
&= ^(X_0) ^+ ^2 + ^[X_0][X_0]^^- 2 ^[X_0 X_0^] + ^[X_0 X_0^]
&= ^2 + ^[X_0 X_0^] ^+ ^[X_0 X_0^] - 2 ^[X_0 X_0^]
&= ^2 + (^- )^[X_0 X_0^] (^- )
&= ^2 + ̍^- ̍_2^2 since the features are drawn from a distribution with mean 0 and identity covariance matrix. The other three expressions follow from the definition of the quantities. 0◻

The Monte Carlo estimation of \(\sigma^2(h_n)\) and \(\gamma(h_n)\) is based on 5,000,000 replications when using deterministic \(\lambda_n\), but on 1,000,000 when \(\lambda_n\) is selected via inner cross-validation due to computational complexity. Based on the Monte Carlo standard errors obtained for \(\sigma^2(h_n)\) and \(\gamma(h_n)\), we applied the Delta method as follows to obtain a standard error for \(r(h_n) = \frac{n \, \cdot \, \gamma(h_n)}{\sigma^2(h_n)}\). We define \(f(x, y) = \frac{n x}{y}\) and we denote by \(M\) the number of Monte Carlo replications used to estimate \(\sigma^2(h_n)\) and \(\gamma(h_n)\). Starting from the Monte Carlo standard errors \(\frac{\sigma_x}{\sqrt{M}}\) of \(\sigma^2(h_n)\) and \(\frac{\sigma_y}{\sqrt{M}}\) of \(\gamma(h_n)\), and using \(\nabla f = (\frac{n}{y}, -\frac{n x}{y^2})\), we get to a standard error for \(r(h_n)\) by computing align f(x, y)^(_x^2, _y^2) f(x, y) = + .

Denoting the Monte Carlo estimates of \(\sigma^2(h_n)\) and \(\gamma(h_n)\) by \(\hat{x}\) and \(\hat{y}\), respectively, the standard error we use for \(r(h_n)\) is then \[\frac{1}{\sqrt{M}}\sqrt{\frac{n^2 \sigma_x^2}{\hat{y}^2} + \frac{n^2 \hat{x}^2 \sigma_y^2}{\hat{y}^4}}.\]

19 Additional Experimental Results↩︎

References↩︎