Control variate selection for Monte Carlo integration

TL;DR

This paper explores the use of regularized control variate selection via Lasso in Monte Carlo integration, improving efficiency and stability when using many control variates.

Contribution

It introduces a Lasso-based regularization method for selecting control variates in Monte Carlo integration, enhancing accuracy and computational stability.

Findings

Lasso regularization improves estimator stability.

Using many control variates increases efficiency.

Concentration inequalities confirm numerical results.

Abstract

Monte Carlo integration with variance reduction by means of control variates can be implemented by the ordinary least squares estimator for the intercept in a multiple linear regression model with the integrand as response and the control variates as covariates. Even without special knowledge on the integrand, significant efficiency gains can be obtained if the control variate space is sufficiently large. Incorporating a large number of control variates in the ordinary least squares procedure may however result in (i) a certain instability of the ordinary least squares estimator and (ii) a possibly prohibitive computation time. Regularizing the ordinary least squares estimator by preselecting appropriate control variates via the Lasso turns out to increase the accuracy without additional computational cost. The findings in the numerical experiment are confirmed by concentration inequalities for the integration error.

Tables

Table 1. Table 1 : Number of control variates m 𝑚 m by degree threshold 𝑑𝑒𝑔 𝑑𝑒𝑔 \mathit{deg} in dimension d 𝑑 d constructed out of tensor products of k 𝑘 k univariate polynomials.

$d$	$k$	Degree threshold ($\mathrm{𝑑𝑒𝑔}$)
		1	3	5	10	12
3	12	3	19	55	285	454
5	10	5	55	251	3 001	6 157
8	3	8	164	1 214	20 993	36 813

Table 2. Table 2 : Sample sizes n 𝑛 n and sub-sample sizes N 𝑁 N together with the range [ c 1 ​ n , c 2 ​ n ] subscript 𝑐 1 𝑛 subscript 𝑐 2 𝑛 [c_{1}\sqrt{n},c_{2}\sqrt{n}] corresponding to the imposed number of selected control variates in LSLASSO.

$n$	$N$	$\lfloor 3\sqrt{n}\rfloor$	$\lfloor 12\sqrt{n}\rfloor$
02 000	0 700	134	0 536
05 000	1 000	212	0 848
10 000	2 000	300	1 200

Table 3. Table 3 : Statistical efficiency for φ 𝜑 \varphi ; see also Figure  1(a) .

$m=$	3	19	55	285	454
OLS	8.42e00	8.56e02	2.12e05	2.49e11	5.27e14
LASSO	8.42e00	8.53e02	6.72e04	7.71e04	7.71e04
LSL	8.42e00	8.58e02	2.10e05	6.26e05	1.37e06
LSLX	8.42e00	8.51e02	2.09e05	2.49e11	2.91e05
QMC	Halton: 8.76e01   Sobol: 3.29e02

Table 4. Table 7 : European dipper capture-recapture data (Marzolin, 1988 ) . The counts in the triangle refer to the number of birds released in a given year and recaptured for the first time in a later year.

Release year	Birds released	Year of recapture: $1981+\cdots$
		1	2	3	4	5	6
1981	22	11	2	0	0	0	0
1982	60		24	1	0	0	0
1983	78			34	2	0	0
1984	80				45	1	2
1985	88					51	0
1986	98						52

Table 5. Table 8 : Capture data: statistical efficiency ( n = 2000 𝑛 2000 n=2000 )

$m=$	90	444	1 062	3 090	5 730
OLS	9.33	20.7	14.7	0.14	0.06
LASSO	9.34	20.3	16.7	14.4	8.57
LSL	9.33	20.4	12.8	8.43	4.60
LSLX	9.33	19.4	19.8	12.9	7.86

Equations

$$P(f)={\int }_{\mathcal{X}}f(x)P(dx).$$

$${\hat{\alpha }}_n^{\mathrm{mc}}(f)=P_n(f)=\frac{1}{n}\sum _{i=1}^{n}f(X_i).$$

$${\beta }^{\star }(f)\in \underset{\beta \in {\mathbb{R}}^m}{\mathrm{argmin}}P[(f-P(f)-{\beta }^{T}h{)}^2]$$

$${\hat{\alpha }}_n^{\mathrm{or}}(f)=P_n[f-{\beta }^{\star }(f{)}^{T}h].$$

$$(P(f),{\beta }^{\star }(f))\in \underset{(\alpha ,\beta )\in \mathbb{R}\times {\mathbb{R}}^m}{\mathrm{argmin}}P[(f-\alpha -{\beta }^{T}h{)}^2].$$

$${\mathcal{R}}_n(\alpha ,\beta )={\Vert f^{(n)}-\alpha \text{mathds}{1}_n-H\beta \Vert }_2^2$$

$$\bigl{(}\hat{\alpha}_{n}^{\mathrm{ols}}(f),\hat{\beta}_{n}^{\mathrm{ols}}(f)\bigr{)}\in\operatornamewithlimits{\arg\min}_{(\alpha,\beta)\in\mathbb{R}\times\mathbb{R}^{m}}\mathcal{R}_{n}(\alpha,\beta)$$

$$H=\bigl{(}h_{j}(X_{i})\bigr{)}_{\begin{subarray}{c}i=1,\ldots,n\\ j=1,\ldots,m\end{subarray}}.$$

$$\{\begin{array}{rl}{\hat{\alpha }}_n^{\mathrm{ols}}(f)& =P_n[f-{\hat{\beta }}_n^{\mathrm{ols}}(f{)}^{T}h],\\ {\hat{\beta }}_n^{\mathrm{ols}}(f)& \in \underset{\beta \in {\mathbb{R}}^m}{\mathrm{argmin}}{\Vert f_c^{(n)}-H_c\beta \Vert }_2^2,\end{array}$$

$$\bigl{(}\hat{\alpha}_{n}^{\mathrm{lasso}}(f),\hat{\beta}_{n}^{\mathrm{lasso}}(f)\bigr{)}\in\operatornamewithlimits{\arg\min}\limits_{(\alpha,\beta)\in\mathbb{R}\times\mathbb{R}^{m}}\frac{1}{2n}\mathcal{R}_{n}(\alpha,\beta)+\lambda\left\lVert{\textstyle\beta}\right\rVert_{1}$$

$$\{\begin{array}{rl}{\hat{\alpha }}_n^{\mathrm{lasso}}(f)& =P_n[f-{\hat{\beta }}_n^{\mathrm{lasso}}(f{)}^{T}h],\\ {\hat{\beta }}_n^{\mathrm{lasso}}(f)& \in \underset{\beta \in {\mathbb{R}}^m}{\mathrm{argmin}}\frac{1}{2n}{\Vert f_c^{(n)}-H_c\beta \Vert }_2^2+\lambda {\Vert \beta \Vert }_1.\end{array}$$

$$\bigl{(}\hat{\alpha}_{n}^{\mathrm{lslasso}}(f),\hat{\beta}_{n}^{\mathrm{lslasso}}(f)\bigr{)}\in\operatornamewithlimits{\arg\min}_{(\alpha,\beta)\in\mathbb{R}\times\mathbb{R}^{\hat{\ell}}}\left\lVert{\textstyle f^{(n)}-\alpha\mathds{1}_{n}-H^{(n)}_{\hat{S}_{N}}\beta}\right\rVert_{2}^{2}$$

$$\left|{\hat{\alpha }}_n^{\mathrm{or}}(f)-P(f)\right|⩽\sqrt{2\log (2/\delta )}\frac{\tau }{\sqrt{n}}.$$

$$B=\underset{x\in \mathcal{X}}{\sup }h(x{)}^{\top }{G}^{-1}h(x)=\underset{x\in \mathcal{X}}{\sup }\hslash (x{)}^{\top }\hslash (x),$$

$$n\geqslant\max\bigl{(}18B\log(4m/\delta),\;75m\log(4/\delta)\bigr{)}$$

$$\begin{gathered} \left|{\hat{\alpha }}_n^{\mathrm{ols}}(f)-P(f)\right|⩽\sqrt{2\log (8/\delta )}\frac{\tau }{\sqrt{n}} \\ +58\sqrt{Bm\log (8m/\delta )\log (4/\delta )}\frac{\tau }{n}. \end{gathered}$$

$$m⩽B⩽\underset{x\in \mathcal{X}}{\sup }{h}^{\top }(x)h(x)/\gamma ⩽m{U}_h^2/\gamma .$$

$$S^{\star }:=\mathrm{supp}({\beta }^{\star }(f)),$$

$$\mathcal{C}(S;\alpha )=\{u\in {\mathbb{R}}^m:{\Vert u_{\overline{S}}\Vert }_1⩽\alpha {\Vert u_S\Vert }_1\}.$$

$$n$$

$$\lambda$$

$$\begin{gathered} \left|{\hat{\alpha }}_n^{\mathrm{lasso}}(f)-P(f)\right|⩽\sqrt{2\log (8/\delta )}\frac{\tau }{\sqrt{n}} \\ +68\lambda {\ell }^{\star }\sqrt{\log (8m/\delta )}\frac{U_h/{\gamma }^{\star }}{\sqrt{n}}. \end{gathered}$$

$$\begin{gathered} \left|{\hat{\alpha }}_n^{\mathrm{lasso}}(f)-P(f)\right|⩽\sqrt{2\log (8/\delta )}\frac{\tau }{\sqrt{n}} \\ +476{\ell }^{\star }\log (8m/\delta )(U_h^2/{\gamma }^{\star })\frac{\tau }{n}. \end{gathered}$$

$$n⩾70({\ell }^{\star }{U}_h^2/{\gamma }^{\star \star }{)}^2\log (10{\ell }^{\star }m/\delta ),$$

$$13U_h\sqrt{\log (10m/\delta )}\frac{\tau }{\sqrt{n}}⩽\lambda ⩽\frac{{\gamma }^{\star \star }}{3\sqrt{{\ell }^{\star }}}\underset{k\in S^{\star }}{\min }|{\beta }_k^{\star }(f)|,$$

$$B^{\star }⩽{\lambda }_{\max }(G_{S^{\star }}^{-1})\underset{x\in \mathcal{X}}{\sup }{h}_{S^{\star }}^T(x)h_{S^{\star }}(x)⩽{\ell }^{\star }{U}_h^2/{\gamma }^{\star \star }.$$

$$N⩾75{{\xi }^{\star }}^2\log (20{\ell }^{\star }m/\delta ),$$

$$13U_h\sqrt{\log (20m/\delta )}\frac{\tau }{\sqrt{N}}⩽\lambda ⩽\frac{{\gamma }^{\star \star }}{3\sqrt{{\ell }^{\star }}}\underset{k\in S^{\star }}{\min }|{\beta }_k^{\star }(f)|,$$

$$\begin{gathered} \left|{\hat{\alpha }}_n^{\mathrm{lslasso}}(f)-P(f)\right|⩽\sqrt{2\log (16/\delta )}\frac{\tau }{\sqrt{n}} \\ +58\sqrt{B^{\star }{\ell }^{\star }\log (16{\ell }^{\star }/\delta )\log (8/\delta )}\frac{\tau }{n}. \end{gathered}$$

$$\begin{gathered} \left|{\hat{\alpha }}_n^{\mathrm{ols}}(f)-P(f)\right|⩽\sqrt{2\log (8/\delta )}\frac{\tau }{\sqrt{n}} \\ +83m\sqrt{\log (8m/\delta )\log (4/\delta )}\frac{\tau }{n} \end{gathered}$$

$$\begin{gathered} \left|{\hat{\alpha }}_n^{\mathrm{lslasso}}(f)-P(f)\right|⩽\sqrt{2\log (16/\delta )}\frac{\tau }{\sqrt{n}} \\ +83{\ell }^{\star }\sqrt{\log (16{\ell }^{\star }/\delta )\log (8/\delta )}\frac{\tau }{n}. \end{gathered}$$

Full text

1 Introduction

Whereas the basic Monte Carlo (MC) estimate of an integral or expectation is given by $(1/n)\sum_{i}f_{i}$, for independent and identically distributed random variables $f_{i}$, the control variates method is based on $(1/n)\sum_{i}(f_{i}+h_{i})$, where the variables $h_{i}$, called control variates, are constructed to have zero expectation. When the controls $h_{i}$ have been selected or estimated properly (based on the samples $f_{i}$), the use of control variates might reduce the variance of the basic MC estimate significantly. The method of control variates, already used frequently to compute prices of financial derivatives Glasserman (2013); Gobet and Labart (2010), has been employed recently in many different fields of Machine Learning and Statistics. Examples include (i) reinforcement learning and more particularly policy gradient methods Jie and Abbeel (2010); Liu et al. (2018) where the score function permits to define many control variates; (ii) inference in complex probabilistic models (Ranganath et al., 2014) where the Stein method allows to define accurate control variates (see e.g., (Oates et al., 2017; Brosse et al., 2018; Belomestny et al., 2020) and the references therein); (iii) gradient based optimization (Wang et al., 2013; Gower et al., 2018), (iv) time series analysis when approximating the characteristic function (Davis et al., 2019), and (v) semi-supervised inference (Zhang et al., 2019).

Suppose that $m\geqslant 1$ control variates are available and $n\geqslant 1$ samples have been generated. Any linear combination of control variates can be used as a particular control variate. In terms of the variance of the estimation error, the optimal linear combination can be estimated based on the empirical risk minimization principle applied to an ordinary least squares (OLS) regression problem [see Eq. (3) below]. This approach, referred to as OLS, is the most common implementation of the control variates method as detailed for instance in (Owen, 2013, Section 8.3) or (Portier and Segers, 2019; South et al., 2018), although other implementations are possible, see Remark 2 below.

Asymptotically, the OLS error is bounded by the MC error and is proportional to the $L_{2}$ approximation error of the integrand in the linear span of control variates (Glynn and Szechtman, 2002). In combination with well-known approximation results in $L_{p}$-spaces (Rudin, 2006), this representation of the OLS error suggests to use an increasing number of control variates. Indeed, in (Portier and Segers, 2019) it is shown that when $m$ grows with $n$, the OLS error rate can be faster than $1/\sqrt{n}$.

However, when based on a large number of control variates, the OLS suffers from two classical problems common for least squares methods: (i) numerical instabilities when the control variates are nearly collinear, and (ii) a computational complexity in $m^{3}+nm^{2}$, which might be prohibitive.

To deal with these two issues, it has been proposed in South et al. (2018) to regularize the OLS estimate by adding a $\ell_{1}$-penalty term in the minimization problem, just as in the LASSO (Tibshirani, 1996). Simulation results in South et al. (2018) show that this approach, referred to as LASSO, provides great improvements in practice. However, those practical findings are not supported by an asymptotic error rate nor by a non-asymptotic error bound.

The main objective of the paper is to provide a non-asymptotic theory for the use of control variates in Monte Carlo simulations. The contributions are as follows.

A new method called LSLASSO is proposed. In the spirit of (Belloni and Chernozhukov, 2013), it consists in selecting the best control variates via the LASSO, using subsampling to decrease the computation time, and then to apply OLS with the selected controls. 2. 2.

Support recovery: the LASSO is shown to select the correct control variates with large probability. 3. 3.

Concentration inequalities are derived for the OLS, LASSO and LSLASSO integration errors. The one for the OLS highlights a compromise between the approximation error of the integrand in the linear span of control variates and the multicollinearities between the control variates. The ones for (LS)LASSO show significant improvements regarding the effects of multicollinearity.

The approach for the proofs combines well known sub-Gaussian concentration inequalities (Boucheron et al., 2013) along with a lower bound for the smallest eigenvalue of an empirical Gram matrix, based on a Chernoff inequality for matrices (Tropp, 2015, Theorem 5.1.1).

The outline of the paper is as follows. Section 2 introduces the theoretical background and the different MC estimates and provides some comments about their practical implementation and some possible alternative approaches. Section 3 contains the statements of the theoretical results. Sections 4 and 5 describe numerical experiments on artificial and real data to illustrate the finite-sample behavior of the methods. Section 6 concludes the main part of the paper with a discussion of avenues for further research. Section A contains some auxiliary results, whereas the proofs of the four theorems stated in Section 3 are given in Sections B to E.

2 Monte Carlo integration and control variates

Background.

Let $f\in L_{2}(P)$ be a square integrable, real-valued function on a probability space $(\mathcal{X},\mathcal{A},P)$ of which we would like to calculate the integral

[TABLE]

The MC estimator of $P(f)$ based on independent random variables $X_{1},\ldots,X_{n}$ taking values in $\mathcal{X}$ and with common distribution $P$ is

[TABLE]

This estimator is unbiased and has variance $n^{-1}\sigma_{0}^{2}(f)$, where $\sigma_{0}^{2}(f)=P[(f-P(f))^{2}]$.

The control variates are functions $h_{1},\ldots,h_{m}\in L_{2}(P)$ with known expectations. Without loss of generality, assume that $P(h_{k})=0$ for all $k\in\{1,\ldots,m\}$. Let $h=(h_{1},\ldots,h_{m})^{T}$ denote the $\mathbb{R}^{m}$-valued function with the $m$ control variates as elements. Let $\mathcal{F}_{m}=\operatorname{Span}\{h_{1},\ldots,h_{m}\}=\{\beta^{T}h:\beta\in\mathbb{R}^{m}\}$ denote the closed linear subspace of $L_{2}(P)$ generated by the control variates.

For any coefficient vector $\beta=(\beta_{1},\ldots,\beta_{m})^{T}\in\mathbb{R}^{m}$, we have $P(f-\beta^{T}h)=P(f)$, so that $P_{n}(f-\beta^{T}h)$ is an unbiased estimator of $P(f)$, with variance $n^{-1}P[(f-P(f)-\beta^{T}h)^{2}]$. Any oracle coefficient

[TABLE]

minimizes the variance. If such a $\beta^{\star}(f)$ would be known, the resulting oracle estimator would be {ceqn}

[TABLE]

By definition, the oracle estimator achieves the minimal variance $n^{-1}\sigma_{m}^{2}(f)$ where $\sigma_{m}^{2}(f)$ is the minimum value of the variance term $\ P[(f-P(f)-\beta^{T}h)^{2}]$ with respect to $\beta$. For any $m^{\prime}\in\{0,1,\ldots,m\}$, if we use only the first $m^{\prime}$ control variates $h_{1},\ldots,h_{m^{\prime}}$, or even none at all in case $m^{\prime}=0$, we have $\sigma_{m}^{2}(f)\leqslant\sigma_{m^{\prime}}^{2}(f)$. In particular, if $\beta^{\star}(f)$ would be known, the use of control variates would always reduce the variance of the basic Monte Carlo estimator.

As $\beta^{\star}(f)^{T}h$ is the $L_{2}(P)$-projection of $f-P(f)$ on the linear vector space $\mathcal{F}_{m}$ and since the control variates are centered, $\beta^{\star}(f)$ satisfies the normal equations $P(hh^{T})\beta^{\star}(f)=P(hf)$. The integral $P(f)$ thus appears as the intercept of a linear regression model with response $f$ and explanatory variables $h_{1},\ldots,h_{m}$, and it can be expressed as {ceqn}

[TABLE]

The empirical risk minimization paradigm applied to the risk function on the right-hand side of (2) will lead to the OLS and LASSO estimates, to be defined further in this section. The same paradigm suggests the use of other regression methods for MC integration such as Principal Component Regression (PCR) or Ridge Regression, which will not be considered in this paper.

Remark 1 (Choice of control variates)

Which control variates work well depends on the problem. In the Black–Scholes model, for instance, an effective control variate for the price of an option is the geometric average of the price series (Glasserman, 2013, Example 4.1.2)). Two generic ways to construct control variates are to be noted. Whenever $P(dx)=w(x)\,Q(dx)$, where $w:\mathcal{X}\to[0,\infty)$ and $Q$ is a probability measure on $(\mathcal{X},\mathcal{A})$, the quantity of interest is $P(f)=Q(wf)$, so that we can use control variates for $wf$ with respect to $Q$. This trick can be useful in combination with importance sampling (Owen and Zhou, 2000). If $P$ has density $p$ with respect to the Lebesgue measure and if we have access to the derivatives of $p$, Stein’s method might be used to build infinitely many control functions (Oates et al., 2017).

Ordinary Least Squares Monte Carlo.

Replacing the distribution $P$ by the sample measure $P_{n}$ in (2), we obtain the OLS estimator $\hat{\alpha}_{n}^{\mathrm{ols}}(f)$ of $P(f)$ as a minimizer of the empirical risk {ceqn}

[TABLE]

given by {ceqn}

[TABLE]

where $\left\lVert{\textstyle\,\cdot\,}\right\rVert_{2}$ denotes the Euclidean norm, $\mathds{1}_{n}=(1,\ldots,1)^{T}\in\mathbb{R}^{n}$, $f^{(n)}=(f(X_{1}),\ldots,f(X_{n}))^{T}\in\mathbb{R}^{n}$ and $H$ is the random $n\times m$ matrix defined by {ceqn}

[TABLE]

The minimization problem in (3) can be expressed using an OLS estimate with centered variables as {ceqn}

[TABLE]

where $f_{c}^{(n)}=f^{(n)}-\mathds{1}_{n}(\mathds{1}_{n}^{T}f^{(n)})/n$ and $H_{c}=H-\mathds{1}_{n}(\mathds{1}_{n}^{T}H)/n$. Indeed, for fixed $\beta\in\mathbb{R}^{m}$, the minimizer over $\alpha\in\mathbb{R}$ of the objective function in (3) is just $P_{n}(f-\beta^{T}h)=P_{n}(f)-\beta^{T}P_{n}(h)$, and since $P_{n}(f)=(\mathds{1}_{n}^{T}f^{(n)})/n$ and $P_{n}(h)=(\mathds{1}_{n}^{T}H)/n$, the equivalence of (3) and (4) follows.

Remark 2 (Variations)

The solution of the linear regression problem (4) involves the empirical covariance matrix defined by $n^{-1}H_{c}^{T}H_{c}=P_{n}(hh^{T})-P_{n}(h)P_{n}(h^{T})$. Using different estimates of the Gram matrix $P(hh^{T})$ leads to alternative control variate MC estimates for $P(f)$ (Glynn and Szechtman, 2002; Portier and Segers, 2019). For fixed $m$ and as $n\to\infty$, all these estimators are consistent and asymptotically normal. The OLS estimator, however, is the only one that can integrate both the constant functions and the control functions without error.

Remark 3 (Invariance)

The OLS estimator does not change if we replace the control variate vector $h$ by $Ah$, where $A$ is an arbitrary invertible $m\times m$ matrix. Provided the control functions are linearly independent, the property of isotropy, i.e., $P(hh^{T})=I_{m}$, can therefore always be enforced by an appropriate linear transformation of the vector of control variates.

Remark 4 (Computation time)

The computation time of the OLS method is of the order $nm^{2}+m^{3}+nt$, where $nm^{2}$ and $m^{3}$ operations are needed for computing and inverting $H_{c}^{T}H_{c}$ respectively and where $t$ stands for the time needed to evaluate $f$. Computational benefits occur when there are multiple integrands, since the OLS estimate can be represented as $w^{T}f^{(n)}$, where the weight vector $w\in\mathbb{R}^{n}$ does not depend on the integrands Portier and Segers (2019). If $q$ integrals need to be evaluated, the computing time becomes $nm^{2}+m^{3}+qnt$, since the matrix $H_{c}^{T}H_{c}$ only depends on the control variates but not on the integrand.

Remark 5 (Variance reduction)

The advantage of using a given set of $m$ control variates over standard MC can be assessed through the value of the residual standard deviation $\sigma_{m}(f)$. In Portier and Segers (2019), bounds for $\sigma_{m}(f)$ are computed in specific examples. For instance, if $\mathcal{X}=[-1,1]^{d}$ and the $h_{k}$ are tensor products of Legendre polynomials, then for any $k$-times continuously differentiable function $f$ it holds that $\sigma_{m}(f)=O(m^{-k/d})$ as $m\to\infty$. This bound emphasizes the benefits of using polynomials when the integrand is regular.

LASSO Monte Carlo.

The LASSO, introduced in Tibshirani (1996), is a regression technique that consists in minimizing the usual least squares loss plus an $\ell_{1}$-penalty term on the vector of regression coefficients. In contrast with OLS, the LASSO usually produces a vector with many zero coefficients, meaning that the corresponding variables are no longer included in the predictive model. The LASSO thus achieves estimation and variable selection at the same time. As the use of control variates in MC integration is linked with regression, the LASSO can take advantage from situations where many control variates are present but not all of them are useful.

The LASSO estimator $\hat{\alpha}_{n}^{\mathrm{lasso}}(f)$ of $P(f)$ follows from adding a $\ell_{1}$-penalization to the objective function in (3). It is formally defined as {ceqn}

[TABLE]

where $\left\lVert{\textstyle\,\cdot\,}\right\rVert_{1}$ denotes the $\ell_{1}$-norm on Euclidean space. By the same argument used to justify the equivalence of (3) and (4), the LASSO can be based on centered variables via {ceqn}

[TABLE]

Remark 6 (Computation)

For the practical implementation of the LASSO, it is commonly recommended to first center and rescale the explanatory variables empirically (Tibshirani et al., 2015, section 2.2). The centering by the sample mean is taken care of in (5). However, for ease of presentation, no empirical rescaling of the control variates is considered in the theoretical analysis. This is in line with the approach proposed in (Tibshirani et al., 2015, Chapter 11). Still, such rescaling is done in the simulation experiments reported in Section 4.

Remark 7 (Computation time)

The LASSO solution is usually computed approximately by cyclical coordinate descent. At each iteration, this algorithm minimizes (5) with respect to a single coordinate, say $\beta_{k}$, while considering other coordinates, $\beta_{(-k)}\in\mathbb{R}^{m-1}$, as constant. This one-dimensional optimization problem has an explicit argmin. Let $H_{c,k}$ be the $k$-th column of $H_{c}$ that has been normalized such that $\|H_{c,k}\|_{2}=1$ (as indicated in the previous remark). The argmin is then simply given by $\eta_{\lambda}(\langle z_{k},H_{c,k}\rangle)$ where $z_{k}=f_{c}^{(n)}-H_{c,(-k)}\beta_{(-k)}$, $H_{c,(-k)}$ is obtained by removing $H_{c,k}$ from $H_{c}$ and $\eta$ is the soft-thresholding function (Tibshirani et al., 2015, Section 2.4, Eq. (2.14)). Since $n$ operations are needed to update $z_{k}$ and the same number is needed to compute the scalar product, the LASSO requires only $nD+nt$ operations, where $D$ stands for the number of iterations conducted in the cyclical coordinate descent and $t$ represents the time needed to evaluate $f$. The value of $D$ is often imposed by a stopping rule within the algorithm but it could also be fixed by the user in order to control the computing time. The selection of the next coordinate $k$ to update can be done cyclically or at random.

LSLASSO Monte Carlo.

The application of ordinary least squares after model selection by the LASSO has been recently studied in Belloni and Chernozhukov (2013). They show, in the setting of nonparametric regression, that OLS post-LASSO, which is also known under the name LSLASSO, performs better than the LASSO in terms of rate of convergence. Motivated by this result we propose to first use the LASSO to select the active variables among a large number of control variates and then to compute the OLS estimate using only the variables selected at the previous stage. We refer to this approach as the LSLASSO. To decrease the computation time when the dimensions involved in the problem, either $n$ or $m$, are large, we recommend to use sub-sampling of size $N$ smaller than $n$ when conducting the first step.

The active set associated to the coefficient $\beta\in\mathbb{R}^{m}$ is $\operatorname{supp}(\beta)=\{j=1,\ldots,m\,:\,\beta_{j}\neq 0\}$. Let $\hat{S}_{N}=\operatorname{supp}(\hat{\beta}_{N}^{\mathrm{lasso}}(f))$ denote the active set of control variates based on the LASSO coefficient vector defined as in (5) but using only the first $N$ random variables $X_{1},\ldots,X_{N}$ generated. The LSLASSO estimate $\hat{\alpha}_{n}^{\mathrm{lslasso}}(f)$ of $P(f)$ is then defined as the OLS estimate in (3) based on the full sample $X_{1},\ldots,X_{n}$ but using only the control variates $h_{j}$ restricted to $j\in\hat{S}_{N}$, that is, {ceqn}

[TABLE]

where $H^{(n)}_{\hat{S}_{N}}$ is the $n\times\hat{\ell}$ matrix $(h_{j}(X_{i}))_{i=1,\ldots,n,\,j\in\hat{S}_{N}}$ and $\hat{\ell}$ is the cardinality of $\hat{S}_{N}$.

Remark 8 (Computation time)

The number of operations needed for the LSLASSO is of the order $ND+n\hat{\ell}^{2}+\hat{\ell}^{3}+nt$, combining the cost of selecting the control variates on the subsample of size $N$ via cyclical coordinate descent as in Remark 7 and running the OLS estimate based on the selected control variates for the full sample of size $n$ as in Remark 4.

3 Non-asymptotic bounds

To derive concentration inequalities for the errors of the estimators proposed in Section 2, we use the notion of sub-Gaussianity as defined for instance in (Boucheron et al., 2013, Section 2.3). Recall that the moment generating function of a centered Gaussian random variable with variance $\sigma^{2}$ is equal to $\lambda\mapsto\exp(\lambda^{2}\sigma^{2}/2)$.

Definition 1

A centered random variable $Y$ is sub-Gaussian with variance factor $\tau^{2}>0$, notation $Y\in\mathcal{G}(\tau^{2})$, if and only if $\log\operatorname{\mathbb{E}}[\exp(\lambda Y)]\leqslant\lambda^{2}\tau^{2}/2$ for all $\lambda\in\mathbb{R}$.

If $Y\in\mathcal{G}(\tau^{2})$, then necessarily $\operatorname{Var}(Y)\leqslant\tau^{2}$ (Boucheron et al., 2013, Exercise 2.16). Chernoff’s inequality provides exponential bounds on the tails of sub-Gaussian random variables. Moreover, the sum of independent sub-Gaussian variables is again sub-Gaussian. Centered, bounded random variables taking values in an interval $[a,b]$ are sub-Gaussian with variance factor at most $(b-a)^{2}/4$ (Boucheron et al., 2013, Lemma 2.2).

The concentration inequalities for the various Monte Carlo methods with control variates will be largely due to the following assumption that requires the residuals to be sub-Gaussian.

Assumption 1** (Sub-Gaussian residuals)**

The residual function $\epsilon=f-P(f)-\beta^{\star}(f)^{T}h$ satisfies $\epsilon\in\mathcal{G}(\tau^{2})$ for some $\tau>0$, that is, $\int_{\mathcal{X}}\exp\{\lambda\epsilon(x)\}\,P(dx)\leqslant\exp(\lambda^{2}\tau^{2}/2)$ for all $\lambda\in\mathbb{R}$.

The estimation error of the oracle estimator in (1) is just $\hat{\alpha}_{n}^{\mathrm{or}}(f)-P(f)=P_{n}(\epsilon)=n^{-1}\sum_{i=1}^{n}\epsilon(X_{i})$. Under Assumption 1, this is a sub-Gaussian variable with variance factor $\tau^{2}/n$. Chernoff’s inequality (Boucheron et al., 2013, p. 25) then implies that for all $\delta\in(0,1)$ and all integer $n\geqslant 1$, with probability at least $1-\delta$, {ceqn}

[TABLE]

This concentration inequality provides a baseline when the best possible control variate in the space $\mathcal{F}_{m}$ is selected. The case $m=0$ also covers the basic MC method: in that case, $\tau^{2}$ is the variance factor of the sub-Gaussian variable $f-P(f)$ on $(\mathcal{X},\mathcal{A},P)$.

Assumption 2** (Bounded control variates)**

The control variates $h_{1},\ldots,h_{m}\in L_{2}(P)$ are uniformly bounded. Put $U_{h}:=\max_{j=1,\ldots,m}\sup_{x\in\mathcal{X}}\left\lvert{\textstyle h_{j}(x)}\right\rvert$.

For a symmetric real matrix $A$, let $\lambda_{\min}(A)$ and $\lambda_{\max}(A)$ denote its smallest and largest eigenvalues, respectively.

Assumption 3** (Linear independence of control variates)**

The control variates $h_{1}$, $\ldots$, $h_{m}\in L_{2}(P)$ are linearly independent. As a consequence, the $m\times m$ Gram matrix $G:=P(hh^{T})$ is positive definite and its smallest eigenvalue $\gamma:=\lambda_{\min}(G)$ is positive.

Consider the ortho-normalized vector of control variates $\hbar=(\hbar_{1},\ldots,\hbar_{m})^{\top}=G^{-1/2}h$ and put {ceqn}

[TABLE]

a finite quantity by Assumptions 2 and 3. The error OLS estimation error is subject to the following concentration bound.

Theorem 3.1 (Concentration inequality for OLS)

Suppose Assumptions 1, 2 and 3 hold. Then for all $\delta\in(0,1)$ and all integer $n$ such that {ceqn}

[TABLE]

we have, with probability at least $1-\delta$,

[TABLE]

Compared to the bound (6) for the oracle estimator, the bound (8) for the OLS estimator has an additional term. This term is due to the additional learning step that is needed to estimate the optimal control variate.

Remark 9 (On the factor $B$)

Defined as the supremum of the leverage function $q_{n}$ in (Portier and Segers, 2019, Eq. (14)), the quantity $B$ plays an important role in our analysis as well as in other regression studies (Hsu et al., 2012; Newey, 1997). Just as the OLS estimate (see Remark 3), the quantity $B$ remains invariant by invertible linear transformation of the control variates. We have {ceqn}

[TABLE]

Remark 10 (On the parameters $\tau$ and $\gamma$)

The parameter $\tau$ in Assumption 1 is by definition an upper bound of the residual variance $\sigma^{2}_{m}(f)$. In many situations, its value is not too far from $\sigma^{2}_{m}(f)$. Hence, $\tau$ should capture the adequacy between the control variate space and the integrand $f$ and should decrease with $m$. The full rank condition expressed in Assumption 3 is not crucial as one could work with the Moore–Penrose inverse when solving (4). More importantly, a large value of the minimal eigenvalue $\gamma$ of the Gram matrix $G$ reflects that the OLS problem is well-conditioned, enhancing numerical stability. As control functions are added, rows and columns are added to $G$ and so $\gamma$ cannot increase. For the Fourier basis in Example 1, we have $\gamma=1$, while for the Legendre polynomials in Example 2, we have $\gamma\simeq 1/m$.

Remark 11 (Link with OLS prediction risk analysis)

The approach taken in the proof of Theorem 3.1 requires to bound what is called the prediction risk, defined as $\|G^{1/2}(\hat{\beta}_{n}^{\mathrm{ols}}(f)-\beta^{\star}(f))\|_{2}$. With probability greater than $1-\delta$, we obtain an upper bound of order $\sqrt{B\tau^{2}\log(m/\delta)/n}$ on the prediction risk. This makes our approach comparable to the one of the recent study (Hsu et al., 2012) where concentration bounds for the OLS prediction risk (and ridge) with random design are established. In contrast to their bound, our bound involves the quantity $B$ which shares the same invariant property as the OLS estimate and we don’t require the noise to be sub-Gaussian conditionally on the covariate but just sub-Gaussian which is weaker.

Remark 12 (Rates)

Consider an asymptotic set-up where the number of control variates $m$ tends to infinity with the Monte Carlo sample size $n$. The OLS method improves upon the basic MC method ($m=0$), which has rate $1/\sqrt{n}$, as soon as $\tau+\tau\sqrt{{mB}\log(m)/{n}}\to 0$. To recover the same order as the one of the oracle estimator $\hat{\alpha}_{n}^{\mathrm{or}}(f)$, which has rate $\tau/\sqrt{n}$, one must have ${mB}\log(m)=O(n)$ as $n\to\infty$, that is, $m$ must not be too large compared to $n$.

Remark 13 (Leverage condition)

Theorem 3.1 may be seen as a non-asymptotic version of the asymptotic results provided in Portier and Segers (2019) in which the leverage condition, $\sup\{h(x)^{T}G^{-1}h(x):x\in\mathcal{X}\}=o(n/m)$, is required to obtain a similar (asymptotic) bound (see Theorem 1 therein) as the one of Theorem 3.1. In the present non-asymptotic version, the leverage condition is expressed through $mB$ when requiring that $18B\log(4m/\delta)\leqslant n$.

LASSO takes advantage of sparse regression models. A regression model is sparse whenever many of the coefficients of the parameter vector $\beta$ are equal to zero, i.e., many of the covariates are useless to predict the output in the presence of the other covariates. The number of elements in the active set of the vector of regression coefficients $\beta^{\star}(f)$, {ceqn}

[TABLE]

is denoted by $\ell^{\star}:=\left\lvert{\textstyle S^{\star}}\right\rvert$ and quantifies the level of sparsity associated to the regression model. To avoid trivialities, we tacitly assume that $S^{\star}$ is non-empty, so $\ell^{\star}\geqslant 1$. The factor $\ell^{\star}$ represents the level of sparsity of $f$ with respect to the control functions and plays an important role in describing the benefits of the LASSO over the OLS. No assumption is made on $\ell^{\star}$, which could be any integer in $\{1,\ldots,m\}$.

We follow the approach presented in (Tibshirani et al., 2015, Section 11.4.1) (see also Bickel et al. (2009); van de Geer and Bühlmann (2009)), in which the analysis of the LASSO is carried out using a restricted eigenvalue condition. For a vector $\beta\in\mathbb{R}^{m}$ and for a non-empty set $S\subset\{1,\ldots,m\}$, write $\beta_{S}=(\beta_{k})_{k\in S}$, seen as a (column) vector in $\mathbb{R}^{\left\lvert{\textstyle S}\right\rvert}$. Define a collection of cones of interest. For $\alpha>0$ and $S\subset\{1,\ldots,m\}$, we set $\overline{S}=\{1,\ldots,m\}\setminus S$ and {ceqn}

[TABLE]

Assumption 4** (Restricted eigenvalue condition)**

There exists $\gamma^{\star}>0$ such that $u^{T}Gu\geqslant\gamma^{\star}\left\lVert{\textstyle{u}}\right\rVert_{2}^{2}$ for all $u\in\mathcal{C}(S^{\star};3)$.

In practice, we do not know the active set $S^{\star}$, so the only way to ensure Assumption 4 is to make sure all control variates $h_{1},\ldots,h_{m}$ are linearly independent. The practical value of the assumption is that $\gamma^{\star}\geqslant\gamma$, yielding sharper bounds below.

Recall that the $\ell_{1}$-penalty of the LASSO is weighted by a regularization parameter $\lambda>0$.

Theorem 3.2 (Concentration inequality for LASSO)

Suppose Assumptions 1, 2 and 4 hold. Introduce $\xi=\ell^{\star}(U_{h}^{2}/\gamma^{\star})$. Then for all $\delta\in(0,1)$ and all integer $n$ such that

[TABLE]

we have, with probability at least $1-\delta$,

[TABLE]

For $\lambda$ equal to the lower bound, we have on the same event

[TABLE]

Remark 14 (LASSO vs OLS)

The benefits of LASSO over OLS can be observed by comparing the bounds in (8) and (10). The total number $m$, of control functions has been replaced by the active number $\ell^{\star}$ of such functions. Further, because $\Gamma_{S^{\star}}=\{u\in\mathbb{R}^{p}\,:\,\|u\|_{2}=1,\,u\in\mathcal{C}(S^{\star};3)\}$ is included in the unit sphere, $\gamma^{\star}=\inf_{u\in\Gamma_{S^{\star}}}u^{T}Gu$ in Assumption 3 is at least as large as the smallest eigenvalue of $G$, $\gamma=\inf_{\|u\|_{2}=1}u^{T}Gu$ in Assumption 4.

The theoretical analysis of the LSLASSO estimator depends on the success of the LASSO-based model selection, i.e., the LASSO needs to correctly recover all the components of the true model. To ensure this selection step, the restricted eigenvalue condition is replaced by the two following ones.

Assumption 5** (Linear independence of active functions)**

The active control variates $h_{k}$, $k\in S^{\star}$, are linearly independent. As a consequence, the $\ell^{\star}\times\ell^{\star}$ Gram matrix $G_{S^{\star}}=P(h_{S^{\star}}h_{S^{\star}}^{T})$ is positive definite and its smallest eigenvalue $\gamma^{\star\star}:=\lambda_{\min}(G_{S^{\star}})$ is strictly positive.

Note that because $\{u\in\mathbb{R}^{p}\,:\,\|u\|_{2}=1,\,\forall k\notin S,\,u_{k}=0\}\subset\Gamma_{S}$ (introduced in remark 14), we have that $\gamma^{\star\star}\geqslant\gamma^{\star}$. Finally, it is required that that the active control functions are orthogonal, in $L_{2}(P)$, to the inactive ones.

Assumption 6** (Orthogonality)**

We have $P(h_{j}h_{k})=0$ for all $j\in\{1,\ldots,m\}\setminus S^{\star}$ and all $k\in S^{\star}$.

Since we do not know $S^{\star}$ in practice, the way to ensure Assumption 6 is by making all control variates orthogonal: $P(h_{j}h_{k})=0$ for all $j,k\in\{1,\ldots,m\}$. The Gram matrices $G$ and $G^{\star}$ are then diagonal. In the absence of zero control variates, Assumptions 3 and 4 are then satisfied as well, with $\gamma^{\star\star}=\min_{k\in S^{\star}}P(h_{k}^{2})\geqslant\min_{k=1,\ldots,m}P(h_{k}^{2})=\gamma>0$.

Theorem 3.3 (Support recovery of LASSO)

Suppose Assumptions 1, 2, 5 and 6 hold. Then for all $\delta\in(0,1)$, all integer $n$ such that {ceqn}

[TABLE]

and all $\lambda$ such that {ceqn}

[TABLE]

it holds that, with probability at least $1-\delta$, the LASSO based solution $\hat{\beta}_{n}^{\mathrm{lasso}}(f)$ is unique and the true active set is recovered, $\operatorname{supp}(\hat{\beta}_{n}^{\mathrm{lasso}}(f))=S^{\star}$.

The upper and lower bounds on $\lambda$ in (11) must not contradict each other, and this effectively implies an additional lower bound on $n$. Define $B^{\star}=\sup_{x\in\mathcal{X}}h_{S^{\star}}^{T}(x)G_{S^{\star}}^{-1}h_{S^{\star}}(x)$ and note that {ceqn}

[TABLE]

Theorem 3.4 (Concentration inequality for LSLASSO)

Suppose Assumptions 1, 2, 5 and 6 hold. Write $\xi^{\star}=\ell^{\star}(U_{h}^{2}/\gamma^{\star\star})$. Then for all $\delta\in(0,1)$ and all integer $N\in\{1,\ldots,n\}$ such that

{ceqn}

[TABLE]

and all $\lambda$ such that {ceqn}

[TABLE]

we have, with probability at least $1-\delta$,

[TABLE]

The logic behind Theorem 3.4 is that, by Theorem 3.3, the active set $\hat{S}_{N}=\operatorname{supp}(\hat{\beta}_{N}^{\mathrm{lasso}}(f))$ identified by means of the subsample of size $N$ is equal to the true active set $S^{\star}=\operatorname{supp}(\beta^{\star}(f))$ with large probability. On the event that the two sets coincide, the LSLASSO estimator is then the same as the OLS estimator based on the active control variates only, and the error bound follows from Theorem 3.1. In practice, it turns out that LSLASSO works well even when the true active set is not identified perfectly. However, to show this formally remains an open problem.

The assumptions and concentration inequalities in our theorems feature explicit rather than generic constants. Although we have worked hard to keep these constants under control [see in particular the proof of Lemma 4 as well as Step 6(ii) in the proof of Theorem 3.3], it is likely that, at the cost of lengthier computations, sharper constants can still be found.

Remark 15 (Bounded control variates)

In Assumption 2, the control variates were assumed to be bounded. Even if this assumption is valid for the two classic families in Examples 1 and 2 below, it might fail when control variates are produced with the Stein’s method as suggested in Remark 1. The boundedness assumption is needed to keep the same variance factor $\tau^{2}$ in the sub-Gaussian property of both variables $\epsilon(X_{1})$ and $\epsilon(X_{1})h(X_{1})$; see, e.g., Step 3.2 in the proof of Theorem 3.1 or Equation (35) in the proof of Theorem 3.2. Avoiding this assumption is thus possible at the price of more specific assumptions on the sub-Gaussianity of $\epsilon(X_{1})h(X_{1})$. Note finally that (different) asymptotic results are valid for unbounded control variates Portier and Segers (2019).

Remark 16 (Overfitting)

Theorems 3.2 and 3.4 advocate the use of the LASSO in favor of the OLS in scenarios where $\ell^{*}$ is smaller than $m$ or in the presence of collinearities in the design matrix making the parameter $\gamma$ close to zero; see also Remark 14. Another notable advantage of the (LS)LASSO and more generally of penalization methods, is the ability to prevent over-fitting. This occurs when the number of control variates $m$ is large compared to the Monte Carlo sample size $n$ or, more generally, when the approximation space is large compared to the sample size. While the theory developed here is unable to address such phenomena, one of the objectives of the numerical experiments conducted in the next section is to empirically demonstrate the superior performance of the LASSO-based methods even in the absence of sparsity.

To illustrate the application of our results in a standard framework, we consider two classic families of control functions, the Fourier basis and the Legendre polynomials.

Example 1 (Fourier basis)

On $\mathcal{X}=[0,1]$ equipped with the uniform distribution $P$, let $h_{j}(x)$ be equal to $\sqrt{2}\cos((j+1)\pi x)$ is $j$ is odd and to $\sqrt{2}\sin(j\pi x)$ is $j$ is even. The Fourier basis is orthonormal so that the Gram matrix is the identity, $G=I_{m}$, and $\gamma=\gamma^{\star}=\gamma^{\star\star}=1$. The cosine and sine functions being bounded by 1, a uniform bound is $U_{h}=\sqrt{2}$, which implies $B\leqslant 2m,B^{\star}\leqslant 2\ell^{\star}$. Under the proper assumptions, we get from Theorems 3.1 and 3.4 that with probability at least $1-\delta$, since $58\sqrt{2}<83$,

[TABLE]

and

[TABLE]

Example 2 (Legendre polynomials)

Suppose that $h_{j}=L_{j}$ is the Legendre polynomial of degree $j\in\{1,\ldots,m\}$. The Legendre polynomials are orthogonal on $\mathcal{X}=[-1,1]$ with respect to the uniform distribution $P$ and satisfy $|L_{j}(x)|\leqslant 1$ for $x\in[-1,1]$ with $L_{j}(1)=1$ and {ceqn}

[TABLE]

The Gram matrix $G=P(hh^{T})$ is diagonal with entries $1/(2j+1)$, so the minimum eigenvalue is $\gamma=1/(2m+1)$ and a uniform bound is $U_{h}=1$. Consequently, $B\leqslant 2m+1$. Similarly, considering only active control variates, we have $U_{h}^{\star}=1$, while the smallest eigenvalue, $\gamma^{\star\star}$, of $G_{S^{\star}}$ satisfies $1/(2m+1)\leqslant\gamma^{\star\star}\leqslant 1/(2\ell^{\star}+1)$. Under suitable assumptions, we get from Theorems 3.1 and 3.4 that with probability at least $1-\delta$,

[TABLE]

[TABLE]

Compared to the Fourier basis, the improvement of LSLASSO over the OLS estimator is not only related to the number of active varables $\ell^{*}$ compared to $m$ but also to the place of the active variables within the set of Legendre polynomials.

4 Numerical illustration

To compare the finite-sample performance of the various control variate methods, we consider synthetic data examples involving the standard integration problem over the unit cube $[0,1]^{d}$. The goal is to compute $\int_{[0,1]^{d}}f(x)\,\mathrm{d}x$. We shall consider various dimensions $d\geqslant 1$, different integrands $f:\mathbb{[}0,1]^{d}\to\mathbb{R}$, and several choices for the Monte Carlo sample size, $n$, and the number of control variates, $m$. We shall focus on difficult situations where $d$ is relatively large compared to $n$. In Section 5, we turn to real data examples in the context of Bayesian inference. For the sake of reproducibility, the data and Python code are available online111https://github.com/RemiLELUC/ControlVariateSelection.git.

Methods in competition.

We consider all the methods presented in Section 2 with two different strategies regarding the sub-sample size used to compute the active set in LSLASSO. The methods in competition are OLS, LASSO, LSLASSO (sub-sample size $N=n$) and LSLASSOX (sub-sample size $N=\lfloor 15\sqrt{n}\rfloor$). The latter choice accelerates the computation in a substantial manner without deteriorating too much the support recovery property of the LASSO. For synthetic data, because the integration domain is the unit cube $[0,1]^{d}$, Quasi-Monte Carlo (QMC) methods (Caflisch, 1998) are suitable for comparison. We run such methods in the experiments with two classical low-discrepancy sets of particles, namely Halton and Sobol sequences.

On the choice of $\lambda$.

In the LASSO-step of LSLASSO(X), the choice of the regularization parameter $\lambda$ is essential since it controls the number of active variables. It is common to tune this parameter using $K$-fold cross-validation at the price of additional computations. This method, presented in general form in Algorithm 1, uses the prediction error of the underlying regression problem as a proxy to calibrate the control variates estimate. In Algorithm 1, the “data” $X$ correspond to the matrix $H$ of observed control variables and the “labels” $y$ to the vector $f^{(n)}$ of observed function values. The method is computationally expensive, partitioning the training set in several folds and solving many regression problems for every value of $\lambda$ in a given grid.

To accelerate the computations, we suggest a new method based on a dichotomic search. Motivated by Eq. (8) and Remark 12, the value of $\lambda$ is tuned such that the number of selected control variates is of the order $\sqrt{n}$, which is the order obtained for $m$ when equating the two terms in (8) with $B=m$. Specifically, we enforce the number of activated control functions to lie in the range $[c_{1}\sqrt{n},c_{2}\sqrt{n}]$ for constants $0<c_{1}<c_{2}$ to be chosen (see below). This choice offers two advantages. On the one hand, the upper bound $c_{2}\sqrt{n}$ ensures that the number of selected control variates is relatively small compared to the sample size $n$, promoting stability and fast computation in the final OLS step. On the other hand, the lower bound $c_{1}\sqrt{n}$ reduces the risk of excluding relevant control variates.

The full procedure for the selection of the regularization parameter using a dichotomic search is described below in Algorithm 2. In all experiments, we set $c_{1}=3$ and $c_{2}=12$. We initialize $\lambda=\lambda_{\infty}$ to be the smallest value of $\lambda$ for which $\hat{\beta}^{\mathrm{lasso}}=0$, that is, $\lambda_{\infty}=\max_{k=1,\ldots,m}|H_{c,k}^{(N)T}f_{c}^{(N)}|/N$, where $H_{c,k}^{(N)}$ stands for the $k$-th column of $H_{c}^{(N)}$, which is the same as the matrix $H_{c}$ but then based on the first $N$ Monte Carlo draws (Tibshirani et al., 2015, Exercise 2.1). Next, we decrease the value of $\lambda$, e.g., by dividing it by two, such as to incorporate more and more control variates. If too many control functions are selected, i.e., more than $c_{2}\sqrt{n}$, we increase the value of $\lambda$ again, e.g., by multiplying it by two, to finally reach the desired range for the number of active variables. In the end, this procedure ensures a straightforward computation of the LSLASSO(X) because the size of the associated linear system remains reasonable. Contrary to $K$-fold cross-validation, it is not necessary to split the data into multiple folds, leading to a reduced computation time.

The pseudo-code of the corresponding LSLASSO(X) method is provided in Algorithm 3. The regression coefficients $\hat{\beta}_{n}^{\mathrm{ols}}$ and $\hat{\beta}_{n}^{\mathrm{lasso}}$ for OLS and LASSO are computed using the Scikit-Learn library (Pedregosa et al., 2011), employing coordinate descent to solve the LASSO problem.

Integrands.

We consider several integrands $f$ on $[0,1]^{d}$: {ceqn}

[TABLE]

and for all $j=1,\ldots,d$, {ceqn}

[TABLE]

All these functions integrate to $1$ on $[0,1]^{d}$. The functions $f_{j}$ and $g_{j}$ are built using tensor products of log-normal and exponential density functions, respectively, and depend on the first $j$ coordinates only. This construction ensures that for small $j$, the integrands $f_{j}$ and $g_{j}$ lend themselves to Monte Carlo integration based on selected control variates. In contrast, the functions $\varphi$, $f_{d}$ and $g_{d}$ represent more difficult situations where all the coordinates are involved and the symmetry of their role makes it harder to select some meaningful control functions. None of the integrands belongs to the linear span of the control variates constructed in the next paragraph.

Control variates.

Multidimensional control functions with respect to the uniform distribution over $[0,1]^{d}$ are easy to construct based on univariate ones. Let $(h_{1},\ldots,h_{k})$ be a vector of one-dimensional control functions, i.e., $\int_{0}^{1}h_{j}(x)\,\mathrm{d}x=0$ for each $j=1,\ldots,k$. Let $h_{0}=1$ denote the constant function equal to one. Without further information on the integrand, the usual way to construct multivariate controls is by forming tensor products of the form

[TABLE]

for a multi-index $\ell=(\ell_{1},\ldots,\ell_{d})$ in $\{0,\ldots,k\}^{d}\setminus\{(0,\ldots,0)\}$, yielding a total number of $(k+1)^{d}-1$ control functions.

A drawback of such a construction is that the number of control functions grows quickly with $k$. Alternative approaches yielding smaller control spaces consist of imposing $\ell_{j}=0$ for all but a small number (one or two, say) of coordinates $j=1,\ldots,d$ or simply picking at random a desired number, say $m$, of indices $\ell=(\ell_{1},\ldots,\ell_{d})$.

In this study, the set of control variates at our disposal is constructed as follows. We consider different settings of dimension $d$ with $k$ univariate control functions in each dimension. For $j\in\{1,\ldots,k\}$, let $h_{j}(x)=L_{j}(2x-1)$ for $x\in[0,1]$, with $L_{j}$ the univariate Legendre polynomial (Legendre function of the first kind) of degree $j$; see Example 2. We have $\int_{0}^{1}h_{j}(x)\,\mathrm{d}x=0$ for all $j=1,\ldots,m$. Because the Legendre polynomials are orthogonal, they provide some numerical stability when inverting the Gram matrix. The multivariate control functions are sorted in ascending order according to the total degree $\sum_{j=1}^{d}\ell_{j}\in\{1,\ldots,,kd\}$ of the polynomial. In the experiments, the number of control functions $m$ is increased by progressively including all polynomials whose total degree is lower than or equal to a fixed threshold $\mathit{deg}$.

Settings.

For the triple $(d,k,n)$ we consider $d\in\{3,5,8\}$, $k\in\{12,10,3\}$ and $n\in\{2\,000,\;5\,000,\;10\,000\}$. For each choice of $(d,k)$, the number of control variates $m$ with a total degree lower than or equal to a fixed threshold $\mathit{deg}$ are given in Table 1. The case $d=8$ represents a difficult situation as the number of points $n$ is relatively small compared to the dimension. For instance, a grid made of only four points in each direction would already comprise $65\,536$ points.

Results.

The different Monte Carlo estimates are compared on the basis of their mean squared error (MSE). Figure 1 presents the boxplots obtained over $100$ replications of the values returned by each of the methods. In Tables 6 to 6, we provide the ratio $MSE(\mathrm{vanilla})/MSE(\cdot)$, the MSE of the vanilla Monte Carlo estimate divided by the MSE for the current method, as a measure of statistical efficiency of the method relative to naive Monte Carlo integration. The four tables correspond to the four panels (a) to (d) in Figure 1. For a given number of control variates $m$, the most efficient method is indicated in bold. For the Lasso-based methods, the results for the $\lambda$ selection based on cross-validation (Algorithm 1) and dichotomic search (Algorithm 2) did not differ much; for the sake of brevity, the figures and the tables report the results associated to the dichotomic search.

Figures 1(a) and 1(b) highlight the success or failure of the OLS estimator depending on the size of $m$ compared to $n$. In Figures 1(c) and 1(d), we consider larger values of $m$ and only compare the Lasso-based methods as it takes too much time to solve the OLS. In all our experiments, the LSLASSOX is the clear winner as it has the highest accuracy in almost all configurations. Moreover, the LSLASSOX can be computed much faster than the LSLASSO: in our implementation, preselecting the control variates based on a smaller subsample led to a reduction op the computation time by a factor between three and twenty.

In Figure 1(a), boxplots of the values returned by each of the methods are provided for $\varphi$ in (14) when $d=3$ and $n=10\,000$. In this situation, where $m$ is small compared to $n$, the OLS performs very well and the LSLASSO procedure selects almost all control variates so it performs as well as OLS. In Figure 1(b), boxplots of the values returned by each of the methods are provided for $g_{3}$ in (16) when $d=5$, $n=2\,000$, and $N=700$. In this case, the OLS estimator starts to break down as soon as the number, $m$, of control variates is of the same order as $n$. It is then necessary to perform some control variate selection, which is succesfully carried out by the LASSO and LSLASSO. Both of these estimators give the best results. Although the number of sample points used in the selection step of LSLASSOX has been reduced compared to the LSLASSO, the stability of the active set is barely affected. Accordingly, the error distributions for LSLASSO and LSLASSOX are quite similar.

Figures 1(c) and 1(d) reveal the benefits of selecting appropriate control variates before applying the OLS estimator. Figure 1(c) covers the function $f_{1}$ in (15) when $d=5$ and $n=5\,000$, while Figure 1(d) deals with the function $g_{4}$ in (16) when $d=8$ and $n=2\,000$. In the latter case, the number of control variates, $m=36\,813$, is huge compared to the sample size $n=2\,000$. However, the Lasso-based methods perform remarkably well in those settings. More precisely, in dimension $d=5$ with the function $f_{1}$, the mean square error of the naive Monte Carlo estimator is of the order $10^{-5}$ whereas the one of the LSLASSOX is of the order $10^{-10}$. Similarly, in dimension $d=8$ with the function $g_{4}$, the mean square error goes down from $10^{-4}$ to $10^{-8}$. Table 6 highlights the benefits of the LSLASSO over the LASSO in difficult situations.

In the recent study South et al. (2018), the authors investigate the use of regularization in computing control variates estimates. They focus on the LASSO and ridge regression and they show, based on several examples, that the LASSO generally outperforms the ridge. In the applications they consider, they found that polynomials with relatively small degrees in each direction ($k$ equal to $2$ and $3$) give the best performance. The examples considered here show a similar pattern as the results do not generally improve beyond degree $k=3$.

5 Bayesian inference

In this section, we compare the different Monte Carlo estimates on Bayesian inference examples. Given some observed data $x$, the goal is to infer the parameter $\theta$ of a statistical model. We have some information through the prior distribution $\pi(\theta)$ and observe the model likelihood $\ell(x|\theta)$. Bayes’ rule gives the posterior distribution as

[TABLE]

The normalizing constant in the denominator is called evidence and is of interest for Bayesian model selection:

[TABLE]

Typically, this integral is analytically intractable. It is also difficult to compute numerically if the dimension $d$ of the parameter space $\Theta$ is large.

We consider the same datasets as in (South et al., 2018): the European dipper capture-recapture data from (Marzolin, 1988) in Section 5.1 and the sonar data from (Gorman and Sejnowski, 1988) in Section 5.2. The dimensions of the integration domains are $d=12$ and $d=61$, respectively.

As in Section 4, we consider multivariate control functions based on univariate orthogonal polynomials by forming tensor products of the form $h_{\ell}(x_{1},\ldots,x_{d})=\prod_{j=1}^{d}h_{\ell_{j}}(x_{j})$, for a multi-index $\ell=(\ell_{1},\ldots,\ell_{d})$ in $\{0,\ldots,k\}^{d}\setminus\{(0,\ldots,0)\}$. In both examples, the dimension $d$ is so large that considering all tensor products is infeasible. Instead, we focus on combinations where $\ell_{j}$ equals [math] for all but one or two coordinates, leading to a total number of $m=kd$ and $m=kd+k^{2}d(d-1)/2$ control variates, respectively.

The different Monte Carlo estimates are compared on the basis of their mean squared errors (MSE). In contrast to Section 4, the true value of the integral is unknown. An estimate of this value, referred to as the gold standard $Z^{\star}$, is obtained by naive Monte Carlo with sample size $n=10^{8}$. The variance of this estimate, computed on $20$ independent runs, is smaller than the variance of all the other considered methods. The different boxplots of Figure 2 show the results obtained over $100$ independent runs of $\hat{Z}/Z^{\star}$ where $\hat{Z}$ is the estimate of the evidence. Tables 11 to 11 provide numerical values for the statistical efficiency $\widehat{MSE}(\mathrm{vanilla})/\widehat{MSE}(\cdot)$. We consider various settings and the parameter configuration is $n\in\{2\,000;\,5\,000\}$ for the Monte Carlo sample size with $N\in\{700;\,1\,000\}$ for the Monte Carlo subsample size for the LSLASSO. The regularization parameter $\lambda$ is chosen via dichotomic search (Algorithm 2).

5.1 European dipper capture-recapture data

The data-set given in Table 7 was collected by (Marzolin, 1988) and describes the annual capture and recapture counts of the bird species Cinclus cinclus, also known as the European dipper, in eastern France from 1981 to 1987. We observe count data $x_{i,j}$ with $i\in\{1,\ldots,I\}$ and $j\in\{i+1,\ldots,J\}$, where $x_{i,j}$ denotes the number of birds released in year $i$ and subsequently recaptured for the first time in year $j$. In the example, we have $I=6$ and $J=7$, where 1981 corresponds to year $i=1$. Also observed is $R_{i}$, the number of marked birds released into the population in year $i$.

Following (Brooks et al., 2000; Nott et al., 2018; South et al., 2018), we consider a Bayesian approach for the Cormack–Jolly–Seber model (Lebreton et al., 1992). The model parameters are $\phi_{i}$, a bird’s survival probability from year $i$ to $(i+1)$ for $i\in\{1,\ldots,I\}$, together with $p_{j}$, the probability of a bird being recaptured in year $j\in\{2,\ldots,J\}$. Let $\nu_{i,j}$ denote the probability that a bird captured and released in year $i$ gets recaptured for the first time in year $j$. Since the bird must survive from year $i$ to year $j$, not be recaptured in years $i+1$ to $j-1$ and then finally be recaptured in year $j$, the probability is modelled as

[TABLE]

The number of birds released at year $i$ that are never recaptured at all is equal to $r_{i}=R_{i}-\sum_{j=i+1}^{J}x_{i,j}$ while the probability that a bird released in year $i$ is never recaptured is $\chi_{i}=1-\sum_{j=i+1}^{J}\nu_{i,j}$. The resulting likelihood is equal to

[TABLE]

where $\theta=(\phi_{1},\ldots,\phi_{6},p_{2},\ldots,p_{7})\in[0,1]^{12}$. The uniform distribution is chosen as prior and we use tensor products of Legendre polynomials with $k=10$ (Example 2) as controls.

The results for the various integration methos are reported in the same way as in Section 4. The boxplots and statistical efficiencies are given in Figures 2(a) and 2(b) and Tables 11 and  11 respectively. Similarly to the synthetic data, Figures 2(a) and 2(b) reveal the success or failure of the OLS on the capture-recapture data when the number of control variates $m$ is larger than the Monte Carlo sample size $n$. The variance goes down as $m$ increases. Tables 11 and 11 show that for $n=2\,000$, the OLS estimate gives the best performance whereas for $n=5\,000$, the LASSO-based methods profit from the large number of available control variates. In this case, the LASSO is most efficient while the LSLASSOX performs similarly but at a reduced computing time.

5.2 Sonar data

The data were collected by (Gorman and Sejnowski, 1988) and are available from the UCI Machine Learning Repository (Asuncion and Newman, 2007). The data matrix $X$ represents $208$ sonar signals, each one composed of $60$ attributes within the binary classification framework. A column of $1$’s is added to the matrix $X$ to represent the intercept so that $X\in\mathbb{R}^{208\times 61}$. The goal is to assess whether the sonar signal bounces off a metal cylinder (label $y=1$) or a roughly cylindrical rock (label $y=-1$). The different covariates represent the energy within particular frequency bands, integrated over a certain period of time. Using the encoding $y\in\{-1,+1\}$ and following a logistic regression model, the resulting log-likelihood is

[TABLE]

where the model coefficient $\theta\in[-1,1]^{61}$ has a uniform prior distribution. We use the family of Legendre polynomials as control functions with $k=20$. The boxplots and statistical efficiencies are presented in Figures 2(c) and 2(d) and Tables 11 and 11, respectively. Once again, the Lasso-based methods, with their selection strategy, are able to benefit from a larger control variates space. The winner of this competition is LSLASSOX as it offers the best performance combined with a smaller computation time compared to the LSLASSO.

6 Conclusion and perspective

The use of high-dimensional control variates with the help of a LASSO-type procedure has been shown to be efficient in order to reduce the variance of the basic Monte Carlo estimate. The method, called LSLASSO(X), that first selects appropriate control variates by the LASSO, possibly on a smaller subsample, and then estimates the control variate coefficients by least squares performs excellently considering the modest computing time required. Several avenues for further research are now discussed.

The construction of control variates by a change of measure (Remark 1) presupposes some knowledge on the underlying integration measure in order to choose an appropriate sampling distribution. For instance, if the support of the sampling measure does not cover the whole integration domain then the method will certainly fail. Adaptive importance sampling (see, e.g., (Owen and Zhou, 2000; Portier and Delyon, 2018)) offers a possible solution, involving online estimates of the appropriate sampling policy and the optimal linear combination of control variates.

Assumption 1 on the sub-Gaussianity of the residuals is key to obtain concentration inequalities. For certain applications, it might by too restrictive, however. In the absence of such an assumption or more generally of suitable bounds on the tails of the residual distribution, other types of results such as almost sure convergence rates might still be pursued.

In the random design setting, the estimators of coefficient vector $\beta^{\star}(f)$ are all biased, even the OLS estimator. The bias may be removed by sample splitting (Avramidis and Wilson, 1993), but at the cost of an increased variance, especially if the number of control variates is large. For the Lasso-based methods, debiasing methods are studied in (Javanmard and Montanari, 2018) and the references therein. The merits of these techniques for Monte Carlo control variate methods remain to be investigated.

We have presented different control variate methods from the point of view of estimation only. Equally important questions are that of model evaluation and Monte Carlo sample size calculation, assessing the accuracy of the estimate. Several ways can be imagined such as sample splitting (e.g., cross-validation) and plug-in estimation of the residual variance $\sigma^{2}(f)$, using for instance the estimated residuals.

Acknowledgments.

The authors are grateful to the Associate Editor and two anonymous Reviewers for their many valuable comments and interesting suggestions. Further, the authors gratefully acknowledge the exceptionally careful proofreading done by Aigerim Zhuman (UCLouvain).

Appendix A Auxiliary results

Lemma 1

*(Sub-Gaussian)

Let $X_{1},\ldots,X_{n}$ be independent and identically distributed random variables in $(\mathcal{X},\mathcal{A})$ with distribution $P$. Let $\varphi_{1},\ldots,\varphi_{p}$ be real-valued functions on $\mathcal{X}$ such that $P(\varphi_{k})=0$ and $\varphi_{k}\in\mathcal{G}(\tau^{2})$ for all $k=1,\ldots,p$. Then for all $\delta>0$, we have with probability at least $1-\delta$, {ceqn}*

[TABLE]

Proof

For each $k=1,\ldots,p$, the centered random variable $\sum_{i=1}^{n}\varphi_{k}(X_{i})$ is sub-Gaussian with variance factor $n\tau^{2}$. By the union bound and by Chernoff’s inequality, we have, for each $t>0$, {ceqn}

[TABLE]

Set $t=\sqrt{2n\tau^{2}\log(2p/\delta)}$ to find the result.

Lemma 2

*(Smallest eigenvalue lower bound)

Let $X_{1},\ldots,X_{n}$ be independent and identically distributed random variables in $(\mathcal{X},\mathcal{A})$ with distribution $P$. Let $g=(g_{1},\ldots,g_{p})^{T}$ in $L_{2}(P)^{p}$ be such that the $p\times p$ Gram matrix $G=P(gg^{T})$ satisfies $\lambda_{\min}(G)>0$. Define the transformation $\tilde{g}=G^{-1/2}g$ and put $B_{\tilde{g}}:=\sup_{x\in\mathcal{X}}\left\lVert{\textstyle\tilde{g}(x)}\right\rVert_{2}^{2}$. Let $\delta,\eta\in(0,1)$. For $\delta\in(0,1)$, the empirical Gram matrix $\hat{G}_{n}=P_{n}(gg^{T})$ satisfies, with probability at least $1-\delta$,*

[TABLE]

Proof

Suppose that the result is true in the special case that $G$ is the identity matrix. In case of a general Gram matrix $G$, we could then apply the result for the special case to the vector of functions $\tilde{g}=G^{-1/2}g$, whose Gram matrix is the identity matrix. We would get that $\lambda_{\min}(P_{n}(\tilde{g}\tilde{g}^{T}))>1-\eta$ with probability at least $1-\delta$. Since $P_{n}(\tilde{g}\tilde{g}^{T})=G^{-1/2}\hat{G}_{n}G^{-1/2}$ and since $u^{T}G^{-1}u\leqslant 1/\lambda_{\min}(G)$ for every unit vector $u\in\mathbb{R}^{p}$, we would have {ceqn}

[TABLE]

It would then follow that

[TABLE]

as required. Hence we only need to show the result for $G=I$, in which case $\tilde{g}=g$.

We apply the matrix Chernoff inequality in (Tropp, 2015, Theorem 5.1.1) to the random matrices $n^{-1}g(X_{i})g(X_{i})^{T}$. These matrices are independent and symmetric with dimension $p\times p$. Their minimum and maximum eigenvalues are between [math] and $L=B_{g}/n$, with $B_{g}=\sup_{x\in\mathcal{X}}\lambda_{\max}(g(x)g(x)^{T})=\sup_{x\in\mathcal{X}}\left\lVert{\textstyle g(x)}\right\rVert^{2}_{2}$. Their sum is equal to $P_{n}(gg^{T})=\hat{G}_{n}$, whose expectation is $G=I$ by assumption. In the notation of the cited theorem, we have $\mu_{\min}=\lambda_{\min}(G)=1$, and thus, by Eq. (5.1.5) in that theorem, we have, for $\eta\in[0,1)$,

[TABLE]

The term in square brackets is bounded above by $\exp(-\eta^{2}/2)$. Indeed, we have, for $\eta\in[0,1)$,

[TABLE]

and {ceqn}

[TABLE]

It follows that

[TABLE]

Solving $p\exp\left(-\frac{\eta^{2}n}{2B_{g}}\right)=\delta$ in $\eta$, we find that, with probability at least $1-\delta$,

[TABLE]

Lemma 3 (Upper bound of moments)

Let $X$ be a random variable such that $\operatorname{\mathbb{E}}(|X|^{2p})\leqslant 2^{p+1}p!$ for every integer $p\geqslant 1$. Then {ceqn}

[TABLE]

in which it is implicitly understood that the series on the left-hand side converges.

Proof

Let $\lambda\in\mathbb{R}$. We split the series in terms with even and odd indices $k$, leading to

[TABLE]

We will bound the series on the odd indices in terms of the series on the even indices.

Since the geometric mean of two nonnegative numbers is bounded by their arithmetic mean, we have, for all $x\geqslant 0$ and all $a>0$, {ceqn}

[TABLE]

Applying the previous inequality to $x=\lambda X$ and scalars $a_{p}>0$ to be chosen later,

[TABLE]

Here, $\mathds{1}$ denotes an indicator function. We obtain

[TABLE]

Define $b_{p}=a_{p}(2p+1)$ and use the hypothesis on $\operatorname{\mathbb{E}}(|X|^{2p})$ to see that, for any constants $b_{p}>0$,

[TABLE]

The objective is to find a constant $c>0$ as small as possible and such that the right-hand side is bounded by $\exp(c\lambda^{2})=1+\sum_{p=1}^{\infty}c^{p}\lambda^{2p}/p!$. Comparing coefficients, this means that we need to determine scalars $b_{p}>0$ and $c>0$ in such a way that for all $p=1,2,\ldots$ {ceqn}

[TABLE]

or, equivalently, {ceqn}

[TABLE]

The case $p=1$ gives {ceqn}

[TABLE]

showing that, with this proof technique, we will always find $c>2$. Setting $b_{p}\equiv b>0$ for all integer $p\geqslant 1$ and $c=2+1/b$, inequality (18) is automatically satisfied, so it remains to find $b>0$ such that forall $p=2,3,\ldots$

[TABLE]

The left-hand side is decreasing in $p$ whereas the right-hand side is increasing in $p$. It is thus sufficient to have the inequality satisfied for $p=2$, i.e., {ceqn}

[TABLE]

Equating both sides leads to a nonlinear equation in $b$ that can be solved numerically, giving the root $b\approx 4.006156$. With $b=4$, inequality (19) is satisfied, as can be checked directly ($91/72\leqslant 81/64$). We conclude that $c=2+1/4=9/4$ is a valid choice.

Note that the series in (17) starts at $k=2$. If also $\operatorname{\mathbb{E}}(X)=0$, the left-hand side in (17) is an upper bound for $\operatorname{\mathbb{E}}(\exp(\lambda X))$, and we obtain the following corollary.

Corollary 1

Let $Z$ be a centered random variable such that

[TABLE]

Then $\log\operatorname{\mathbb{E}}(\exp(\lambda Z))\leqslant 9\lambda^{2}/4$ for all $\lambda\in\mathbb{R}$, i.e., $Z\in\mathcal{G}(9/2)$.

Lemma 4

Let $(X,Y)$ be a pair of uncorrelated random variables. If $X\in\mathcal{G}(\nu)$ and $|Y|\leqslant\kappa$ for some $\nu>0$ and $\kappa>0$, then $XY\in\mathcal{G}((9/2)\kappa^{2}\nu)$.

Proof

The random variable $X/\sqrt{\nu}$ is sub-Gaussian with variance factor $1$. As on page 25 in Boucheron et al. (2013), this implies that $\operatorname{\mathbb{P}}(|X/\sqrt{\nu}|>t)\leqslant 2\exp(-t^{2}/2)$ for all $t\geqslant 0$ and thus $\operatorname{\mathbb{E}}[|X/\sqrt{\nu}|^{2p}]\leqslant 2^{p+1}p!$ for all integer $p\geqslant 1$ (see (Boucheron et al., 2013, Theorem 2.1)).

Let $Z=XY/(\sqrt{\nu}\kappa)$. Since $X$ is centered and $X$ and $Y$ are uncorrelated, $XY$ is centred too, and therefore also $Z$. From the previous paragraph, we have $\operatorname{\mathbb{E}}(|Z|^{2p})\leqslant\operatorname{\mathbb{E}}(|X/\sqrt{\nu}|^{2p})\leqslant 2^{p+1}p!$ for all integer $p\geqslant 1$. Corollary 1 gives for all $\lambda\in\mathbb{R}$ that $\log\operatorname{\mathbb{E}}(\exp(\lambda Z))\leqslant 9\lambda^{2}/4$, from which {ceqn}

[TABLE]

Lemma 5 (Upper bound for norm-subGaussian random vector)

Let $X$ be a $d$-dimensional random vector with zero-mean and such that $\operatorname{\mathbb{P}}(\left\lVert{\textstyle X}\right\rVert_{2}\geqslant t)\leqslant 2\exp\left(-t^{2}/(2\sigma^{2})\right)$ for all $t\geqslant 0$. Then the random matrix $Y$ defined by {ceqn}

[TABLE]

satisfies $\operatorname{\mathbb{E}}(\exp(\theta Y))\preceq\exp(c\theta^{2}\sigma^{2})I$ for any $\theta\in\mathbb{R}$, with $c=9/4$, where $I$ denotes the identity matrix.

Proof

The non-zero eigenvalues of $Y$ are $\left\lVert{\textstyle X}\right\rVert$ and $-\left\lVert{\textstyle X}\right\rVert$. The non-zero eigenvalues of $Y^{k}$ are thus $\left\lVert{\textstyle X}\right\rVert^{k}$ and $(-\left\lVert{\textstyle X}\right\rVert)^{k}$ for integer $k\geqslant 1$. It follows that $Y^{k}\preceq\left\lVert{\textstyle X}\right\rVert^{k}I$ for all integer $k\geqslant 1$, and therefore also $\operatorname{\mathbb{E}}(Y^{k})\preceq\operatorname{\mathbb{E}}(\left\lVert{\textstyle X}\right\rVert^{k})I$ for all integer $k\geqslant 1$. Furthermore, the operator norm of $Y^{k}$ is bounded by $\left\lVert{\textstyle Y^{k}}\right\rVert\leqslant\left\lVert{\textstyle X}\right\rVert^{k}$.

Since $\operatorname{\mathbb{E}}(Y)=0$, we get, for any $\theta\in\mathbb{R}$,

[TABLE]

where $\xi=\left\lVert{\textstyle X}\right\rVert/\sigma$. The first series converges in operator norm since $\left\lVert{\textstyle\operatorname{\mathbb{E}}(Y^{k})}\right\rVert\leqslant\operatorname{\mathbb{E}}(\left\lVert{\textstyle Y^{k}}\right\rVert)\leqslant\operatorname{\mathbb{E}}(\left\lVert{\textstyle X}\right\rVert^{k})$.

By assumption, $\operatorname{\mathbb{P}}(\xi>t)=\operatorname{\mathbb{P}}(\left\lVert{\textstyle X}\right\rVert\geqslant\sigma t)\leqslant 2e^{-t^{2}/2}$ for all $t\geqslant 0$ and thus $\operatorname{\mathbb{E}}(|\xi|^{2p})\leqslant 2^{p+1}p!$ for all integer $p\geqslant 1$ But then we can apply Lemma 3 with $\lambda=\theta\sigma$ and $X=\xi$, completing the proof.

The following result is a special case of Jin et al. (2019, Corollary 7). Our contribution is to make the constant $c$ in the cited result explicit. In passing, we correct an inaccuracy in the proof of Jin et al. (2019, Lemma 4), in which it was incorrectly claimed that the odd moments of a certain random matrix Y as in our Lemma 5 are all zero.

Lemma 6 (Hoeffding inequality for norm-subGaussian random vectors)

Let the $d$-dimensional random vectors $Z_{1},\ldots,Z_{n}$ be independent, have mean zero, and satisfy {ceqn}

[TABLE]

for some $\sigma>0$. Then for any $\delta>0$, with probability at least $1-\delta$, we have

[TABLE]

Proof

Given Corollary 7 in Jin et al. (2019), the only thing to prove is that their constant can be set equal to $3$. Their Corollary 7 follows from their Lemma 6 in which it is shown that when the matrix $Y$ defined in (20) satisfies

[TABLE]

then we have for any $\theta>0$, with probability at least $(1-\delta)$, {ceqn}

[TABLE]

Taking $\theta=\sqrt{\log(2d/\delta)/(cn\sigma^{2})}$ yields {ceqn}

[TABLE]

and we conclude with Lemma 5 ($c=9/4,2\sqrt{c}=3$).

Appendix B Proof of Theorem 3.1

The proof is organized as follows. We first provide an upper bound on the error (Step 1). This bound involves the norm of the error made on the rescaled coefficients and is controlled in Step 2. Then (Step 3), we construct an event that has probability at least $1-\delta$ on which we can control the terms that appear in the upper bound of Step 2. Collecting all the inequalities, we will arrive at the stated bound (Step 4).

Step 1. — Since $f=P(f)+\beta^{\star}(f)^{T}h+\epsilon$, the oracle estimate of $P(f)$, which uses the unknown, optimal coefficient vector $\beta^{\star}(f)$, is

[TABLE]

The difference between the OLS and oracle estimates is

[TABLE]

Let $G=P(hh^{\top})$ be the $m\times m$ Gram matrix. By assumption, $G$ is positive definite. Write {ceqn}

[TABLE]

The estimation error of the OLS estimator can thus be decomposed as {ceqn}

[TABLE]

By the triangle and Cauchy–Schwarz inequalities, {ceqn}

[TABLE]

Step 2. — We will show that, if $\lambda_{\min}(P_{n}(\hbar\hbar^{\top}))>\left\lVert{\textstyle P_{n}(\hbar)}\right\rVert_{2}^{2}$, then {ceqn}

[TABLE]

and thus, by (22),

[TABLE]

Step 2.1 — Considered the column-centered $n\times m$ design matrices {ceqn}

[TABLE]

Since $\bar{H}^{\top}\mathds{1}_{n}=nP_{n}(\hbar)$, we have {ceqn}

[TABLE]

As a consequence, for $u\in\mathbb{R}^{m}$, {ceqn}

[TABLE]

by the Cauchy–Schwarz inequality. In particular, $u^{\top}\bar{H}_{c}^{\top}\bar{H}_{c}u$ is non-zero for non-zero $u\in\mathbb{R}^{m}$, so that $\bar{H}_{c}^{\top}\bar{H}_{c}$ is invertible, and so is the matrix

[TABLE]

Also, the smallest eigenvalue of $\bar{H}_{c}^{\top}\bar{H}_{c}$ is bounded from below by

[TABLE]

The largest eigenvalue of the inverse matrix $(\bar{H}_{c}^{\top}\bar{H}_{c})^{-1}$ is then bounded from above by {ceqn}

[TABLE]

Step 2.2. — Write $\epsilon_{c}^{(n)}=(\epsilon(X_{i})-P_{n}(\epsilon))_{i=1}^{n}$ for the centered vector of error terms. Recall $f_{c}^{(n)}=(f(X_{i})-P_{n}(f))_{i=1}^{n}$, the centered vector of samples from the integrand. As $f=P(f)+h^{\top}\beta^{*}(f)+\epsilon$, we have

[TABLE]

From the characterization (4) of the OLS estimate of the coefficient vector and since $H_{c}^{\top}H_{c}$ is invertible, {ceqn}

[TABLE]

We obtain {ceqn}

[TABLE]

Step 2.3. — We combine the results from Steps 2.1 and 2.2. From the upper bound (25) and the identity (26), we obtain {ceqn}

[TABLE]

Finally, as $\bar{H}_{c}=(\hbar_{j}(X_{i})-P_{n}(\hbar_{j}))_{i,j}$, we find {ceqn}

[TABLE]

Equation (23) follows.

Step 3. — In view of (24), we need to ensure that $\left\lvert{\textstyle P_{n}(\epsilon)}\right\rvert$, $\left\lVert{\textstyle P_{n}(\hbar)}\right\rVert$ and $\left\lVert{\textstyle P_{n}(\hbar\epsilon)}\right\rVert_{2}$ are small and that $\lambda_{\min}(P_{n}(\hbar\hbar^{\top}))$ is large. Let $\delta>0$. We construct an event with probability at least $1-\delta$ on which four inequalities hold simultaneously. Recall $B=\sup_{x\in\mathcal{X}}\left\lVert{\textstyle\hbar(x)}\right\rVert_{2}^{2}$, defined in (7).

Step 3.1. — Because $\epsilon\in\mathcal{G}(\tau^{2})$, Chernoff’s inequality (or Lemma 1 with $p=1$) implies that with probability at least $1-\delta/4$, {ceqn}

[TABLE]

Step 3.2. — For the term $\left\lVert{\textstyle\sum_{i=1}^{n}\hbar(X_{i})}\right\rVert_{2}$, we apply the vector Bernstein bound in (Hsu et al., 2014, Lemma 9). On the one hand $\sup_{x\in\mathcal{X}}\left\lVert{\textstyle\hbar(x)}\right\rVert_{2}\leqslant\sqrt{B}$ and on the other hand

[TABLE]

The cited vector Bernstein bound gives {ceqn}

[TABLE]

Setting $t=\log(4/\delta)$, we find that, with probability at least $1-\delta/4$, we have

[TABLE]

Since $\log(4/\delta)\geqslant\log(4)$, we have {ceqn}

[TABLE]

and thus {ceqn}

[TABLE]

The condition on $n$ easily implies that

[TABLE]

and thus {ceqn}

[TABLE]

Step 3.3. — To control $\left\lVert{\textstyle\sum_{i=1}^{n}\hbar(X_{i})\epsilon(X_{i})}\right\rVert_{2}$, we apply Lemma 6 with $Z_{i}=\hbar(X_{i})\epsilon(X_{i})$. The random vectors $\hbar(X_{i})\epsilon(X_{i})$ for $i=1,\ldots,n$ are independent and identically distributed and have mean zero. Since $\left\lVert{\textstyle\hbar(X_{i})}\right\rVert_{2}\leqslant\sqrt{B}$ by (7) and since $\epsilon\in\mathcal{G}(\tau^{2})$ by Assumption 1, we have, for all $t>0$, {ceqn}

[TABLE]

and (21) holds with $\sigma^{2}=B\tau^{2}$. Lemma 6 then implies that, with probability at least $1-\delta/4$ and $c=3$ that {ceqn}

[TABLE]

Step 3.4. — Recall the $n\times m$ matrix $H=(h_{j}(X_{i}))_{i,j}$ and put

[TABLE]

The empirical Gram matrix of the vector $\hbar=(\hbar_{1},\ldots,\hbar_{m})^{\top}\in L_{2}(P)^{m}$ based on the sample $X_{1},\ldots,X_{n}$ is

[TABLE]

We apply Lemma 2 with $g=\tilde{g}=\hbar$, $p=m$, and $\delta$ replaced by $\delta/4$. We find that, with probability at least $1-\delta/4$, {ceqn}

[TABLE]

Since $P_{n}(\hbar\hbar^{\top})=n^{-1}\bar{H}^{\top}\bar{H}$, it follows that {ceqn}

[TABLE]

as the assumption on $n$ implies that $2Bn^{-1}\log(4m/\delta)\leqslant 1/9$.

By the union bound, the inequalities (27), (28), (29), and (30) hold simultaneously on an event with probability at least $1-\delta$. For the remainder of the proof, we work on this event, denoted by $E$.

Step 4. — We combine the bound (24) on the estimation error with the bounds valid on the event $E$ constructed in Step 3. By (31), we have

[TABLE]

since the assumption on $n$ implies that $25mn^{-1}\log(4/\delta)\leqslant 1/3$. As $B\geqslant m\geqslant 1$, we have

{ceqn}

[TABLE]

since, by assumption, $n\geqslant 75m\log(4/\delta)$ which implies that $\sqrt{n^{-1}\log(4/\delta)}\leqslant 1/(5\sqrt{3})$. We find {ceqn}

[TABLE]

and the value $c=3$ gives $15(c+\sqrt{2/3})\approx 57.2<58$ which is the bound stated in Theorem 3.1. ∎

Appendix C Proof of Theorem 3.2

For a vector $\beta\in\mathbb{R}^{m}$ and for a non-empty set $S\subset\{1,\ldots,m\}$, write $\beta_{S}=(\beta_{k})_{k\in S}$. For any matrix $A\in\mathbb{R}^{n\times m}$ and $k\in\{1,\ldots,m\}$, let $A_{k}$ denote its $k$-th column and if $S=\{k_{1},\ldots,k_{\ell}\}\subset\{1,\ldots,m\}$ with $k_{1}<\ldots<k_{\ell}$, write $A_{S}=(A_{k_{1}},\ldots,A_{k_{\ell}})\in\mathbb{R}^{n\times\ell}$.

The proof is organized in a similar way as the one of Theorem 3.1. We first provide an initial upper bound on the error (Step 1). Then we construct an event that (Step 2) has probability at least $1-\delta$ and (Steps 3, 4, 5) on which we can control each of the terms of the previous upper bound. The combination of all steps to deduce the final statement is made clear in Step 6.

Step 1. — As in the proof of Theorem 3.1, with $\hat{\beta}_{n}^{\mathrm{ols}}(f)$ replaced by $\hat{\beta}_{n}^{\mathrm{lasso}}(f)$, the estimation error of the LASSO estimator can be decomposed as

[TABLE]

Writing $\hat{u}=\hat{\beta}_{n}^{\mathrm{lasso}}(f)-\beta^{\star}(f)$, we get, by the triangle and Hölder inequalities, {ceqn}

[TABLE]

Step 2. — Let $\delta>0$. We construct an event, $E$, with probability at least $1-\delta$ on which four inequalities, namely (33), (34), (35) and (36), hold simultaneously.

• •

Since $\epsilon\in\mathcal{G}(\tau^{2})$, we can apply Lemma 1 with $p=1$ to get that, with probability at least $1-\delta/4$, {ceqn}

[TABLE]

• •

In view of (Boucheron et al., 2013, Lemma 2.2) and Assumption 2, we have $h_{k}\in\mathcal{G}(U_{h}^{2})$ for all $k=1,\ldots,m$. Hence we can apply Lemma 1 with $p=m$ to get that, with probability at least $1-\delta/4$, {ceqn}

[TABLE]

• •

By virtue of Assumptions 1 and 2, we can apply Lemma 4 to find $h_{k}\epsilon\in\mathcal{G}(C\tau^{2}U_{h}^{2})$ with $C=9/2$. Hence we can apply Lemma 1 to get that, with probability at least $1-\delta/4$, {ceqn}

[TABLE]

• •

In view of (Boucheron et al., 2013, Lemma 2.2) and Assumptions 2 and 6, we have $h_{k}h_{l}-P(h_{k}h_{l})\in\mathcal{G}(U_{h}^{4})$ for all $k,l\in\{1,\ldots,m\}$. Hence we can apply Lemma 1 with $p=m^{2}$ to get that, with probability at least $1-\delta/4$, {ceqn}

[TABLE]

Denote by $\Delta=(P_{n}-P)\{hh^{T}\}$. Because by assumption $2(\ell^{\star}/\gamma^{\star})\sqrt{2U_{h}^{4}\log(8m^{2}/\delta)}\leqslant\sqrt{n}$, we have that {ceqn}

[TABLE]

Remark that {ceqn}

[TABLE]

Then, following (Bickel et al., 2009, equation (3.3)), use the inequality $|u^{T}\Delta u|\leqslant\left\lVert{\textstyle u}\right\rVert_{1}^{2}\max_{1\leqslant k,l\leqslant m}|\Delta_{k,l}|$, to obtain that, with probability $1-\delta/4$, for all $u\in\mathcal{C}(S^{\star};3)$, {ceqn}

[TABLE]

It follows that with probability at least $1-\delta/4$, {ceqn}

[TABLE]

Step 3. — We claim that, on the event $E$, we have {ceqn}

[TABLE]

We have

[TABLE]

and thus, {ceqn}

[TABLE]

We treat both terms on the right-hand side. On the one hand, we just have obtained a lower bound for the first term. On the other hand, in view of (34) and because $\left\lVert{\textstyle u}\right\rVert_{1}^{2}\leqslant 16\left\lVert{\textstyle u_{S^{\star}}}\right\rVert_{1}^{2}\leqslant 16\ell^{\star}\left\lVert{\textstyle u}\right\rVert_{2}^{2}$, we have {ceqn}

[TABLE]

as $n\geqslant(16\times 8)\ell^{\star}(U_{h}^{2}/\gamma^{\star})\log(8m/\delta)$ by assumption. In combination with (36), we find {ceqn}

[TABLE]

Step 4. — We claim that, on the event $E$, we have

{ceqn}

[TABLE]

Indeed, on the left-hand side in (38) we have in virtue of (33), (34) and (35), {ceqn}

[TABLE]

Since $\ell^{\star}\geqslant 1$ and $\ell^{\star}U_{h}^{2}\geqslant\sum_{k\in S^{\star}}P(h_{k}^{2})\geqslant\gamma^{\star}$, the assumed lower bound on $n$ implies that $n\geqslant 128\log(8/\delta)$. As $C=9/2$, the factor $\sqrt{2C}(1+\sqrt{2\log(8/\delta)/(Cn)})$ is bounded by $3+\sqrt{2}/8$ and we get (38).

Step 5. — Recall $\hat{u}=\hat{\beta}_{n}^{\mathrm{lasso}}(f)-\beta^{\star}(f)$. We claim that, on the event $E$, we have {ceqn}

[TABLE]

To prove this result, we shall rely on the following lemma.

Lemma 7

If $n\lambda\geqslant 2\left\lVert{\textstyle H_{c}^{T}\epsilon_{c}^{(n)}}\right\rVert_{\infty}$ then, writing $\hat{u}=\hat{\beta}_{n}^{\mathrm{lasso}}(f)-\beta^{\star}(f)$, we have $\hat{u}\in\mathcal{C}(S^{\star};3)$ and {ceqn}

[TABLE]

Proof

This is just a reformulation of the reasoning on p. 298 in (Tibshirani et al., 2015) with a slightly sharper upper bound. The vector $\hat{\nu}$ at the right-hand side of their Eq. (11.23) can be replaced by $\hat{\nu}_{S}$. For the sake of completeness, we provide the details.

In the proof we use the short-cuts $\beta^{\star}=\beta^{\star}(f)$ and $\hat{\beta}_{n}^{\mathrm{lasso}}=\hat{\beta}_{n}^{\mathrm{lasso}}(f)$. Recall $\epsilon_{c}^{(n)}=f_{c}^{(n)}-H_{c}\beta^{\star}(f)$ and define {ceqn}

[TABLE]

Because $G(\hat{u})\leqslant G(0)$, we have {ceqn}

[TABLE]

From the triangle inequality {ceqn}

[TABLE]

implying that {ceqn}

[TABLE]

From Hölder’s inequality, we get {ceqn}

[TABLE]

which leads to {ceqn}

[TABLE]

Consequently, because $\left\lVert{\textstyle H_{c}^{T}\epsilon_{c}^{(n)}}\right\rVert_{\infty}/n\leqslant\lambda/2$ by assumption, we obtain {ceqn}

[TABLE]

The right-hand side must be nonnegative, whence $\left\lVert{\textstyle\hat{u}_{\overline{S^{\star}}}}\right\rVert_{1}\leqslant 3\left\lVert{\textstyle\hat{u}_{S^{\star}}}\right\rVert_{1}$, i.e., $\hat{u}\in\mathcal{C}(S;3)$. The bound in (40) follows as well.

On the event $E$, the conclusion of Lemma 7 is valid because the bound on $\left\lVert{\textstyle H_{c}^{T}\epsilon_{c}^{(n)}}\right\rVert_{\infty}$ in (38) and the assumption on $\lambda$ in Theorem 3.2 together imply that $\lambda\geqslant 2\left\lVert{\textstyle H_{c}^{T}\epsilon_{c}^{(n)}}\right\rVert_{\infty}/n$. The cone property of Lemma 7 yields $\hat{u}\in\mathcal{C}(S^{\star};3)$ so that {ceqn}

[TABLE]

Thanks to (37) and Lemma 7, and since $\left\lvert{\textstyle S^{\star}}\right\rvert=\ell^{\star}$, we get {ceqn}

[TABLE]

It follows that $\left\lVert{\textstyle\hat{u}_{S^{\star}}}\right\rVert_{1}\leqslant 12\ell^{\star}\lambda/\gamma^{\star}$. In combination with (41), we find (39).

Step 6. — Equation (32) gave a bound on the estimation error involving three terms. On the event $E$, these terms were shown to be bounded in (33), (34), and (39). It follows that, on $E$, we finally have {ceqn}

[TABLE]

Divide by $n$ and use $48\sqrt{2}<68$ to obtain (9). ∎

Appendix D Proof of Theorem 3.3

Recall that $S^{\star}=\{j=1,\ldots,m:\beta_{j}^{\star}(f)\neq 0\}$ with $\ell^{\star}=\left\lvert{\textstyle S^{\star}}\right\rvert$ and that $\overline{S^{\star}}=\{1,\ldots,m\}\setminus S^{\star}$. Further, $H_{c,S^{\star}}$ is the $n\times\ell^{\star}$ matrix having columns $H_{c,k}$ for $k\in S^{\star}$, where $H_{c,k}$ is the $k$-th column of $H_{c}$.

Step 1. — We first establish some (non-probabilistic) properties of $\hat{\beta}_{n}^{\mathrm{lasso}}(f)$. To this end, we consider the linear regression of the non-active control variates on the active ones: for $k\in\overline{S^{\star}}=\{j=1,\ldots,m:\beta_{j}^{\star}(f)=0\}$, this produces the coefficient vector {ceqn}

[TABLE]

Further, we consider the OLS oracle estimate $\hat{\beta}_{n}^{\star}$, which is the OLS estimator based upon the active control variables only, i.e., {ceqn}

[TABLE]

Our assumptions will imply that, with large probability, $H_{c,S^{\star}}$ has rank $\ell^{*}$, in which case {ceqn}

[TABLE]

The following lemma provides a number of (non-probabilistic) properties of $\hat{\beta}_{n}^{\mathrm{lasso}}(f)$, given certain conditions on $H_{c}$ and $\epsilon_{c}^{(n)}$. Recall that a norm $\left\lVert{\textstyle\,\cdot\,}\right\rVert$ on $\mathbb{R}^{p}$ induces a matrix norm on $\mathbb{R}^{p\times p}$ via $\left\lVert{\textstyle A}\right\rVert=\sup\{\left\lVert{\textstyle Au}\right\rVert:u\in\mathbb{R}^{p},\left\lVert{\textstyle u}\right\rVert=1\}$ for $A\in\mathbb{R}^{p\times p}$.

Lemma 8

If $H_{c,S^{\star}}$ has rank $\ell^{\star}$ and if there exists $\kappa\in(0,1]$ such that {ceqn}

[TABLE]

then the minimizer $\hat{\beta}_{n}^{\mathrm{lasso}}(f)$ in (5) is unique, with support $\operatorname{supp}(\hat{\beta}_{n}^{\mathrm{lasso}}(f))\subset S^{\star}$, and it satisfies

[TABLE]

Proof

The proof of the previous result is actually contained in Tibshirani et al. (2015). The uniqueness of the LASSO solution and the property that it does not select inactive covariates follows directly from the proof of their Theorem 11.3. The only difference is that, in our case, the inequality (43) is an assumption whereas in Tibshirani et al. (2015) it is a property of the Gaussian fixed design model. The approach in Tibshirani et al. (2015) is based upon checking the strict dual feasibility condition. The bound (44) is Eq. (11.37) in Tibshirani et al. (2015).

We slightly modify Lemma 8 to make the conditions (42) and (43) easier to check and to make the bound (44) easier to use.

Lemma 9

If there exists $\nu>0$ such that {ceqn}

[TABLE]

and if there exists $\kappa\in(0,1]$ such that {ceqn}

[TABLE]

then the minimizer $\hat{\beta}_{n}^{\mathrm{lasso}}(f)$ in (5) is unique, with support satisfying $\operatorname{supp}(\hat{\beta}_{n}^{\mathrm{lasso}}(f))\subset S^{\star}$, and it holds that {ceqn}

[TABLE]

Proof

By (45), the smallest eigenvalue of the $\ell^{\star}\times\ell^{\star}$ matrix $H_{c,S^{\star}}^{T}H_{c,S^{\star}}$ is positive, so that it is invertible and $H_{c,S^{\star}}$ has rank $\ell^{\star}$.

We show that (46) implies (42). For each $k\in\overline{S^{\star}}$, the vector $\hat{\theta}_{n}^{(k)}$ has length $\ell^{\star}$, so that

[TABLE]

Because $\hat{\theta}_{n}^{(k)}$ is an OLS estimate, using that the largest eigenvalue of $(H_{c,S^{\star}}^{T}H_{c,S^{\star}})^{-1}$ being bounded from above by $(n\nu)^{-1}$, we obtain

[TABLE]

Since $\left\lVert{\textstyle x}\right\rVert_{2}\leqslant\sqrt{m}\,\left\lVert{\textstyle x}\right\rVert_{\infty}$ for $x\in\mathbb{R}^{m}$, we can conclude that {ceqn}

[TABLE]

Combining the two bounds, we find that (46) indeed implies (42).

Next we show that (47) implies (43). For $k\in\overline{S^{\star}}$, we have {ceqn}

[TABLE]

Using (42) and (47) we deduce (43).

The conditions of Lemma 8 have been verified, and so its conclusion holds. We simplify the two terms in the upper bound (44). First, we use that {ceqn}

[TABLE]

Second, for any matrix $A\in\mathbb{R}^{p\times p}$, we have $\left\lVert{\textstyle A}\right\rVert_{\infty}\leqslant\sqrt{p}\left\lVert{\textstyle A}\right\rVert_{2}$ (e.g., (Horn and Johnson, 2012, page 365)), and this we apply to $(H_{c,S^{\star}}^{T}H_{c,S^{\star}})^{-1}$. In this way, the upper bound in (44) is dominated by

[TABLE]

since the largest eigenvalue of $(H_{c,S^{\star}}^{T}H_{c,S^{\star}})^{-1}$ is at most $(n\nu)^{-1}$. Use (47) to further simplify the right-hand side, yielding (48).

Step 2. — Let $\delta\in(0,1)$ and $n=1,2,\ldots$. In a similar way as in the proof of Theorem 3.1, we construct an event of probability at least $1-\delta$. This time, we need five inequalities to hold simultaneously.

• •

Because $\epsilon\in\mathcal{G}(\tau^{2})$, with probability at least $1-\delta/5$, {ceqn}

[TABLE]

• •

In view of (Boucheron et al., 2013, Lemma 2.2) and Assumption 2, we have $h_{k}\in\mathcal{G}(U_{h}^{2})$ for all $k=1,\ldots,m$. Hence we can apply Lemma 1 with $p=m$ to get that, with probability at least $1-\delta/5$, {ceqn}

[TABLE]

• •

By virtue of Assumptions 1 and 2, we can apply Lemma 4 to have $h_{k}\epsilon\in\mathcal{G}(CU_{h}^{2}\tau^{2})$, where $C=9/2$. Hence we can apply Lemma 1 to get that, with probability at least $1-\delta/5$, {ceqn}

[TABLE]

• •

Recall that $B^{\star}=\sup_{x\in\mathcal{X}}h_{S^{\star}}^{T}(x)G_{S^{\star}}^{-1}h_{S^{\star}}(x)$ with {ceqn}

[TABLE]

The assumption on $n$ easily implies that $n\geqslant 8B^{\star}\log(5\ell^{\star}/\delta)$. Applying Lemma 2 with $p=\ell^{\star}$, $g=h_{S^{\star}}$, and $\delta$ replaced by $\delta/5$, we find that, with probability at least $1-\delta/5$, {ceqn}

[TABLE]

• •

Finally, because $\left\lvert{\textstyle h_{j}(x)}\right\rvert\leqslant U_{h}$ for all $x\in\mathcal{X}$ and $j\in\{1,\ldots,m\}$ and because $P(h_{k}h_{j})=0$ for all $(k,j)\in\overline{S^{\star}}\times S^{\star}$, we have $h_{k}h_{j}\in\mathcal{G}(U_{h}^{4})$ for such $k$ and $j$, and thus, with probability at least $1-\delta/5$, {ceqn}

[TABLE]

By the union bound, the event, say $E$, on which (49), (50), (51), (52) and (53) are satisfied simultaneously has probability at least $1-\delta$. We work on the event $E$ for the rest of the proof.

Step 3. — On the event $E$, we have {ceqn}

[TABLE]

where $\alpha\in(0,1/2)$ is an absolute constant whose value will be fixed in Step 6(ii). We have

[TABLE]

and thus, by the Cauchy–Schwarz inequality and by (52), {ceqn}

[TABLE]

In view of (50), we have

[TABLE]

We thus get

[TABLE]

A sufficient condition for (54) is thus that the term in square brackets is at least $\alpha$, i.e.,

[TABLE]

Since $\ell^{\star}\geqslant 1$ and $U_{h}^{2}\geqslant\gamma^{\star\star}$, a condition of the form {ceqn}

[TABLE]

is thus sufficient, with much to spare, provided $\rho>2/(1/2-\alpha)$. In Step 6(ii), we will choose $\alpha$ in such a way that the constant $\rho=70$ appearing in the statement of the theorem is sufficient.

Step 4. — On the event $E$, we have

[TABLE]

Indeed, denote $A=\overline{S^{\star}}\times S^{\star}$, in virtue of (50) and (53), the left-hand side is bounded by {ceqn}

[TABLE]

which is (56).

Step 5. — On the event $E$, we have

[TABLE]

The proof is the same as the first part of the one (38).

Step 6. — We will verify that on the event $E$, the three assumptions of Lemma 9 are satisfied with $\kappa=1/2$ and $\nu=\alpha\gamma^{\star\star}$, with $\alpha$ as in Step 3.

1. (i)

Eq. (45) with $\nu=\alpha\gamma^{\star\star}$ is just (54). 2. (ii)

Eq. (46) with $\nu=\alpha\gamma^{\star\star}$ and $\kappa=1/2$ follows from (56) provided we have

[TABLE]

To check whether this is satisfied, we will make use of the elementary inequality222The convex parabola $x\mapsto ax^{2}-bx-c$ has zeroes at $x_{\pm}=(b\pm\sqrt{b^{2}+4ac})/(2a)$, and $x_{-}<0<x_{+}<\sqrt{b^{2}+4ac}/a$. {ceqn}

[TABLE]

with $x=\sqrt{n}$ and {ceqn}

[TABLE]

Sufficient is that $n=x^{2}$ is bounded from below by $(b^{2}+4ac)/a^{2}=(b/a)^{2}+4c/a$, which is

[TABLE]

But $\ell^{\star}\geqslant 1$ and $\gamma^{\star\star}\leqslant(1/\ell^{\star})\sum_{j\in S^{\star}}P(h_{j}^{2})\leqslant U_{h}^{2}$, so that a sufficient condition is that

[TABLE]

The constant $\rho$ in (55) must thus be such that

[TABLE]

The minimum of the right-hand side as a function of $\alpha\in(0,1/2)$ occurs at $\alpha=\sqrt{2}/3$ and is equal to $2/(1/2-\sqrt{2}/3)\approx 69.9$. Taking $\rho=70$ as in the assumption on $n$ is thus sufficient. 3. (iii)

Eq. (47) with $\kappa=1/2$ follows from (57), since

[TABLE]

by the assumed lower bound on $\lambda$. Indeed, since $\ell^{\star}\geqslant 1$ and $U_{h}^{2}\geqslant\gamma^{\star\star}$, the assumed lower bounds on $n$ imply that $n\geqslant 70\log(10/\delta)$, so that $2\sqrt{\log(10/\delta)/n}$ is bounded by $2/\sqrt{70}$; recall $C=9/2$. Since $4\cdot(3+2/\sqrt{70})\approx 12.9$, the assumed lower bound for $\lambda$ suffices.

Step 7. — By the previous step, the conclusions of Lemma 9 with $\kappa=1/2$ and $\nu=\alpha\gamma^{\star\star}$ hold on the event $E$, where $\alpha=\sqrt{2}/3$ was specified in Step 6(ii). The minimizer $\hat{\beta}_{n}^{\mathrm{lasso}}$ in (5) is thus unique and we have $\operatorname{supp}(\hat{\beta}_{n}^{\mathrm{lasso}}(f))\subset S^{\star}$.

To show the reverse inclusion, we need to verify that $\left\lvert{\textstyle\hat{\beta}_{n,k}^{\mathrm{lasso}}(f)}\right\rvert>0$ for all $k\in S^{\star}$. To this end, we apply (48) with $\kappa=1/2$ and $\nu=\alpha\gamma^{\star\star}$, which becomes

[TABLE]

For any $k\in S^{\star}$, we thus have

[TABLE]

But for $\alpha=\sqrt{2}/3$, we have approximately $5/(4\alpha)\approx 2.65$. Since $\min_{j\in S^{\star}}\left\lvert{\textstyle\beta_{j}^{\star}(f)}\right\rvert>3\sqrt{\ell^{\star}}\lambda/\gamma^{\star\star}$ by the assumed upper bound for $\lambda$, we find $\left\lvert{\textstyle\hat{\beta}_{n,k}^{\mathrm{lasso}}(f)}\right\rvert>0$, as required. ∎

Appendix E Proof of Theorem 3.4

Recall that the LSLASSO estimator is defined as an OLS estimate computed on the active variables selected by the LASSO based on a subsample of size $N\in\{1,\ldots,n\}$. Let $\hat{\beta}_{N}^{\mathrm{lasso}}(f)$ denote the LASSO coefficient vector in (5) based on the subsample $X_{1},\ldots,X_{N}$ and let

[TABLE]

denote the estimated active set of $\hat{\ell}=|\hat{S}_{N}|$ control variates. The LSLASSO estimate $\hat{\alpha}_{n}^{\mathrm{lslasso}}(f)$ based on the full sample $X_{1},\ldots,X_{n}$ is defined as the OLS estimator based on the control variates $h_{k}$ for $k\in\hat{S}_{N}$: writing $H_{\hat{S}_{N}}$ for the $n\times\hat{\ell}$ matrix with columns $(h_{k}(X_{i}))_{i=1}^{n}$ with $k\in\hat{S}_{N}$, we have {ceqn}

[TABLE]

Therefore, we can derive a concentration inequality by combining the support recovery property (Theorem 3.3) along with the concentration inequality for the OLS estimate (Theorem 3.1) using only the active control variates.

Let $\delta>0$ and $n\geqslant 1$. We construct an event with probability at least $1-\delta$ on which the support recovery property and the concentration inequality for the OLS estimate hold simultaneously. Recall that $S^{\star}=\operatorname{supp}(\beta^{\star}(f))$ is the true set of $\ell^{\star}=\left\lvert{\textstyle S^{\star}}\right\rvert$ active control variables.

• •

Thanks to Theorem 3.3, with probability at least $1-\delta/2$, {ceqn}

[TABLE]

Indeed, the conditions on $N$ and $\lambda$ in Theorem 3.4 are such that we can apply Theorem 3.3 with $n$ and $\delta$ replaced by $N$ and $\delta/2$, respectively.

• •

Thanks to Theorem 3.1, with probability at least $1-\delta/2$,

[TABLE]

where for any $S\subset\{1,\ldots,m\}$, $\hat{\alpha}_{n}^{\mathrm{ols}}(f,h_{S})$ is the OLS estimate of $P(f)$ based on the control variates $h_{S}$. Indeed, we apply Theorem 3.1 with $h$ and $\delta$ replaced by $h_{S^{\star}}$ and $\delta/2$, respectively. The required lower bound on $n$ is now $n\geqslant\max\left(18B^{\star}\log(8\ell^{\star}/\delta),75\ell^{\star}\log(8/\delta)\right)$. By assumption we have $N\geqslant 75[\ell^{\star}(U_{h}^{2}/\gamma^{\star\star})]^{2}\log(20\ell^{\star}/\delta)$. The required lower bound is already satisfied for $N$, and thus certainly by $n$.

By the union bound, the event on which (58) and (59) are satisfied simultaneously has probability at least $1-\delta$. On this event, we can, by definition of $\hat{\alpha}_{n}^{\mathrm{lslasso}}(f)$ and by (58), write the integration error as {ceqn}

[TABLE]

But the right-hand side is bounded by (59), yielding (13), as required.