Optimal Nuisance Function Tuning for Estimating a Doubly Robust Functional under Proportional Asymptotics

We greatly appreciate the reviewer’s feedback and suggestions, and would like to address them as follows:

1. ECC, not general functional:

We chose ECC as the main parameter, as in the literature it is often considered a starting point for analyzing more general doubly robust functionals (see refs. [7,10,14,15,19,25] of our paper). ECC is studied extensively in both fixed-dimensional nonparametric and ultra-high-dimensional setup under sparsity constraints (e.g., ref [19]). Moreover, as noted in those references, deriving necessary and sufficient ways of tuning nuisance functions for optimal doubly robust functional estimation has been carried out in the nonparametric fixed $p$ regime -- often using the ECC functional as a template. In contrast, we work in a more challenging proportional asymptotic regime ($p/n \to c \in (0, \infty)$), where dimension can exceed sample size and no sparsity is assumed. Here, there are no results for the necessary and sufficient ways of tuning nuisance functions. We view ECC as a foundational step for extending this theory to more general functionals, as has been achieved in other asymptotic regimes. That said, we appreciate your feedback and will revise the paper title to avoid confusion.

2. Assumption (2.1) and asymptotic normality:

We would first like to point out, for the majority of our results (Theorem 3.1 - 3.3, which shows that the estimators are $\sqrt{n}$-consistent), we do not need asymptotic normality of the covariates or errors; all we need is: (i) $X_{ij}$ are i.i.d. with mean $0$ and variance $1$ and finite $(4 + \delta)$ moment, and (ii) $(\epsilon, \mu) \overset{i.i.d.}{\sim} F$ where $F$ with mean $0$, variance $1$, correlation $\rho$, and finite $(2 + \delta)$ moment for some $\delta > 0$.

We believe the results can be extended to settings where $X_i$ is a subgaussian random variable in $\mathbb{R}^p$ with bounded sub-gaussian norm using longer Martingale arguments, but we leave it as future work, as it doesn't change the main message of the paper and allows us to keep our proofs relatively compact.

For exact asymptotic variance, we use the normality of the covariates to shorten the length of the proof significantly. One may once again use Martingale arguments along with deterministic equivalents for suitable functionals of reolvents of large covariance matrices, but it would yield a much longer proof without changing the message.

As for the linearity assumption in regression, one could extend the analysis to feature space linear models or consider kernel ridge regression for nuisance functions, leveraging tools from [1]. While the overall framework would remain similar, the extension would require handling more involved deterministic equivalents of the associated Gram matrices, leading to significantly more complex computations. Given the current length and scope of our paper, we leave this promising direction for future work.

Lastly, while we do not provide asymptotic normality, we note that the main crux of the proof shares the same complexity as that of obtaining precise asymptotic variance as we describe above.

3. Practical implementation of (3.11) and a data-driven optimal estimator:

The analytic expression of the variances indicates that given any $(\lambda_1, \lambda_2)$, the limiting variance is a deterministic function of three unknowns: $(||\alpha_0||_2^2, ||\beta_0||_2^2, \alpha_0^\top \beta_0)$. The final goal is to choose $(\lambda_1, \lambda_2)$ which minimizes these variance expressions. This also echoes the main point of our paper: the prediction-optimal tuning parameters, i.e., the value of $\lambda_1$ (resp. $\lambda_2$) that yields the minimum prediction risk of $\alpha_0$ (resp. $\beta_0$), may not yield an estimator of ECC with minimum asymptotic variance. This contrasts with the conventional wisdom in double machine learning (DML), where nuisance parameters are typically estimated by minimizing prediction error, a point which is further highlighted in our simulation section. We next present two practical strategies for consistently estimating the limiting variance and selecting optimal tuning parameters:

(i) consistent estimation of unknown quantities One may consistently estimate $(||\alpha_0||\_2^2, ||\beta_0||\_2^2, \alpha\_0^\top \beta\_0)$ and plug them into the variance expressions. To estimate $||\alpha_0||\_2^2$, one can first compute a simple plugged-in estimator $||\hat \alpha(\lambda)||\_2^2$, ($\hat \alpha(\lambda)$ is a ridge estimator of $\alpha\_0$) with tuning parameter $\lambda$, and then apply bias correction to obtain a consistent estimator. Similar procedures apply to $||\beta_0||\_2^2$ and $\alpha\_0^\top \beta\_0$.

(ii) Bootstrapping Alternatively, a parametric bootstrap approach can be employed: we first simulate noise terms $(\varepsilon, \mu)$ from a bivariate Gaussian distribution with mean $0$, variance $1$, and correlation $\hat \theta$ and add them to signals (with $\hat \alpha(\lambda_1)$ and $\hat \beta(\lambda_2)$) to generate synthetic responses, from which a bootstrap estimate $\hat \theta^{\rm Boot}$ is computed. Repeating this procedure multiple times (say $B$ times) yields a collection of bootstrap estimates, say $\hat \theta\_j^{\rm Boot}$ for $1 \le j \le B$, whose empirical variance can then be used as an estimate of the variance of $\hat \theta$. Finally, to ensure consistency of this variance estimator, a bias-correction step is needed, for which we need to estimate $(||\alpha_0||_2^2, ||\beta_0||_2^2, \alpha_0^\top \beta_0)$ consistently as before. This method circumvents the need to evaluate the complex analytical form of the variance. We would be happy to provide further details in final version.

Finally, we evaluate the estimated variance over a grid of $(\lambda_1, \lambda_2)$ values and select the pair minimizing estimated variance. This would give us an estimator of $\hat \theta$, which is $\sqrt{n}$-consistent and has minimum asymptotic variance among all estimators in our proposed class. This aids in inference via giving a confidence interval with a small width.

Simulation To validate, we run simulations for $\hat{\theta}^{{db, INT}}_{{3sp}}$ using bootstrap. In particular, we compute the ratio of the bootstrap estimator of the limiting standard deviation (s.d.) to that of the gold-standard oracle Monte-Carlo s.d. (obtained from repeated draws from the true data-generating process). We set $p = 250, n = 500, \rho = 0.5, B = 10000$ and a grid of 100 $\lambda$ values from $0.05$ to $10$. We have the summary of this ratio (mean = $1.024$):

This shows that the bootstrap estimator closely matches the true asymptotic variance, and we have already shown (see Figure 2,4, and 6 in the appendix) that the bias is negligible. Thus, the data-driven procedure yields valid coverage.

4. Debiasing function g:

The debiasing functions $g$ are defined using $\triangleq$, e.g., in line 184 $g^{\rm INT}\_{\rm 3sp}(\lambda_1, \lambda_2)= \left( \int \frac{x \, dF_{\rm MP}(x) }{x+\lambda_1}\right) \left( \int \frac{x \, dF_{\rm MP}(x) }{x + \lambda_2}\right),$ and in line 198 $g_{1, \rm 2sp}^{\rm INT}(\lambda_1, \lambda_2) = \int \frac{x^2 \ dF_{\rm MP}(x)}{(x + \lambda_1)(x + \lambda_2)} \,,$ where $F_{\rm MP}$ is Marchenko-Pasteur (MP)distribution, the limiting spectral distribution of eigenvalues of sample covariance matrices. In practice, we calculated by 10000 Monte Carlo samples of MP distribution.

5. Optimality of ridge estimator and sharp asymptotic analysis:

In parameter estimation, optimality typically refers to either i) rate optimality in the minimax sense, and ii) semi-parametric efficiency. As our estimator is proven to be $\sqrt{n}$-consistent, it achieves the former, as it attains the gold-standard parametric estimation rate.

In classical fixed-dimensional settings, semiparametric efficiency typically refers to the smallest achievable asymptotic variance among all regular estimators. In contrast, our analysis operates under the proportional asymptotic regime (e.g., see ref. [2] in our paper), where the covariate dimension $p$ grows proportionally larger than the sample size $n$, i.e., $p/n \to c \in (1, \infty)$, without assuming sparsity in $E[Y|X]$ or $E[A|X]$.

In this high-dimensional context, foundational concepts like regularity of estimators and Hajek–LeCam local minimax theory are not yet well developed. Thus, rather than pursuing classical efficiency bounds, we focus on identifying the best estimator within a given estimation strategy. Specifically, we ask: given a debiasing approach with regularized nuisance estimation and a fixed sample-splitting scheme (e.g., 1/2 or 1/3 splits), which configuration yields the lowest asymptotic variance?

Can that be achieved by merely a priori choosing the level of regularization from a prediction optimality perspective -- as is the case for DML estimators under fixed $p$ or ultra-high dimensional sparse regimes? Our notion of sharp optimality is grounded in this perspective. We argue that when ridge regularization is used to estimate the nuisance functions, the regularization parameter should be chosen by minimizing the asymptotic variance of $\theta_0$ directly, rather than the prediction optimal one (as advocated by DML). That said, it remains possible that alternative regularization schemes, such as LASSO, may yield even lower variance. However, identifying the minimum possible achievable variance in this regime is highly nontrivial and is unknown even for simple linear regression. Therefore, this remains an important open problem for future investigation. We would be happy to include this discussion in the final version and are open to incorporating any other suggestions you may have.

We greatly appreciate the reviewer’s feedback and suggestions, and would like to address them as follows.

1. Linearity assumption and ridge regression:

We agree that the linear model assumption may appear restrictive at first glance. However, we would like to emphasize that our analysis includes the challenging proportional asymptotic regime where the number of covariates $p$ is larger than the sample size $n$, and $p/n \to c \in (1, \infty)$. This regime is fundamentally different from the standard high-dimensional setting, where valid inference typically relies on some form of sparsity assumption. Theoretical understanding of the proportional regime is still in its early stages. In fact, results on simple linear regression [1] and AIPW estimators under the linear outcome regression and propensity score model and $p < n$ (ref. [23] of our paper) have only recently been established. From this perspective, focusing on the linear model is not merely a simplifying assumption, but rather a necessary and meaningful first step toward developing a deeper understanding of more complex models in this regime.

A natural next step is to extend the linear regression setting to regression in a general feature space, such as kernel ridge regression. Since our analysis relies on deterministic equivalents of functionals of Gram matrices, the techniques can be carried over, provided similar results hold for Gram matrices in the feature space setting.

Regarding the ridge penalty, we believe that similar analyses can be conducted in the presence of other convex penalties (e.g., $\ell_1$ norm), borrowing tools from Approximate Message Passing (AMP) theory. Debiasing methods based on AMP have been studied in the context of linear models [1] and doubly robust estimators such as AIPW under the regime $p < n$ of ref. [23]. These developments suggest that similar techniques could be adapted to our setting.

2. Gaussian assumption:

Thank you for raising this concern. We would first like to point out that for the majority of our results (Theorems 3.1 - 3.3, which show that the estimators are $\sqrt{n}$-consistent), we do not need normality of the covariates or the errors. The following two conditions are sufficient for these three estimators: (i) $X_{ij}$ are i.i.d. with mean $0$ and variance $1$ and finite $(4 + \delta)$ moment, and (ii) $(\epsilon, \mu) \overset{i.i.d.}{\sim} F$ where $F$ with mean $0$, variance $1$, correlation $\rho$, and has finite $(2 + \delta)$ moment for some $\delta > 0$.

Regarding allowing dependence between the covariates, we believe that our results can be extended to the settings where $X_i$ is a subgaussian random variable in $\mathbb{R}^p$ with bounded sub-gaussian norm using more modern techniques. However, we leave this as future work, as we do not believe it will add anything substantial to the main message of the paper, and the current analysis is already quite complicated.

Theorem 3.4 indeed relies on the assumption of Gaussian covariates and error distribution to derive the exact form of the limiting variance. That said, we believe these assumptions could potentially be relaxed by appealing to suitable universality results. Nonetheless, doing so would require substantial theoretical development and a careful analysis to rigorously establish the necessary universality results, which is a challenging task in random matrix theory.

We added additional simulations below to show robustness against the Gaussian assumption. See our response to ``Additional simulations" below.

3. Minimaxity:

In this paper, we establish that ECC can be estimated at the $\sqrt{n}$ rate, which is known to be minimax optimal for estimating finite-dimensional parameters under standard assumptions. Consequently, our estimator is minimax rate optimal. Obtaining sharp efficiency bounds beyond rate minimaxity for functional estimation under proportional asymptotics remains an open direction. Therefore, we explore optimal estimation under a given strategy of estimation, i.e., ridge tuning, sample splitting, and given choice of estimators such as INT, NR, or DR.

4. Confidence interval coverage:

In the analyses in the current paper, the bias of the debiased estimators is negligible compared to their standard errors. We also discuss standard error estimation in the following bullet point, which we found to be within 5% of the true standard error. Together, these suggest near-nominal confidence interval coverage. However, we would be happy to report how sensitive confidence intervals are to variance estimation in the revised paper.

5. Practical methods for selecting tuning parameters:

Thank you for raising this concern. While the expressions for the limiting variances may appear complex, they are in fact smooth deterministic functions of three unknown quantities: $(||\alpha_0||_2^2, ||\beta_0||_2^2, \alpha_0^\top \beta_0)$. This allows practitioners to adopt one of two practical strategies, as outlined below:

(i) consistent estimation of unknown quantities One may consistently estimate $(||\alpha_0||\_2^2, ||\beta_0||\_2^2, \alpha\_0^\top \beta\_0)$ and plug them into the variance expressions. To estimate $||\alpha_0||\_2^2$, one can first compute a simple plugged-in estimator $||\hat \alpha(\lambda)||\_2^2$, ($\hat \alpha(\lambda)$ is a ridge estimator of $\alpha\_0$) with tuning parameter $\lambda$, and then apply bias correction to obtain a consistent estimator. Similar procedures apply to $||\beta_0||\_2^2$ and $\alpha\_0^\top \beta\_0$.

(ii) Bootstrapping Alternatively, a parametric bootstrap approach can be employed: we first simulate noise terms $(\varepsilon, \mu)$ from a bivariate Gaussian distribution with mean $0$, variance $1$, and correlation $\hat \theta$ and add them to signals (with $\hat \alpha(\lambda_1)$ and $\hat \beta(\lambda_2)$) to generate synthetic responses, from which a bootstrap estimate $\hat \theta^{\rm Boot}$ is computed. Repeating this procedure multiple times (say $B$ times) yields a collection of bootstrap estimates, say $\hat \theta\_j^{\rm Boot}$ for $1 \le j \le B$, whose empirical variance can then be used as an estimate of the variance of $\hat \theta$. Finally, to ensure consistency of this variance estimator, a bias-correction step is needed, for which we need to estimate $(||\alpha_0||_2^2, ||\beta_0||_2^2, \alpha_0^\top \beta_0)$ consistently as before. This method circumvents the need to evaluate the complex analytical form of the variance. We would be happy to provide further details in the final version.

Finally, we evaluate the estimated variance over a grid of $(\lambda_1, \lambda_2)$ values and select the pair minimizing estimated variance. This would give us an estimator of $\hat \theta$, which is $\sqrt{n}$-consistent and has minimum asymptotic variance among all estimators in our proposed class. This aids in inference via giving a confidence interval with a small width. Due to space constraints, we omit the technical details here, but we would be happy to elaborate upon request and add a discussion in the final version of the paper.

6. Additional simulations

Parametric Bootstrap: We run simulations for $\hat{\theta}^{{\rm INT, db}}_{{\rm 3sp}}$ using bootstrap. In particular, we compute the ratio of the bootstrap estimator of the limiting standard deviation (s.d.) to that of the gold-standard oracle Monte-Carlo s.d. (obtained from repeated draws from the true data-generating process). We set $p = 250, n = 500, \rho = 0.5, B = 10000$ and and a grid of 100 $\lambda$ values from $0.05$ to $10$. Below is the summary of this ratio (mean = $1.024$) over $\lambda$:

This shows that the bootstrap estimator closely matches the true asymptotic variance across different $\lambda$, and we have already shown (see Figure 2,4, and 6 in the appendix) that bias is negligible. Thus, our data-driven procedure yields valid coverage.

Non-gaussian error: We perform simulations with non-Gaussian errors for $\hat{\theta}^{{\rm INT, db}}_{{\rm 3sp}}$, where we now generate errors from a bivariate Laplace distribution with mean $0$, marginal variances $1$, and correlation $\rho=0.5$, same as the Gaussian setting in the paper. The de-biased estimators had negligible bias in this case, and the following table provides a summary of the ratio (Laplace to Gaussian) of s.d. over 10000 Monte Carlo iterations:

This implies our results continue to hold beyond Gaussianity. We would be happy to elaborate upon request and add this to the final version of the paper.

7. Variance blowup:

The variance blow-up is primarily driven by the bias correction factor. For instance, consider the definition of $\hat{\theta}\_{2sp}^{\rm INT, db}$ in Equation (3.4). It is immediate from the definition that the variance of this estimator diverges when $\left(1 - \frac{g\_{2, 2sp}^{\rm INT}(\lambda_1, \lambda_2)}{g\_{1, 2sp}^{\rm INT}(\lambda_1, \lambda_2)}\right) \approx 0 \,.$ This scenario arises in our simulations when $\lambda_1 = \lambda_2 = \lambda$ is around 1.5 (Top left picture of Figure 1), where $g\_{1, 2sp}^{\rm INT}(\lambda, \lambda) \approx g\_{2, 2sp}^{\rm INT}(\lambda, \lambda)$, causing the variance to explode. However, this phenomenon does not affect our method. Since our goal is to select $(\lambda_1, \lambda_2)$ that minimizes the asymptotic variance, such values where the variance is unstable are naturally excluded from consideration.

[1] Hastie, T., Montanari, A., Rosset, S., & Tibshirani, R. J. (2022). Surprises in high-dimensional ridgeless least squares interpolation. Annals of statistics, 50(2), 949.

We greatly appreciate the reviewer’s feedback and suggestions, and would like to address them as follows.

1. The sample-splitting approach appears to sacrifice the sample efficiency of DML:

Thank you for raising this important point. You are absolutely right that sample splitting generally sacrifices statistical efficiency, and cross-fitting is a widely adopted strategy to address this issue. However, our analysis is conducted under the proportional asymptotic regime where $p/n \to c \in (1, \infty)$. In this setting, a rigorous analysis of cross-fitting becomes substantially more challenging. To the best of our knowledge, recent work (ref. [23] of our paper) analyzed cross-fitting in the context of AIPW estimators, but only in the regime $p/n \to c \in (0, 1)$, where the dimensionality is smaller than the sample size. The main technical difficulty is that, in classical fixed-$p$ settings (or in high-dimensional settings with sparsity constraints), the estimators obtained from different sample splits are only weakly correlated and the correlation is asymptotically negligible. However, in the proportional regime, this correlation is often non-negligible, rendering the analysis of cross-fitting significantly more challenging. Moreover, whether cross-fitting effectively reduces variance, i.e., whether the non-negligible correlation between estimators from different splits is negative, can depend on several factors, including the sample splitting scheme, the type of estimator used, and the specific parameter settings. In fact, recent work (ref. [23] of our paper) has observed a striking phenomenon: the correlation between estimators from different splits can be positive and non-vanishing, leading to increased variance of the cross-fitted estimator. This highlights the nuanced and context-dependent behavior of cross-fitting in the proportional asymptotic regime. While we believe it may be possible to adapt the techniques from [23] and extend them carefully to this setting, doing so is non-trivial and would likely require substantial new insights. As such, we leave this as an important direction for future work. We expect that a complete theoretical treatment of cross-fitting under the regime $p/n \to c \in (1, \infty)$ would constitute a contribution in its own right.

2. Extension beyond linear model and ridge:

Thank you for raising this interesting point! While, as a first step, we focus on a linear data-generating model, we agree that extending the analysis to more expressive non-linear models—such as kernel ridge regression (KRR) or random feature (RF) models is an important direction for future work. Recent advances (e.g., [1]) have analyzed KRR and RF estimators under the proportional asymptotic regime, leveraging tools from random matrix theory. Therefore, we believe that it is possible to construct a consistent estimator of $\theta_0$ when the conditional mean function is a polynomial. However, extending these results to enable $\sqrt{n}$-consistent estimation and drawing a valid inference remains challenging. In particular, debiasing non-linear estimators in the proportional regime is a highly non-trivial task, and to the best of our knowledge, no prior work has addressed this issue, even in the context of simple regression. We believe this represents a promising but technically demanding research avenue for future work.

If we replace the squared error loss with another loss function (e.g., Huber loss), it is still possible to leverage similar tools from random matrix theory (RMT) to arrive at the desired conclusions. This is supported by existing work on general M-estimation problems in high dimensions using RMT techniques (e.g., [2]). However, if the penalty function is changed to the $\ell_1$ norm or other convex penalties, the analysis becomes more involved. In such settings, the tools from Approximate Message Passing (AMP), Convex Gaussian minmax theorem (CGMT), or Stein's identity become useful to obtain a $\sqrt{n}$-consistent and asymptotically normal (CAN) estimator of $\theta_0$. Debiasing methods based on these techniques have been studied in the context of linear models (ref. [25] of our paper) and doubly robust estimators such as AIPW under the regime $p < n$ (ref. [23] of our paper). These developments suggest that similar techniques could be adapted to our setting. We leave this direction as a promising avenue for future investigation.

3. Asymptotic variance:

The analytic expression of the variances indicates that given any $(\lambda_1, \lambda_2)$, the limiting variance is a deterministic function of three unknowns: $(||\alpha_0||_2^2, ||\beta_0||_2^2, \alpha_0^\top \beta_0)$. The final goal is to choose $(\lambda_1, \lambda_2)$ which minimizes these variance expressions. This also echoes the main point of our paper: the prediction-optimal tuning parameters, i.e., the value of $\lambda_1$ (resp. $\lambda_2$) that yields the minimum prediction risk of $\alpha_0$ (resp. $\beta_0$), may not yield an estimator of ECC with minimum asymptotic variance. This contrasts with the conventional wisdom in double machine learning (DML), where nuisance parameters are typically estimated by minimizing prediction error, a point which is further highlighted in our simulation section. We next present two practical strategies for consistently estimating the limiting variance and selecting optimal tuning parameters:

(i) consistent estimation of unknown quantities One may consistently estimate $(||\alpha_0||\_2^2, ||\beta_0||\_2^2, \alpha\_0^\top \beta\_0)$ and plug them into the variance expressions. To estimate $||\alpha_0||\_2^2$, one can first compute a simple plugged-in estimator $||\hat \alpha(\lambda)||\_2^2$, ($\hat \alpha(\lambda)$ is a ridge estimator of $\alpha\_0$) with tuning parameter $\lambda$, and then apply bias correction to obtain a consistent estimator. Similar procedures apply to $||\beta_0||\_2^2$ and $\alpha\_0^\top \beta\_0$.

(ii) Bootstrapping Alternatively, a parametric bootstrap approach can be employed: we first simulate noise terms $(\varepsilon, \mu)$ from a bivariate Gaussian distribution with mean $0$, variance $1$, and correlation $\hat \theta$ and add them to signals (with $\hat \alpha(\lambda_1)$ and $\hat \beta(\lambda_2)$) to generate synthetic responses, from which a bootstrap estimate $\hat \theta^{\rm Boot}$ is computed. Repeating this procedure multiple times (say $B$ times) yields a collection of bootstrap estimates, say $\hat \theta\_j^{\rm Boot}$ for $1 \le j \le B$, whose empirical variance can then be used as an estimate of the variance of $\hat \theta$. Finally, to ensure consistency of this variance estimator, a bias-correction step is needed, for which we need to estimate $(||\alpha_0||_2^2, ||\beta_0||_2^2, \alpha_0^\top \beta_0)$ consistently as before. This method circumvents the need to evaluate the complex analytical form of the variance. We would be happy to provide further details in final version.

Finally, we evaluate the estimated variance over a grid of $(\lambda_1, \lambda_2)$ values and select the pair minimizing estimated variance. This would give us an estimator of $\hat \theta$, which is $\sqrt{n}$-consistent and has minimum asymptotic variance among all estimators in our proposed class. This aids in inference via giving a confidence interval with a small width.

Our simulation suggests that the bootstrap estimator works well, but we cannot add it due to space constraints. Kindly see our response to Comment 3 of Reviewer FNJX and Comment 6 of AoFP for the details.

4. Variance blowup:

The variance blow-up is primarily driven by the bias correction factor. For instance, consider the definition of $\hat{\theta}\_{2sp}^{\rm INT, db}$ in Equation (3.4). It is immediate from the definition that the variance of this estimator diverges when $\left(1 - \frac{g\_{2, 2sp}^{\rm INT}(\lambda_1, \lambda_2)}{g\_{1, 2sp}^{\rm INT}(\lambda_1, \lambda_2)}\right) \approx 0 \,.$ This scenario arises in our simulations when $\lambda_1 = \lambda_2 = \lambda$ is around 1.5 (Top left picture of Figure 1), where $g\_{1, 2sp}^{\rm INT}(\lambda, \lambda) \approx g\_{2, 2sp}^{\rm INT}(\lambda, \lambda)$, causing the variance to explode. However, this phenomenon does not affect our method. Since our goal is to select $(\lambda_1, \lambda_2)$ that minimizes the asymptotic variance, such values where the variance is unstable are naturally excluded from consideration.

[1] Spectrum of inner-product kernel matrices in the polynomial regime and multiple descent phenomenon in kernel ridge regression. Theodor Misiakiewicz.

[2] Asymptotics For High Dimensional Regression M-Estimates: Fixed Design Results. Lihua Lei, Peter J. Bickel, Noureddine El Karoui.

Thank you very much for your kind words! If you have any other comments/concerns during the discussion period, kindly do let us know. We would be happy to answer them.