Two Sobolev Spaces for Stein-Kernelized Monte Carlo

We want to estimate an expectation under a target distribution, using samples drawn from another distribution. Function evaluations may also contain simulation noise. The aim is to correct the sampling bias while using the structure of the integrand to reduce Monte Carlo error.

Henry Lam and Haofeng Zhang’s doubly robust Stein-kernelized estimator addresses this problem by combining a learned control functional with sample weights chosen through a Stein kernel. We will first develop their construction and explain its double robustness.

A similar idea of using different spaces for approximation and integration appears in Oates and Girolami (2016) for QMC control functionals. Here, we apply this idea to improve the DRSK convergence rate.

1. The problem and Lam–Zhang's estimator

What is observed, and what is the target?

Suppose a simulator produces independent observations \((X_i,F_i)\), where

The variable \(X\) is the input we want to integrate over; \(Y\) represents randomness inside the function evaluation. Define the conditional mean

Our target is its expectation under a specified distribution \(\pi\):

For example, \(X\) could describe a model parameter and \(F\) the output of a stochastic simulation at that parameter. The available runs use a convenient sampling distribution \(q\), while the scientific question asks for an average under \(\pi\). Here and below, \(\pi\) and \(q\) denote the distributions of \(X\), and \(dx\) is the volume measure on its state space \(M\).

The setting has three useful pieces of structure:

• The target density is known up to normalization. We can evaluate an unnormalized density \(\gamma\) with \(\pi=\gamma/Z\), and its score \(\nabla\log\pi=\nabla\log\gamma\). The construction uses this score to build Stein kernels.
• The conditional observation law is shared. Both the target and sampling settings use the same distribution \(p(dy\mid x)\) for \(Y\) given \(X=x\). Thus \(\bar f(x)\) is the same function in both settings. This is the covariate-shift assumption.
• The samples cover the target. The density ratio \(\rho=\pi/q\) exists, and the analysis uses \(\rho\le B\). Coverage supplies data wherever the target assigns mass. The weighting algorithm uses a kernel matrix rather than evaluations of \(q\).

An ordinary sample average estimates \(q(\bar f)\). Our goal is \(\pi(\bar f)\). In addition to this distribution shift, the observations contain the zero-mean conditional fluctuation

Lam and Zhang combine two tools to address these two sources of error: control functionals learn the predictable part of \(F\), and black-box importance sampling balances the samples against the target.

Control functionals: learn a function with a known integral

Start with a Stein kernel \(k_\pi\) and its reproducing kernel Hilbert space \(\mathcal H_\pi\). Every function in this space has zero target mean:

Use the first \(m\) observations to fit \(s_m=c_m+g_m\), with \(g_m\in\mathcal H_\pi\), by kernel ridge regression in \(\mathbb R\oplus\mathcal H_\pi\). The regression kernel is \(1+k_\pi\), and the fitted function has the known integral \(\pi(s_m)=c_m\).

Now adjust each observation by

This preserves the target expectation. Its conditional mean is \(\bar f(x)-s_m(x)+\pi(s_m)\), which becomes nearly constant as the fit improves. Monte Carlo then averages a smaller residual. This is the control-functional mechanism for variance reduction.

Black-box importance sampling: balance a class of test functions

Write \(Z_i=X_{m+i}\) for the remaining \(N=n-m\) inputs. Choose nonnegative weights with \(\sum_iw_i=1\) by minimizing the Stein discrepancy

Because the Stein functions have zero target mean, this objective measures how closely the weighted samples reproduce their target integrals. The weights are obtained from a quadratic program, with a cap \(w_i\le B_0/N\) to control their concentration. They are chosen to integrate the test class accurately; estimating the density ratio point by point is unnecessary.

DRSK: apply the weights to the adjusted observations

Lam and Zhang’s doubly robust Stein-kernelized estimator combines the two steps:

The training split supplies \(s_m\); the independent second split supplies the weights and residual observations. In their construction, regression uses \(1+k_\pi\) and weighting uses \(k_\pi\). The common Stein space connects the two steps in the error analysis.

2. Why doubly robust?

There are two routes to an accurate estimate. A good fit makes the conditional mean of the adjusted observation nearly constant, so the answer depends less on where the samples lie. Good weights integrate the remaining residual accurately, so they can compensate for an imperfect fit. For example, if \(s_m=\bar f\), only the conditional noise remains; if \(s_m=0\), the construction becomes the BBIS estimate of the original observations.

The two routes appear as a product in the integration-error bound. To see this directly, take noiseless observations, so \(F_i=f(X_i)\) and \(\bar f=f\). We use this smooth-integration setting for the Sobolev construction in the rest of the post.

Write \(\widehat\nu_w=\sum_iw_i\delta_{Z_i}\) for the weighted empirical measure. The estimation error has a useful exact form:

Center the regression residual under the target:

Since both \(\widehat\nu_w\) and \(\pi\) have mass one, centering leaves the error unchanged. If \(h_m\) belongs to the Stein RKHS used for the weights, the reproducing property gives

Cauchy–Schwarz then separates the two tasks:

This explains the double robustness: either factor can make the product small when the other remains controlled. Learning and weighting can also improve together.

With noisy observations, the same calculation applies to \(\bar f-s_m\), and the estimation error has the additional term \(\sum_i\widehat w_i\varepsilon_{m+i}\). Its conditional mean is zero. With bounded weights, its mean squared contribution is \(O(M_0/N)\), where \(M_0=\mathbb E_q[\varepsilon^2]<\infty\) measures the observation noise. The smooth integration term and the conditional-noise term therefore have separate roles.

Under their source and coverage assumptions, Lam and Zhang obtain the bound

The source exponent \(r\) describes how well the conditional mean can be approximated by the regression kernel; \(C_f\) contains its source norm. For \(r>1/2\), the smooth term decays faster than the canonical Monte Carlo MSE rate \(n^{-1}\). Their stronger spectral assumptions yield a further \(o(1)\) factor on that smooth term. These results make the residual–discrepancy pairing a useful starting point for designing an estimator.

We now apply this two-space viewpoint to the DRSK residual–discrepancy bound. For an integer \(s\ge0\), the Sobolev norm in \(H^s\) controls the function and its square-integrable derivatives through order \(s\); fractional orders describe intermediate regularity. Increasing \(s\) puts a stronger penalty on rapidly varying components of a function.

Kernel ridge regression will learn in a smoother Sobolev space \(H^\tau\), while weighting uses a Stein space of order \(a<\tau\). An interpolation inequality connects the two. The following spectrum explains how to choose a norm that both controls the residual and accommodates point samples.

3. A spectrum of norms

Why does the choice of norm matter? A signed measure \(\eta\) acts on a test function \(g\) through \(\eta(g)=\int g\,d\eta\). A dual norm measures how large this action can be over a specified class of test functions:

The test class determines which differences between distributions remain visible. The following spectrum moves from local density information toward averages of regular functions.

At the density end, point masses matter. An empirical measure is a finite sum of Dirac masses. Against a smooth target density on a positive-dimensional manifold, their difference has no \(L^2(dx)\) density. Correspondingly, evaluation at a point is not a bounded functional on \(L^2\): an increasingly narrow bump can have small \(L^2\) norm and a large value at that point.

Total variation also sees the finite sample support exactly. If \(A=\{Z_1,\ldots,Z_N\}\) and \(\pi\) has a density, then \(\widehat\nu_w(A)=1\) and \(\pi(A)=0\). With the convention \(d_{\mathrm{TV}}(\mu,\nu)=\sup_A\lvert\mu(A)-\nu(A)\rvert\), the distance is one.

Geometric and smooth tests average over that local structure. The Wasserstein distance \(W_1\) compares measures through 1-Lipschitz test functions. Nearby points have nearby function values, and an i.i.d. empirical measure converges to its sampling distribution in this metric on a compact space. The convergence rates reflect the dimension and geometry of the distribution; Weed and Bach develop this dependence in detail.

Sobolev tests provide another way to control regularity. On a compact \(d\)-dimensional manifold, the embedding

makes point evaluation continuous. Thus \(\delta_x\) belongs to the dual space \(H^{-a}(M)\). An empirical measure can now be paired with the residual using a finite, useful norm.

Within this Sobolev scale, increasing \(a\) imposes more smoothness on the test functions and gives a weaker dual norm. For \(a_2>a_1>d/2\),

The same choice has a complementary effect on the residual: reducing \(a\) makes its \(H^a\) norm easier to control. This is the balance we want to exploit.

The regression estimates explain why. For a suitable smoothness condition and regularization, the KRR residual satisfies the two endpoint bounds

An intermediate \(H^a\) norm inherits part of the faster \(L^2\) rate. We choose \(a\) above the point-evaluation threshold and below the fitting smoothness:

For each fixed such \(a\), the weight argument below supplies an \(O(N^{-1})\) bound for the squared discrepancy. The gain then comes from transferring the regression accuracy to \(H^a\).

4. Two spaces, one estimator

We now specify the setting and the estimator. Let \(M\) be a compact, connected, smooth Riemannian manifold without boundary, of dimension \(d\ge1\). The target and sampling densities \(\pi,q\) are smooth and strictly positive. Let \(f:M\to\mathbb R\) be smooth, and assume a bounded density ratio

Fix \(d/2<a<\tau\). For \(b\in\{a,\tau\}\), use a Stein kernel \(k_{\pi,b}\) whose RKHS is norm-equivalent to the mean-zero Sobolev space:

The quotient identifies functions that differ by a constant. Its norm is

We will construct these kernels in the next section. Their diagonals are bounded, and adjoining the constant kernel gives

The approximation assumption is expressed through the integral operator

For the main result, assume the source condition \(f\in\operatorname{Range}(T_{\tau,q})\): there is a \(v\in L^2(q)\) with \(f=T_{\tau,q}v\). This is the condition that yields the two KRR endpoint rates above.

The condition has a useful spectral interpretation. In a complete orthonormal eigenbasis of \(L^2(q)\), write \(T_{\tau,q}v_j=\sigma_j v_j\) with \(\sigma_j>0\) and \(f=\sum_j f_j v_j\). The source condition is then

It requires sufficiently little energy in directions with small kernel eigenvalues. This connects the integrand’s regularity to the particular fitting kernel: choosing \(\tau\) also changes the eigenvalues against which its coefficients are measured.

Both steps are finite-dimensional computations. If \(K_{\mathrm{cf}}=(k_{\mathrm{cf}}(X_i,X_j))_{i,j\le m}\) and \(\mathbf f=(f(X_1),\ldots,f(X_m))^\top\), then

The weighting step is a convex quadratic program with a positive semidefinite kernel matrix.

Rate comparison. In this noiseless setting, with the same fitting kernel and source condition \(r=1\), using the same space for both steps (\(a=\tau\)) gives the MSE bound \(O(n^{-3/2})\). Choosing \(d/2<a<\tau\) improves this bound to \(O(n^{-2+a/(2\tau)})\); for example, \(\tau=2a\) gives \(O(n^{-7/4})\).

5. Constructing the Sobolev–Stein kernels

The second-order Riemann–Stein operator used by Barp, Oates, Porcu, and Girolami is

It depends on the target through its log-density gradient. Integration by parts on the boundaryless manifold gives

Choose a scalar base kernel \(K_b\) with \(\mathcal H(K_b)\simeq H^{b+2}(M)\) and the differentiability needed for the following expression. Apply the operator to both kernel arguments:

The resulting kernel generates the image of the base RKHS under \(\mathcal L_\pi\). Elliptic regularity identifies that image with the mean-zero functions in \(H^b(M)\): the operator maps \(H^{b+2}\) continuously into \(H^b\), and the equation \(\mathcal L_\pi u=h\) has a solution, unique up to a constant, for mean-zero \(h\in H^b\). The solution estimate gives the equivalence of the induced RKHS norm and the Sobolev quotient norm.

This construction explains both useful properties at once. Kernel sections have known zero target mean, and their finite evaluations can be computed by differentiating a base kernel and using the target score.

A spectral picture

The two-derivative shift is especially transparent in the eigenbasis of the weighted Laplacian. Define

with an orthonormal eigenbasis in \(L^2(\pi)\). Connectedness gives \(\mu_0=0\) with constant \(\phi_0=1\); for \(j\ge1\), \(\mu_j>0\) and \(\pi(\phi_j)=0\). Elliptic regularity gives

A corresponding pair of kernels is

Applying \(\mathcal L_\pi\) to each argument multiplies the nonconstant term by \((-\mu_j)^2\) and removes the constant. The base kernel has Sobolev order \(b+2\); its Steinized kernel has order \(b\) and zero target mean.

For a concrete geometric base kernel, smoothly embed \(M\) into a Euclidean space and restrict a Matérn kernel to that embedded manifold. A Matérn smoothness parameter \(\nu=b+2-d/2\) yields the required Sobolev order \(b+2\) on \(M\). The restriction results of Fuselier and Wright provide this construction. The resulting kernel is then Steinized through the differential formula above.

6. The proof

The proof combines Lam and Zhang’s KRR estimates and their capped-BBIS comparison (Lemma 29) with Sobolev interpolation. The interpolation step measures the regression residual in the lower-order space. The weighting argument then controls the discrepancy of the optimized weights by comparing them with normalized importance weights.

Throughout the proof, \(C\) may change from line to line and depends only on the fixed quantities in the theorem.

Step 1. Pair the centered residual with the weighted measure

Recall \(e_m=f-s_m\) and \(h_m=e_m-\pi(e_m)\). Since \(e_m\in H^\tau(M)\subset H^a(M)\), its centered representative belongs to \(\mathcal H_{\pi,a}\). The reproducing property gives

The fitted function depends on the first sample split; the weights depend on the second. Hence

Step 2. Interpolate the two regression bounds

Under \(f\in\operatorname{Range}(T_{\tau,q})\) and \(\lambda\asymp m^{-1/2}\), the KRR estimates in Corollaries 20 and 21 of Lam and Zhang give

Here the source exponent is \(r=1\), and equivalence of the regression RKHS and \(H^\tau\) converts the second bound to Sobolev notation.

For each realization of the training sample, the Sobolev interpolation theorem (Bergh and Löfström, 1976, Chapter 6) gives

Because \(q\) is continuous and strictly positive on compact \(M\), it has a positive lower bound. Thus \(\|e_m\|_{L^2(dx)}\le C\|e_m\|_{L^2(q)}\). Also \([h_m]=[e_m]\), so

Take expectations and apply Hölder with exponents \(1/(1-a/\tau)\) and \(\tau/a\):

This is the step where the gap between \(a\) and \(\tau\) enters the rate.

Step 3. Bound the optimized discrepancy

Let

Sobolev embedding and compactness give a uniform bound on point evaluations. Equivalently, for some finite \(\kappa_a\),

Stein’s identity gives \(\mathbb E_\pi\Psi_a(X)=0\). Changing measure yields the useful centered feature:

Now define normalized oracle weights on the second sample split:

These weights serve as a comparison in the proof. The algorithm computes \(\widehat w\) from the kernel matrix.

Consider the event

On this event,

The oracle weights are therefore feasible. By optimality of \(\widehat w\),

The norms here are in \(\mathcal H_{\pi,a}\). Independence of the samples and the zero mean in (5) remove the cross terms:

Consequently,

To handle the complementary event, use \(0\le\rho\le B\) and \(\mathbb E_q\rho=1\). Hoeffding’s inequality gives

For every feasible weight vector, nonnegativity and normalization give

It follows that

Combining (6) and (7), and absorbing the exponentially small term into the constant, yields

Step 4. Combine the two factors

Substitute (3) and (8) into (1):

For \(m\asymp N\asymp n\), this becomes

The proof uses the smoother space to learn the integrand, the lower-order space to pair the residual with the sample measure, and interpolation to connect them.

7. What the two smoothness levels control

The exponent is easy to read once the roles of \(a\) and \(\tau\) are separated.

The weighting level \(a\) makes point samples usable. The condition \(a>d/2\) gives bounded point evaluations and a finite kernel diagonal. It is exactly what the discrepancy proof needs. Reducing \(a\) toward this threshold also improves the residual exponent for fixed \(\tau\).

The fitting level \(\tau\) supplies interpolation room. Under the main source condition, the two regression endpoints are \(m^{-1}\) in squared \(L^2\) norm and \(m^{-1/2}\) in squared \(H^\tau\) norm. Measuring the residual at level \(a\) yields the intermediate exponent \(-1+a/(2\tau)\).

For fixed \(a\), increasing the fitting smoothness \(\tau\) brings the exponent closer to \(-2\) when the integrand satisfies the source condition for the resulting kernel. Each choice has its own kernel constants and source norm.

Finite regularity

The same argument also shows how the integrand’s regularity enters. More generally, retain the geometric and kernel assumptions and take an integrand satisfying, for a fixed \(\tau\),

With \(\lambda\asymp m^{-1/2}\), the two KRR bounds become

Interpolation gives

The weight estimate (8) is unchanged. Hence

This expression captures the choice of fitting smoothness: increasing \(\tau\) improves the interpolation factor, while the source exponent \(r\) records how well the particular integrand fits that kernel. For fixed \(r\) and the same fitting kernel, the gain over the same-space bound is \(\tfrac12-\tfrac{a}{2\tau}\) in the exponent.

The useful design freedom is to assign a space to each task. Fit where the function’s structure can be learned, balance where point samples define a controlled discrepancy, and use the relation between those spaces to bound the integral.

8. Beyond Sobolev: a Stein-generated Hilbert scale

The same analysis extends from Sobolev spaces to a Stein-generated Hilbert scale. A stable coefficient representation takes the place of Sobolev interpolation, and a bounded kernel diagonal supplies the point-evaluation control previously obtained from \(a>d/2\). The discrepancy proof carries over directly. We construct the two Stein RKHSs from shared features, derive the coefficient inequality, and combine it with the regression bounds.

We keep the independent sample splits, the bounded density ratio \(\rho=d\pi/dq\le B\), and the weight cap \(B_0/N\) with \(B_0>2B\). The state space is now denoted by \(\mathcal X\). It carries a second-order Stein operator \(\mathcal L_\pi\) satisfying

for potentials \(u\) in its Stein class. The Riemann–Stein operator from Section 5 has these properties.

Build the spaces from Stein features

Choose scalar potentials \(u_j\) in the Stein class and form the features

Assume these features form a Riesz basis of a closed subspace \(V\subset\{h\in L^2(q):\pi(h)=0\}\). This means that they give unique, stable coordinates: there are constants \(0<c_R\le C_R<\infty\) such that, for every finite coefficient sequence,

Let \(\omega_j\uparrow\infty\) be positive weights. For \(s\ge0\), define

Increasing \(s\) puts more weight on the coefficients with large \(\omega_j\). These nested spaces form a Hilbert scale. At level zero, the Riesz inequalities identify the coefficient norm with the \(L^2(q)\) norm on \(V\).

Choose \(0<a<\tau\) and assume bounded point evaluations at the two levels used by the algorithm:

Indeed, Cauchy–Schwarz gives

Thus \(\mathcal H_s^{\mathrm{St}}\) is an RKHS with kernel

Every function in this space has zero target mean. To see this for an infinite expansion, the partial sums converge in \(L^2(q)\), and the bounded density ratio makes \(\pi(h)=\mathbb E_q[\rho h]\) continuous in that norm. In particular, \(\pi(k_s^{\mathrm{St}}(x,\cdot))=0\).

Compute the kernels through Steinization

The potentials also provide scalar base kernels:

Assume these series define kernels in the Stein class in each variable, with enough regularity to apply the two Stein operators term by term. Choose the base kernels and operator so that they can be evaluated from the unnormalized target density, without its normalizing constant. Then

This supplies the same computational interface as before. KRR uses \(1+k_\tau^{\mathrm{St}}\), whose RKHS is \(\mathbb R\oplus\mathcal H_\tau^{\mathrm{St}}\) and whose fitted functions have known integrals. BBIS uses \(k_a^{\mathrm{St}}\) and minimizes the quadratic form

The feature representation describes the spaces; the differential formula supplies kernel evaluations without the target normalizing constant.

Interpolate the coefficients

The key inequality follows directly from the weighted coefficients. Set \(t=a/\tau\in(0,1)\). For \(h=\sum_jd_j\psi_j\in\mathcal H_\tau^{\mathrm{St}}\), Hölder’s inequality and the lower Riesz bound give

Now take the centered residual \(h_m=(f-s_m)-\pi(f-s_m)\). Suppose it belongs to \(\mathcal H_\tau^{\mathrm{St}}\) almost surely and satisfies the two endpoint bounds

Apply the coefficient inequality to each realization of the training sample, then apply Hölder to the expectation:

Pair the residual with the discrepancy

Because the weights sum to one, the error of the noiseless estimator is still

Cauchy–Schwarz and the independent sample splits therefore give

The discrepancy proof in Section 6 applies to the feature map \(\Psi_a(x)=k_a^{\mathrm{St}}(x,\cdot)\). Its two inputs are \(\mathbb E_q[\rho(X)\Psi_a(X)]=0\) and \(\|\Psi_a(x)\|^2\le\kappa_a^2\). Comparing with the normalized oracle weights on \(E_N=\{N^{-1}\sum_i\rho(Z_i)\ge1/2\}\), and bounding its complement, yields

Combining the two estimates gives the general bound:

Insert the KRR rates

For a concrete rate, define the regression integral operator

In the noiseless setting, assume

The same KRR estimates used earlier give

Centering preserves these rates. Since \(\mathbb E_q\rho^2\le B\),

In the direct-sum RKHS, taking the target mean extracts the constant component. Removing it gives

Thus we can take \(b_m\lesssim m^{-r}\) and \(c_m\lesssim m^{-r+1/2}\). Substitution gives

With comparable split sizes, the rate is \(n^{-1-r+a/(2\tau)}\). Setting \(r=1\) recovers \(n^{-2+a/(2\tau)}\) in the general Stein scale.

For noisy observations, use the centered residual of \(\bar f-s_m\) and its corresponding regression bounds \(b_m,c_m\). Conditional mean-zero noise removes the cross term with the integration error. Independence across the second split and the cap \(\widehat w_i\le B_0/N\) give

The resulting bound is \(C N^{-1}b_m^{1-a/\tau}c_m^{a/\tau}+B_0^2M_0/N\).

Recover the Sobolev construction

Return to the compact manifold and the weighted-Laplacian eigenpairs from Section 5. For \(j\ge1\), choose

The eigenfunctions are orthonormal in \(L^2(\pi)\). Smooth positive densities on compact \(M\) make the \(L^2(\pi)\) and \(L^2(q)\) norms equivalent, so they form the required Riesz basis of \(V=\{h\in L^2(q):\pi(h)=0\}\). The scale and its base kernel become

Sobolev embedding supplies the bounded diagonal when \(a>d/2\), and the potentials account for the two-derivative shift in the base kernel. The coefficient proof shows how the same estimator and rate extend to other Stein-generated scales with stable coordinates, bounded kernel diagonals, and a computable Steinization.

References

• Henry Lam and Haofeng Zhang. Doubly Robust Stein-Kernelized Monte Carlo Estimator: Simultaneous Bias-Variance Reduction and Supercanonical Convergence. Journal of Machine Learning Research, 24(85):1–58, 2023. The DRSK estimator, KRR bounds, and capped-BBIS comparison used in the proof.
• Chris Oates and Mark Girolami. Control Functionals for Quasi-Monte Carlo Integration. Proceedings of AISTATS, PMLR 51:56–65, 2016. Separate approximation and integration spaces for control functionals, with residual–discrepancy error bounds.
• Alessandro Barp, Chris J. Oates, Emilio Porcu, and Mark Girolami. A Riemann–Stein Kernel Method. Bernoulli, 28(4):2181–2208, 2022. The second-order Stein construction on a manifold.
• Edward Fuselier and Grady B. Wright. Scattered Data Interpolation on Embedded Submanifolds with Restricted Positive Definite Kernels: Sobolev Error Estimates. SIAM Journal on Numerical Analysis, 50(3):1753–1776, 2012. Sobolev kernels obtained by restriction to an embedded manifold.
• Jöran Bergh and Jörgen Löfström. Interpolation Spaces: An Introduction. Springer, 1976. Chapter 6: interpolation of Sobolev and Besov spaces.
• Jonathan Weed and Francis Bach. Sharp Asymptotic and Finite-Sample Rates of Convergence of Empirical Measures in Wasserstein Distance. Bernoulli, 25(4A):2620–2648, 2019. The geometric side of the spectrum of norms.

@misc{zhu2026sobolevstein,
  author = {Zhu, Youheng},
  title = {Two {Sobolev} Spaces for {Stein}-Kernelized {Monte Carlo}},
  year = {2026},
  howpublished = {Research blog},
  url = {https://zhuyouheng.github.io/notes/two-space-stein/}
}