To improve the accuracy of Monte Carlo estimation of expectations, a set of zero-mean feature functions, known as control variates, can be used. They can be used as feature functions for linear regression of the target function, and we can obtain an unbiased and variance-reduced estimate using its residual. One known way to construct such functions is a method using an equality called Stein’s identity, but these functions are not sufficient for the case where the target distribution is multimodal. We propose a different approach to constructing these zero-mean functions based on distribution approximation and the density ratio. We demonstrate that combining the functions constructed by these two strategies can effectively reduce the estimation variance for a bimodal distribution.

1 Introduction↩︎

The Monte Carlo (MC) method [1], [2] is one of the common tools in data science for handling uncertainty. It repeatedly simulates the system of interest and its uncertainty, and then aggregates the results to obtain its general properties. Specifically, in this study, we consider the problem of computing the expectation of a function over a probability distribution. The expectation to be computed is expressed as \[\begin{align} I=\int_{\mathbb{R}^D} f(\boldsymbol{x})p(\boldsymbol{x})d\boldsymbol{x}, \end{align}\] where \(p(\boldsymbol{x})\) represents the target probability density function in the \(D\)-dimensional space \(\mathbb{R}^D\), and \(f\) represents the target function whose expectation is to be computed. In the MC method, we draw random samples \(\boldsymbol{x}^{(1)}, \ldots, \boldsymbol{x}^{(T)}\) from the distribution \(p(\boldsymbol{x})\) using a method such as Markov Chain Monte Carlo (MCMC) [3], [4], and estimate \(I\) by \[\begin{align} \tilde{I}\equiv \frac{1}{T}\sum_{t=1}^{T}f(\boldsymbol{x}^{(t)}). \end{align}\] If \(\boldsymbol{x}^{(1)}, \ldots, \boldsymbol{x}^{(T)}\) are independent, the estimate is unbiased (\(\mathbb{E}[\tilde{I}]=I\)) and its variance is computed as \(V[\tilde{I}] = V_{\boldsymbol{x}\sim p}[f(\boldsymbol{x})]/T\). Here, \(\mathbb{E}[\cdot]\) and \(V[\cdot]\) represent the expectation and variance, respectively, and their subscripts denote the distribution that the random variable follows.

2 Variance Reduction of MC estimation by Control Variates↩︎

The accuracy of the estimation is governed by its variance. A simple way to reduce the estimation variance is to increase the number of samples \(T\). However, other methods are also known as more direct ways to reduce the variance.

Let us assume that we have a set of functions \(\varphi_m(\boldsymbol{x})\) whose expectations \(\mu_m=\mathbb{E}_{\boldsymbol{x}\sim p}\!\left[\varphi_m(\boldsymbol{x})\right]\) are known in advance, and whose linear combination approximates the target function as \[\begin{align} f(\boldsymbol{x}) = \Bigg({\sum_{m}b_m \varphi_m(\boldsymbol{x})}\Bigg) + C + \varepsilon(\boldsymbol{x}),\label{eq:linear-model} \end{align}\tag{1}\] where \(m\) runs over the index of the functions. Then, using its residual \(\varepsilon(\boldsymbol{x})\), we can estimate \(\tilde{I}\) as \[\begin{align} \tilde{I}&\equiv \Bigg( {\sum_{m}b_m\mu_m+ C} \Bigg) + \frac{1}{T}\sum_{t=1}^{T}\varepsilon(\boldsymbol{x}^{(t)})\nonumber\\ &=\frac{1}{T}\sum_{t=1}^{T}\Big(f(\boldsymbol{x}^{(t)}) - \sum_{m}b_m (\varphi_m(\boldsymbol{x}^{(t)}) - \mu_m)\Big).\label{eq:cv-estimation} \end{align}\tag{2}\] This estimation is still unbiased, and its variance is \(V[\tilde{I}] = V_{\boldsymbol{x}\sim p}[\varepsilon(\boldsymbol{x})]/T\), which is improved from the previous case when the model parameters \(b_m\) are properly chosen.

The parameters can be obtained by a linear regression. Specifically, we divide the samples into two disjoint sets \(\mathcal{X}_1\) and \(\mathcal{X}_2\), perform the linear regression only with \(\mathcal{X}_1\), and compute the estimate (Eq. (2 )) using \(\mathcal{X}_2\). Then, the estimation is still unbiased, and we can obtain the variance reduction. This is known as a method of control variates (CV) [1], [2], [5], [6]. It is also worth noting that the zero-mean functions can also be used for the design of the herding dynamics [7]–[9] to improve the sample generation.

We study how these feature functions are constructed, because the performance of this method is determined by how well the linear model approximates the target. We especially consider constructing zero-mean functions that satisfy \(\mathbb{E}_{\boldsymbol{x}\sim p}\!\left[f(\boldsymbol{x})\right]=0\), because those with known nonzero expectations \(\mu_m\neq 0\) are easily reduced to zero-mean functions by \(\varphi_m(\boldsymbol{x})\leftarrow \varphi_m(\boldsymbol{x})-\mu_m\).

3 Zero-mean Feature Functions Using Stein’s Identity↩︎

It is known in the literature that \(\mathbb{E}_{\boldsymbol{x}\sim p}\!\left[\varphi(\boldsymbol{x})\right]=0\) holds for \[\begin{align} \varphi(\boldsymbol{x}) = \nabla \log p(\boldsymbol{x}) \cdot \boldsymbol{\psi}(\boldsymbol{x}) + \nabla\cdot \boldsymbol{\psi}(\boldsymbol{x}), \label{eq:stein} \end{align}\tag{3}\] where \(\boldsymbol{\psi}\) is a vector-valued function under a suitable assumption. This property is also known as Stein’s identity and can be used to construct the zero-mean function \(\varphi_m\) by designing \(\boldsymbol{\psi}_m\) [6], [9], [10]. To prove the relation, we can use the following integral over a smooth region \(\Omega \subset\mathbb{R}^D\): \[\begin{align} &\int_\Omega p(\boldsymbol{x})\varphi(\boldsymbol{x}) d\boldsymbol{x}\nonumber\\ &=\int_\Omega p(\boldsymbol{x}) \nabla \log p(\boldsymbol{x}) \cdot \boldsymbol{\psi}(\boldsymbol{x}) d\boldsymbol{x}+ \int_\Omega p(\boldsymbol{x}) \nabla\cdot \boldsymbol{\psi}(\boldsymbol{x}) d\boldsymbol{x}\nonumber\\ &=\int_\Omega \nabla \cdot \big(p(\boldsymbol{x})\boldsymbol{\psi}(\boldsymbol{x})\big)=\int_{\partial \Omega} p(\boldsymbol{x})\boldsymbol{\psi}(\boldsymbol{x})^\top\boldsymbol{n}(\boldsymbol{x}) dS, \end{align}\] where \(\boldsymbol{n}(\boldsymbol{x})\) denotes the normal vector on the boundary \(\partial \Omega\), and the last equality follows from the divergence theorem. We can derive the zero-mean property under a decay assumption on \(p(\boldsymbol{x})\) as \(\|\boldsymbol{x}\|\to\infty\).

However, when the distribution \(p\) has multiple separated modes and the region \(\Omega\) covers only one of them, the integral can also be very small because \(p(\boldsymbol{x})\) in the last expression becomes small over the boundary \(\partial\Omega\). When the target function \(f\) has different expectation values across modes, but all feature functions have nearly zero expectations within each mode, accurately representing this discrepancy in the linear model (Eq. (1 )) requires very large coefficients \(b_m\), which leads to instability in the CV estimator Eq. (2 ).

4 Zero-mean Feature Functions Based on Probability Density Ratio↩︎

To avoid this instability, we need another class of zero-mean functions that have nonzero expectations for each mode at the same time. To this end, we propose a method based on density approximation and the probability density ratio as follows.

Let us assume a reference probability distribution \(R(\boldsymbol{x})\) and a zero-mean function \(\tilde{\varphi}(\boldsymbol{x})\) under \(R\); \(\mathbb{E}_{\boldsymbol{x}\sim R}\!\left[\tilde{\varphi}(\boldsymbol{x})\right]=0\). Then, we can define the function \(\varphi(\boldsymbol{x})\) with the zero-mean property \(\mathbb{E}_{\boldsymbol{x}\sim p}\!\left[\varphi(\boldsymbol{x})\right]=0\) as \[\begin{align} \varphi(\boldsymbol{x}) &= w(\boldsymbol{x})\tilde{\varphi}(\boldsymbol{x})\\ w(\boldsymbol{x})&\propto\frac{R(\boldsymbol{x})}{p(\boldsymbol{x})}, \end{align}\] where the property follows directly from the assumptions. Importantly, changing the constant factor of the coefficient function \(w(\boldsymbol{x})\) does not affect the zero-mean property, so this method can be applied to the case where the target distribution \(p(\boldsymbol{x})\) can only be computed up to a constant factor.

Let us further assume that the reference is a mixture of probability density functions \(r_k(\boldsymbol{x})\), each approximating a mode of the target, as \[\begin{align} R(\boldsymbol{x}) \propto \sum_k \alpha_k r_k(\boldsymbol{x}) \end{align}\] and construct the function \(\tilde{\varphi}\) as \[\begin{align} \tilde{\varphi}(\boldsymbol{x}) = \sum_k c_k \frac{r_k(\boldsymbol{x})}{R(\boldsymbol{x})}, \label{eq:feature-basis} \end{align}\tag{4}\] where \(k\) runs over the indices of the mixture components. Then, we can control the expectation of \(\tilde{\varphi}\) over each mode by the parameters \(c_k\). At the same time, we can easily keep the overall expectation to zero, which is computed as \(\mathbb{E}_{\boldsymbol{x}\sim R(\boldsymbol{x})}\!\left[\tilde{\varphi}(\boldsymbol{x})\right]=\sum_k c_k\) because of the assumption that \(r_k\) are density functions. Since we assume that \(r_k\) approximates each mode, the weight \(w(\boldsymbol{x})\) becomes nonzero in the corresponding mode’s region. Then, the nonzero expectation for each mode is inherited by the obtained function \(\varphi\).

5 Numerical Examples↩︎

We demonstrate how the proposed zero-mean feature functions contribute to variance reduction through numerical examples with a bimodal distribution. Let us consider the two-dimensional distribution for \(\boldsymbol{x}=(x_1, x_2)\), as illustrated in Fig. 1. This is based on an energy function \(E(x_1, x_2)\), and its density function \(p(x_1, x_2)\) is defined as \[\begin{align} p(x_1, x_2)&\propto \exp\Big( {-E(x_1, x_2)}\Big),\label{eq:example-distribution}\\ E(x_1,x_2)&=3e(x_1)+3e(x_2)+x_1x_2,\\ e(x)&=x^4 - 2x^2, \end{align}\tag{5}\] where \(e(x)\) represents a double-well potential for a single variable, which has an energy bump at the origin.

After generating \(\mathcal{X}_1\) and \(\mathcal{X}_2\), each consisting of 1000 samples, by standard MCMC, we randomly construct \(M_{\mathrm{S}}\) feature functions using Stein’s identity, as described in Section 3. We define the vector-valued function \(\boldsymbol{\psi}_i\) as \[\begin{align} \boldsymbol{\psi}_i(\boldsymbol{x}) = \sin{\Big(2\pi\big(\boldsymbol{k}_i^\top(\boldsymbol{x}/\sigma)+\theta_i\big)\Big)}\,\boldsymbol{d}_i , \end{align}\] where \(\boldsymbol{d}_i, \boldsymbol{k}_i\in\mathbb{R}^2\) and \(\theta_i\in[0, 2\pi)\) are randomly drawn parameters and \(\sigma\) is the scale parameter. Specifically, \(\boldsymbol{d}_i\) is normalized to \(\|\boldsymbol{d}_i\|=1\) so that its direction is uniformly distributed, \(\boldsymbol{k}_i\) is drawn from the multivariate standard normal distribution with zero mean vector and identity covariance matrix, and \(\theta_i\) follows the uniform distribution on \([0, 2\pi)\). Then, we construct a zero-mean function \(\varphi_i\) for each \(i\) by Eq. (3 ) using \(\boldsymbol{\psi}_i\) as \(\boldsymbol{\psi}\).

We also constructed \(M_{\mathrm{R}}\) ratio-based feature functions following the procedure described in Section 4. We used a reference distribution that has \(M_{\mathrm{R}}\) mixture components with equal weights \(\alpha_k = 1\), and the \(M_{\mathrm{R}}\) base functions defined by Eq. (4 ), each obtained by setting one different coefficient to \(c_k = M_{\mathrm{R}}- 1\) and all the others to \(c_k = -1\). The mode approximations are performed by finding the parameters of the Gaussian distribution \(r_k\) such that it approximates the target density for the points \(\mathcal{C}_k\) in the corresponding mode while keeping the ratio \(r_k(\boldsymbol{x})/p(\boldsymbol{x})\) from exploding for the points \(\mathcal{C}_0\) in the remainder of the region.

Specifically, we solved an optimization problem to minimize \(\delta_k\) that satisfies \[\begin{align} \big|{\log r_k(\boldsymbol{x}) - \log p(\boldsymbol{x})}\big| &\le \delta_k \end{align}\] for all \(\boldsymbol{x}\in\mathcal{C}_k\) and \[\begin{align} \log r_k(\boldsymbol{x}) - \log p(\boldsymbol{x}) &\le \delta_k \end{align}\] for all \(\boldsymbol{x}\in\mathcal{C}_0\). The point sets \(\mathcal{C}_k\) are obtained by clustering \(\mathcal{X}_1\) into \(M_{\mathrm{R}}\) clusters, and \(\mathcal{C}_0\) is obtained by MCMC for the broader distribution whose density is proportional to \((p(\boldsymbol{x}))^{1/10}\). The optimization problem can be expressed and solved in the form of linear programming.

Figs. 2 () and 2 () show these two classes of zero-mean feature functions. The plots show that the functions based on Stein’s identity have both positive and negative values in each mode, while the ratio-based functions have consistent signs within each mode. Using the target function defined as \(f(x_1, x_2)=x_1\), we evaluated how the variance of the residual \(\varepsilon(\boldsymbol{x})\) of Eq. (1 ) is reduced by the proposed method. Fig. 3 shows the predicted values using three combinations of the function classes. When the ratio-based functions are combined with the functions based on Stein’s identity, the linear model approximates the target function well, so that the residual variance is effectively reduced.

We plot the residual variance against the number of feature functions used for the linear model in Fig. 4, varying the parameters as \(\sigma \in\{0.1, 0.2, 0.5, 1.0\}\), \(M_{\mathrm{S}}= \{4, 8, 12, 25, 50, 100\}\), and \(M_{\mathrm{R}}\in \{4, 8\}\). Regardless of these parameter choices, the combination of the two classes of functions leads to a quantitatively better linear approximation that will achieve a smaller variance for estimating the expectation \(I\).

6 Discussion↩︎

In this paper, we proposed a ratio-based construction of zero-mean feature functions for variance reduction in MC estimation of expectations. Since our numerical demonstration is limited to a synthetic two-dimensional case, how the performance is retained in more realistic cases in higher-dimensional spaces remains a topic for future work. In addition, the procedure for choosing the reference mixture \(R\) is not limited to the one used in the numerical example, and there is room for improvement.

As mentioned in the earlier section, the zero-mean feature functions can be used to design the herding dynamics [7]–[9] for improving sample generation for the MC method. Applying the proposed method to this design and evaluating its performance is also an interesting direction for future exploration.

Acknowledgments↩︎

This work was supported by JSPS KAKENHI Grant Number 25K21300, JST ALCA-Next Grant No. JPMJAN23F2, JST Moonshot R&D Grant No. JPMJMS2021, JST CREST Grant No. JPMJCR25R1, and a project, JPNP14004, commissioned by the New Energy and Industrial Technology Development Organization (NEDO).

References↩︎