{'fraction_non_alphanumeric': 0.04829329962073325, 'fraction_numerical': 0.010071639275179099, 'mean_word_length': 4.721070395371263, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Neural Control Variates for Variance Reduction

In statistics and machine learning, approximation of an intractable integration is often achieved by using the unbiased Monte Carlo estimator, but the variances of the estimation are generally high in many applications. Control variates approaches are well-known to reduce the variance of the estimation. These control variates are typically constructed by employing predefined parametric functions or polynomials, determined by using those samples drawn from the relevant distributions. Instead, we propose to construct those control variates by learning neural networks to handle the cases when test functions are complex. In many applications, obtaining a large number of samples for Monte Carlo estimation is expensive, which may result in overfitting when training a neural network. We thus further propose to employ auxiliary random variables induced by the original ones to extend data samples for training the neural networks. We apply the proposed control variates with augmented variables to thermodynamic integration and reinforcement learning. Experimental results demonstrate that our method can achieve significant variance reduction compared with other alternatives.


Introduction
Most of modern machine learning and statistical approaches focus on modelling complex data, where manipulating high-dimensional and multi-modal probability distributions is of great importance for model inference and learning. Under this circumstance, evaluating the expectation of certain function f (θ) with respect to a probability distribution p(θ) is ubiquitous, µ = E θ∼p(θ) [f (θ)] = f (θ)p(θ)dθ, where the random variable of interest θ ∈ R D is typically high-dimensional.
However, in complex models, the integration is often analytically intractable. This drives the development of sophisticated Monte Carlo methods to facilitate efficient computation [15]. The Monte Carlo method is naturally employed to approximate the expectation, i.e., µ ≈ 1 where {θ i } n i=1 are samples drawn from the distribution p(θ). According to the central limit theorem, this estimator converges to µ at the rate O(1/ √ n). For high-dimensional and complex models, when p(θ) is difficult to sample from [11] or the test function f is expensive to evaluate [5], a "large-n" estimation is prohibited computationally. This directly leads to a high-variance estimator. Therefore, with a limited number of samples, how to reduce the variance of Monte Carlo estimations emerges as an essential issue for its practical use.
Along this line, various variance reduction methods have been introduced in the literature of statistics and numerical analysis. One category aims to develop appropriate samplers for variance reduction, including importance sampling and its variants [2], stratified sampling techniques [16], multi-level Monte Carlo [4] and other sophisticated methods based on Markov chain Monte Carlo (MCMC) [15]. Another category of variance reduction methods is called control variates [1,10,14,13,9]. These methods take advantage of random variables with known expectation values, which are negatively correlated with the test function under consideration. Control variates techniques can fully employ the available samples to reduce the variance, which is popular due to its efficiency and effectiveness.
However, existing control variates approaches have limitations. Firsly, most existing methods use a linear or quadratic form to represent the control function [10,14]. Although these control functions have closed forms, the representation power of them is very limited particularly when the test function f (θ) is complex and non-linear. Control functionals were proposed recently to tackle this problem [13]. However, these estimators may significantly suffer from a curse of dimensionality [12]. Secondly, when the available samples are scarce, optimizing the control variates only based on a small number of samples might overfit, which means that it is difficult to generalize on the samples obtained later. These restrictions limit their practical performance.
In order to overcome these issues, in this paper we propose to employ neural networks to represent the control functions, utilising the capability of a neural network to represent a complex test function. This method is named as "Neural Control Variates" (NCV). To accommodate the overfitting issue when available training sample size is small, we augment an auxiliary variable to the original one such that the extra conditional information could be used, abbreviated as NCV with augmented variables (NCVA). Our methods are particularly suitable for the cases when the sample space is high-dimensional or the samples from p(θ) is hard to obtain. We demonstrate the effectiveness of our methods in thermodynamic integration and policy optimization in reinforcement learning. We show that the new control variates methods achieved the best performance comparing to the state-of-the-art methods in literatures.
(1) Thus the desired expectation can be replaced by that of the auxiliary function It is possible to obtain a variance-reduced Monte Carlo estimator by selecting or optimizing g(θ) so that the variance Intuitively, variance reduction can be achieved when g(θ) is negatively correlated with f (θ) under p(θ), since much of the randomness "cancels out" in the auxiliary functionf (θ).
The selection of an appropriate form of control function g(θ) is crucial for the performance of variance reduction. A tractable class of so called zero-variance control variates was proposed in [10]. Those control variates are expressed as a function of the gradient of the log-density, ∇ θ log p(θ), i.e. the score function s(θ). Concretely, it has the following form where the gradient operator ∇ θ = [∂/∂θ 1 , . . . , ∂/∂θ D ] T , the Laplace operator ∆ θ = D i=1 ∂ 2 /∂θ 2 i , and "·" denotes the inner product. The function Q(θ) is often referred to as the trial function. The target is now to find a trail function so that g(θ) and f (θ) are negatively correlated. This could thus reduce the variance of the Monte Carlo estimation. As the trail function could be arbitrary under those mild conditions given in [10], a prametric function could be used for Q(θ). For example, when Q(θ) = a T θ, which is a first degree polynomial function, the auxiliary function becomes f (θ) = f (θ) + a T s(θ) as was proposed in [10]. The optimal choice of the parameter a that minimizes the variance off . Obviously, the representation power of these polynomials is limited, and therefore control functionals have been proposed recently where the trail function is stochastic. For example, the trail function could be a kernel function [13]. In order for using these control variates, we firstly estimate these required parameters in the trail function by using some training samples {θ i } n i=1 . Then the learned control variates can be used for test samples.
However, there are several drawbacks of the current zero-variance techniques: 1) Dilemma between effectiveness and efficiency. Although increasing the order of polynomial could potentially increase the representation power and reduce more variance, the number the parameters needs to be learned would grow exponentially. As pointed out by [10], when quadratic polynomials are used, Q(θ) = a T θ+θ T Bx/2, the number of parameters will be D(D+3)/2. Thus, finding the optimal coefficients requires dealing with Σ ss which is a matrix of dimension of order D 2 . Similar issue occurs when employing the control functionals. This makes the use of high order polynomials computationally expensive when faced with high-dimensional sampling spaces. 2) Poor generalization with small sample size. With small sizes of training samples and complex p(θ), the learned control variates could potentially overfit the training samples, i.e. generalize poorly over new samples. This is because a small size of training sample might be insufficient for representing the full distributional information of p(θ).
These limitations motivate our development of neural control variates and its extension with augmented variables, which will be elaborated below.

Neural Control Variates and Its Extension
Firstly, we focus on alleviating the dilemma between effectiveness and efficiency on designing control variates in high-dimensional sample space. To this end, we propose neural control variates (NCV) where the trail function is represented by a neural network. Equipped with neural networks, their excellent capability of representing complex functions and overcoming the curse of dimensionality can be fully employed in high-dimensional scenarios [7].
Instead of relying on the control variates (3) introduced in [10], we use the following Stein control variates based on Stein identity [19,13], where Φ(θ) is the trial function. Note that in order for E[g(θ)] = 0, we admit those assumptions made in [10,13] which could be easily satisfied under some mild conditions. Compared with Eq (3), Stein control variates is preferred due to its computational advantages since evaluating the second order derivatives of the trial function is avoided. Note that when the trial function Q(θ) is a linear or quadratic polynomial, the Stein trial function Φ(θ) is constant or linear, respectively.
We now represent the trial function Φ(θ) by a neural network Φ(θ; w) parameterized by the weights w. The control function becomes g(θ; w) = ∇ θ · Φ(θ; w) + Φ(θ; w) · ∇ θ log p(θ). In order for variance reduction, we solve the following optimization problem which does not have a closed-form in general. Instead, the following optimization problem will be solved where {θ i } n i=1 are samples drawn from p(θ). Consequently, standard back-propagation techniques and stochastic gradient descent (SGD) can be adopted to obtain the optimal network parameters.

Neural Control Variates with Augmented Variables
When the distribution p(θ) is high-dimensional and multi-modal, it would be very expensive to draw samples, which are required for training those control variates. This typically produces a rather small number of samples which is not sufficient for training control variates. When using control variates methods, the learned parameters for control variates can easily overfit over the training samples, and thus could not generalize well on new samples drawn from p(θ). To deal with this issue, we augment an extra variable δ ∈ R D δ to θ to establish a joint distribution p(θ, δ). This neural control variates with augmented variable (NCVA) method could be helpful to remedy overfitting. Obtain a mini-batch of samples with with size n B n from {θ} n i=1 and evaluate their scores Evaluate the stochastic gradient of J NCVA(w) based on samples (θ i , δ ij ), i = 1, . . . , n B ; j = 1, . . . , M ; Update w using SGD. 7: until Convergence 8: Return the learned parameter w * . 9: 10: //Testing Concretely, our proposal is inspired by the following observation. Given the augmented control function g(θ, δ) satisfying E (θ,δ)∼p(θ,δ) [g(θ, δ)] = 0, we have With the joint distribution p(θ, δ) = p(θ)q(δ|θ), the equation above can be expanded further as follows where the augmented score function can be written as The equations (8)-(10) facilitate expanding the training sample size by adding augmented variable samples from the conditional distribution q(δ|θ). This trick could essentially help to reduce overfitting. Intuitively, the influence of augmented variable δ for training Φ are delivered by the change of the score function presented in equation (10).
The selection of the conditional distribution q(δ|θ) depends on the application of interest and how well the samples approximate the true density p(θ). For simplicity, we obtain the samples of δ by firstly sampling ν ∼ Gamma(α, 1), and then δ = νθ, where the shape parameter ν can be tuned by user. We use this simple scheme for all the applications in Section 4. As for how to choose the optimal form of q(δ|θ) for best generalization, we leave it as our future work.
After parameterizing the augmented trial function by a neural network Φ(θ, δ; w), the following objective function will be minimized, where the "effective control variates" E δ∼q(δ|θi) [g(θ i , δ; w)] can be estimated with M Monte Carlo samples by an online fashion. We present the training and test procedures of NCVA in Algorithm 1.

Experiments
To evaluate our proposed methods, i.e., NCV and NCVA, we apply them to a synthetic problem and two real scenarios, which are thermodynamic integration for Bayesian model evidence evaluation, and policy optimization in reinforcement learning. For comparison purposes, control functionals (CF) [13] and polynomial control variates [10,14] are also applied to these problems. The performance of the trained control variates are measured by the variance reduction ratio on the test data set, i.e., . We used fully connected neural networks to represent the trial function in all the experiments. We found that for the experiments presented in the following, a medium-sized network is empirically sufficient to achieve good performance. More details on network architectures are provided in the Appendix.

Synthetic Problem
To illustrate the advantage of NCV on dealing with high-dimensional problems over other methods, we consider to approximate the expectation of f (θ) = sin(π/D D i=1 θ i ) where θ ∈ R D which is a mixture of Gaussians, i.e., p(θ) = 0.5N (−1, I) + 0.5N (1, I). Figure 1 shows the variance reduction ratio on test data (N = 50) with respect to varying the number of training samples and the dimensions. In both cases, we can observe that NCV outperforms linear, quadratic and control functionals control variates. Particularly, when increasing the dimensions of θ, NCV can still achieve much lower variance reduction ratio compared with CF.

Thermodynamic Integral for Bayesian Model Evidence Evaluation
In Bayesian analysis, data y is assumed to have been generated under a collection of putative models, {M i }. To compare these candidate models, the Bayesian model evidence is constructed as p(y|M i ) = p(y|θ, M i )p(θ|M i )dθ where θ are the parameters associated with model M i . Unfortunately, for most of the models of interest, this integral is unavailable in closed form. Thus many techniques were proposed to approximate the model evidence. Thermodynamic integration (TI) [3] is among the most promising approach to estimate the evidence. This approach is derived from the standard thermodynamic identity, log p(y) = 1 0 E p(θ|y,t) [log p(y|θ)]dt where p(θ|y, t) ∝ p(y|θ) t p(θ), t ∈ [0, 1] is called power posterior. Note that we have dropped the model indicator M i for simplicity. Here t is known as an inverse temperature parameter. In many cases, the posterior expectaion E p(θ|y,t) log(p(y|θ)) can't be analytically computed, thus the Monte Carlo integration is applied. However, Monte Carlo integration often suffer large variance when sample size is not large enough.
In [14], the zero-variance control variates (3) were used to reduce the variance for TI, so that where {θ t i } N i=1 are drawn from the posterior p(θ|y, t). In [14], the trial function Q t (θ) was assumed as a linear or quadratic function, which corresponds to a constant or linear function for Φ(θ) in Stein control variates 2 . These methods achieved excellent performance for simple models [14]. However, they are struggling in some scenarios, for example, a negative example which is Goodwin Oscillator is given in [14]. Note that Goodwin Oscillator is a nonlinear dynamical system, where the form of f (·) is provided in the Appendix. Assuming within only a subset of time points {s i } N i=0 , the solution of (13), i.e. x(s i , θ), is observed under Gaussian noise ε(s) ∼ N (0, σ 2 I), where σ 2 denotes the variance of the noise. That means the observation y(s i ) = x(s i ) + ε(s i ). Then we have the likelihood p(y|θ, The expectation of the log likelihood under the power posterior, i.e., E p(θ|y,t) log p(y|θ), needs to be evaluated. In [14], the authors demonstrated the failure of polynomial-type of control function since the log-likelihood surface is highly multi-modal and there is much weaker canonical correlation between the scores and the log posterior.
In practice, sampling from Goodwin Oscillator is difficult and computationally expensive since simulating the underlying ordinary differential equation is extremely time-consuming. This directly leads to the situation that the available training samples for control variates are not sufficient. We show in the following that the proposed NCVA can be employed to remedy this issue. To illustrate the benefits of NCVA, we compared it with other methods with various sizes of training samples and temperatures. For comparison purposes we evaluated the variance reduction ratios on both training and test sets (500 samples for test). The experiment settings are the same as those in [14].  Figure 2 shows the experimental results when applying different types of control variates. It can be easily observed that the linear and quadratic methods could hardly reduce the variance of the Goodwin Oscillator model on test set. Both neural control variates and control functionals perform well when the training size is larger than 500. However, when the sample size is small, CF and NCV suffer from overfitting. Comparing to all other methods, NCVA successfully achieved the lowest variance reduction ratio when the sample size was small.

Variance Reduction for Policy Optimization
In this application, the proposed NCVA is employed to reduce the variance of policy gradient in reinforcement learning. Policy gradient methods have achieved great success in reinforcement learning problems [17,22,6,18]. It assumes the action a is selected from a conditional probability distribution π(a|s; θ) given a state s and parameterized by θ, typically named as policy. The goal of the policy gradient method is to obtain the optimal parameters by maximizing the expected cumulative reward, J(θ) = E s,a∼π(s)π(a|s) [r(s, a)].
According to [20], the gradient of J(θ) can be written as where Q π (s, a) denotes the expected reward under the policy π starting from state s and action a. When training a policy, a Monte Carlo estimation of the policy gradient (15) is computed using the trajectory {(s t , a t , r t )} n t=1∇ whereQ π (s t , a t ) is the empirical estimate of Q π (s, a), and the trajectory {(s t , a t , r t )} n t=1 is obtained by simulating the environment under the current policy π. However the Monte Carlo policy gradient estimator often suffers a large-variance issue when the number of trajectories is small, which makes the training process inefficient.
Neural Stein control variates were adopted to reduce the variance of the policy gradients estimator in [8], which can be described as follows. Given a policy π(a|s), Stein's identity with respect to the action a is E π(a|s) [∇ a log π(a|s)φ(s, a) + ∇ a φ(s, a)] = 0, ∀s.
where φ(s, a) is the trial function of the Stein control variate. However, (17) can't be directly applied in policy gradient form (15), since the gradient in (15) is taken with respect to the policy parameters θ, while the gradient in Stein identity (17) is taken with respect to the action a. To connect ∇ a log π(a|s) and ∇ θ log π(a|s), the policy π(a|s) is reparameterized to adapt the Stein identity to the policy gradients framework, i.e. assuming the action a is generated by a mapping a = f (s, ξ), where ξ is a random noise drawn from some distribution independent of θ and s. Then we have Combining (18) and (15), the policy gradient with NCV could be derived: Now we introduce the augmented variable method to adapt the reparametrized (18): where the augmented variable b is generated by b = g(a, ζ), the noise ζ drawn from a pre-defined distribution independently of s, θ, ξ. Thus we have ∇ θfθ (s, ξ, ζ) = (∇ θ f (s, ξ), ∇ θ g(f θ (s, ξ), ζ)) T . Besides, ∇ θ log π(a, b|s) = ∇ θ log π(a|s) + ∇ θ log π(b|a, s) = ∇ θ log π(a|s) + ∇ θ log π(ζ) = ∇ θ log π(a|s).
Thus the augmented form of the gradient estimator could be rewritten as: where the control variate is decomposed as φ(s, a, b) =V (s) + ψ(s, a, b), the termV (s) is a parametric function approximation of the value function, andÂ(s t , a t ) =Q π (s t , a t ) −V (s t ) is an estimator of the advantage function [20,22].
In [9], the authors evaluated their reparametrized Stein CV method combining proximal policy optimization (PPO, [18] in continuous control environments from the OpenAI Gym benchmark using the MuJoCo physics simulator. We compare the performance of their method and the augmented version on three different MuJoCo environments, Halfcheetah-v1, Walker2d-v1 and Hopper-v1. Note that in policy gradient for RL applications, the batch size N is typically around several thousands, e.g. for continuous control problems, N = 5000 or 10 4 . This number of sample size is computationally expensive for kernel-based methods, such as control functionals, which are nearly intractable for RL problems. Therefore, we do not consider CFs in this application. Figure 3 shows that when the sample size is relatively small (N = 5000), the reparametrized Stein CVs are not stable and may meet "barrier" to train the policy and obtain higher reward efficiently.
Favorably, the usage of NCVA in policy gradient is efficient and stable. For the large sample size (N = 10 4 ), though [9] obtains good performance as stated, our NCVA still outperforms.

Conclusion & Future work
We have developed the zero-variance control variate methods for Monte Carlo variance reduction by incorporating neural networks into Stein control variates. The neural control variates can deal with complex test function with high-dimensional distribution. Besides, we proposed the augmented variables approach to achieve better generalization for small sample size cases. We demonstrated their effectiveness on two challenging Monte Carlo integration tasks. In the future, we will mainly focus on exploring the theoretical justification on the effectiveness of the augmented variables method. A wiser scheme for constructing the auxiliary random variable also need to be developed.

A Experimental Settings
A.1 Synthetic Problem The following kernel function used in [13] is employed in control functionals in our experiments where α1 = 0.1 and α2 = 1 which were selected by using cross-validation. The trial function Φ in NCV is a neural network, the architecture of which has two fully connected layers. Each layer has 40 units, and the activation function is RELU. Adam was used to train the network with learning rate 0.008.
The control functionals control variates used the kernel with parameters (23) with α1 = 0.01, α2 = 5 which were chosen by using cross-validation. The neural network used for the trial functions in NCV and NCVA has two fully connected layers, each of which has 40 units. Adam is used to train the network with learning rate 0.05. For NCVA, the augmented variable was generated by δ = ξθ, where ξ ∼ Γ(10, 1). All samples were augmented 100 times.

A.3 Policy Gradients with Stein Control Variates
We reimplement the proxy policy optimization procedure with stein control varaites ( [9]) in three mujoco tasks( [21]). The neural network for the trial function in NCV and NCVA has three MLP layers, each of which has 100 units. Adam is used to train the network with learning rate 0.0005. The policy was assumed a Gaussian a ∼ N (µ(s), Σ(s)) where µ(s) and diagonal Σ(s) were both constructed as a MLP with three layers. In NCVA, all samples are augmened 100 times. The augmented variable is generated by δ = ξθ, where ξ ∼ Γ(5, 1).