“…In the Bayesian context, Assaraf & Caffarel (1999); Mira et al (2013) and Oates et al (2017) addressed this challenge by using f n = c n + Lg n where c n ∈ R, g n is a user-chosen parametric or nonparametric function and L is an operator, for example the Langevin Stein operator (Stein, 1972;Gorham & Mackey, 2015), that depends on p through its gradient and satisfies (Lg n )(x)p(x)dx = 0 under regularity conditions (see Lemma 1). Convergence of I CV (f ) to I(f ) has been studied under (strong) regularity conditions and, in particular (i) if g n is chosen parametrically, then in general lim inf σ(f − f n ) 2 > 0 so that, even if asymptotic variance is reduced, convergence rates are unaffected; (ii) if g n is chosen in an appropriate non-parametric manner then lim sup σ(f − f n ) 2 = 0 and a smaller number O( −2+δ ), 0 < δ < 2, of calls to f , p and its gradient are required to achieve an approximation error of O P ( ) for the integral (see Oates et al, 2019;Mijatović & Vogrinc, 2018;Barp et al, 2018;Belomestny et al, 2017Belomestny et al, , 2019Belomestny et al, , 2020. In the parametric case Lg n is called a control variate while in the non-parametric case it is called a control functional.…”