Control Functionals for Quasi-Monte Carlo Integration

0.9 MB·English

by  Chris. J. Oates

Control Functionals for Quasi-Monte Carlo Integration
Chris. J. Oates Mark Girolami
University of Technology Sydney and University of Warwick and
The Alan Turing Institute for Data Science
6 Abstract points are ‘well-spaced’ in a precise sense defined be-
1
low. This contrasts with the Monte Carlo (MC) ap-
0
2 Quasi-Monte Carlo (QMC) methods are be- proach whereby the design points are sampled inde-
pendently from a uniform distribution over X. MC
r ing adopted in statistical applications due to
p integrationachievesarootmeansquareerror(RMSE)
theincreasinglychallengingnatureofnumer-
A convergence rate of O(N−1/2) whereas QMC integra-
ical integrals that are now routinely encoun-
tion can in principle achieve a rate O(N−α/d) on spe-
1 tered. For integrands with d-dimensions and
cific geometric sequences {N }∞ [24]. It is known
derivatives of order α, an optimal QMC rule n n=1

] converges at a best-possible rate O(N−α/d). that this rate is best-possible [19] and explicit algo-
O
rithms to generate design points that attain this rate
However, in applications the value of α can
C arenowavailableformany(butnotall)valuesofα[6].
beunknownand/orarate-optimalQMCrule
. Challenging integration problems are common in con-
t can be unavailable. Standard practice is to
a temporary statistics, for example when computing ex-
st employ αL-optimal QMC where the lower pectations,marginalprobabilitydensitiesornormalis-
bound α ≤ α is known, but in general this
[ L ingconstants,andQMCmethodsarethereforegaining
does not exploit the full power of QMC. One
7 importance in statistical applications [12, 17, 35].
solution is to trade-off numerical integration
v
with functional approximation. This strat- Contrarytotheabovetheoreticalconsiderations,rate-
9
7 egy is explored herein and shown to be well- optimal QMC is often not employed in practice. This
3 suited to modern statistical computation. A ismainlyduetothreereasons;either(R1)thesmooth-
3 challenging application to robotic arm data nessparameterαisunknown,(R2)theredoesnotcur-
0
demonstrates a substantial variance reduc- rently exist an explicit QMC rule that is rate-optimal
.
1 tion in predictions for mechanical torques. for functions of smoothness α, or (R3) it is simply
0 more convenient to employ a basic QMC rule based
5
onaweakersmoothnessassumptionα <α,asimple-
1 L
1 Introduction mented in standard software. In each situation there
:
v is a gap between theory and practice that, as we show
i
X Consider a Lebesgue-integrable test function f :X → in this paper, can be bridged using functional approx-
r R defined on a bounded measurable subspace X ⊆Rd imation.
a (d∈N)withsquareintegrablederivativesoforderα>
Previous work on variance reduction techniques for
0ineachvariable. Ourfocusisnumericalcomputation
QMC includes [1], who considered modified impor-
(cid:82)
of the integral I[f] := f(x)dx. The Quasi-Monte
X tance sampling strategies, and [14], who considered
Carlo (QMC) approach is based on an approximation
constructing control variates for QMC. Neither ap-
proach improved the asymptotic error rate, though in
N
Q[f;x1:N]:= 1 (cid:88)f(xn) some cases the QMC error was reduced by a constant
N factor. Interestingly, [14] reports some quite negative
n=1
resultsforcontrolvariatestrategiesinthissetting,be-
where the (possibly random) design points x1:N = cause the objective being minimised by QMC is not
{x1,...,xN} ⊂ X have low discrepancy; that is, the equivalent to the MC variance that is minimised by
controlvariates. [33]demonstratesvariancereductions
Appearing in Proceedings of the 19th International Con- in QMC are possible using additive approximations,
ference on Artificial Intelligence and Statistics (AISTATS) though again the asymptotics were unchanged.
2016, Cadiz, Spain. JMLR: W&CP volume 41. Copyright
2016 by the authors. This paper studies a general approach to variance re-
Control Functionals for Quasi-Monte Carlo Integration
duction for QMC rules, building on kernel methods (MMD), as studied extensively in the kernel methods
and recent work in the Monte Carlo setting due to literature [4, 30].
[22, 31]. The mathematics that underpins our work
Quadrature is naturally studied in RKHS because
comes from the functional approximation literature.
thereexistsaclosed-formexpressionfortheworstcase
Thistakestheformofa‘controlfunctional’ψ :X →R
error in terms of the kernel K, which facilitates the
thatsatisfies(i)ψ integratestozero,(ii)f−ψ ismore
principled selection of design points [8]:
amenabletoQMCmethodsthanf, inaprecisesense.
(cid:90) (cid:90)
Thegeneralapproachthatweexploreistoreplacethe
e (x1:N)2 = K(x,y)dxdy
integrandf byf−ψandtargettheQMCobjectivedi- H
[0,1]d
rectly. This can lead to accelerated asymptotics. The
maincontributionofthispaperistoexplorethisstrat- − 2 (cid:88)N (cid:90) K(xn,y)dy
egy in the settings (R1-3) above. Theoretical analysis N n=1 [0,1]d
of convergence rates is provided, along with empirical
N
1 (cid:88)
results and a challenging application to robotics. We + K(xn,xm) (3)
begin by presenting some background on QMC the- N2
m,n=1
ory below, before describing the methodology in more
The mainstream QMC literature supposes H is a
detail.
Sobolev space of known order α (defined below). In
thissetting,O(N−α/d)isthebest-possiblerateforthe
2 Background
worst case error when x1:N are chosen deterministi-
cally and O(N−α/d−1/2) is the best-possible RMSE
QMC is naturally studied in reproducing kernel when x1:N are allowed to be random [19]. We will re-
Hilbert spaces (RKHS; [8]). Below we draw connec- fer to QMC rules that achieve these optimal rates as
tions with kernel methods, that are themselves natu- ‘α-QMC rules’.
rally studied in RKHS.
This paper focuses on improving performance in the
Notation. We work in a Hilbert space H, consist- situation where a (sub-optimal) α -QMC rule is used
L
ing of measurable functions f : X → R. For simplic- tointegrateatestfunctionofsmoothnessα>α . For
L
ityofpresentationweassumeH includestheconstant reasons(R1-3),thisscenarioiscommonlyencountered
functions. We follow the mainstream QMC literature instatisticalapplications. IncontrasttoQMC[8](and
by taking X = [0,1]d, equipped with the Euclidean kernel methods that aim to minimise the MMD [2]),
norm (cid:107)x(cid:107) := ((cid:80)d x2)1/2. Denote the scalar prod- the rate constant (cid:107)f(cid:107) is the primary target of our
i=1 i H
uct and norm on H by (cid:104)·,·(cid:105) and (cid:107)·(cid:107) respectively. methodology below.
H H
Suppose further that H is a RKHS with kernel K :
[0,1]d×[0,1]d →R;thatis,K satisfies(i)K(·,x)∈H 3 Methodology
for all x ∈ [0,1]d and (ii) f(x) = (cid:104)f,K(·,x)(cid:105) for all
H
f ∈ H and all x ∈ [0,1]d. K is assumed to be non-
Control Functionals for QMC.Theapproachthat
trivial, i.e. K (cid:54)=0.
we pursue in this paper aims to construct a Lebesgue-
Quadrature Error Analysis. The quadrature integrable functional ψ :[0,1]d →R that satisfies
methods that we focus on aim to minimise the ‘worst
I[ψ]=0. (4)
case’ integration error which, for design points x1:N
and Hilbert space H, is defined to be When x has the interpretation of a random variable,
ψ(x) is classically known as a ‘control variate’ [14].
eH(x1:N):= sup (cid:12)(cid:12)Q[f;x1:N]−I[f](cid:12)(cid:12) (1) When ψ itself is estimated, we follow [22] and refer to
(cid:107)f(cid:107)H≤1 theentiremappingψ asa‘controlfunctional’(CF).In
the CF approach to estimation, the test function f is
where the supremum is taken over all test functions
replaced by f −ψ; it is hoped that the latter is more
f belonging to the unit ball in H. It follows from
amenabletonumericalintegration. ClearlyI[f−ψ]=
linearity that, for any function f ∈H, the integration
I[f]. In this paper we construct a CF ψ based on
error obeys N
a tractable approximation f to f. (The dependence
N
(cid:12)(cid:12)Q[f;x1:N]−I[f](cid:12)(cid:12)≤eH(x1:N)(cid:107)f(cid:107)H. (2) on N will be explained below.) It is required that the
integral I[f ] is available in closed-form. We then set
N
The worst case error e (x1:N) is the usual target
H ψ (x)=f (x)−I[f ] (5)
of QMC innovation, with x1:N chosen to (approxi- N N N
mately, asymptotically) minimise e (x1:N) [8]. Note so that ψ satisfies Eqn. 4. For this to make sense
H N
that Eqn. 1 is also the ‘maximum mean discrepancy’ mathematically, it must be the case that f ∈H and
N
Chris. J. Oates, Mark Girolami
this informs our method of approximation (the con- ‘p-Sobolev space of order k’ to be
stant function with value I[f ] belongs to H by as-
N Wk,p := {f :[0,1]d →R|Daf exists and
sumption). Intuitively, a good CF ψ will provide a
N
closeapproximationtofluctuationsofthetestfunction Daf ∈Lp([0,1]d),∀a∈Nd0 with |a|≤k}.
f, so that the functional difference f − ψ become
N Here Daf denotes the weak (or ‘distributional’)
increasingly ‘flat’ and thus more amenable to QMC
derivative of f; the reader is referred to the above
methods. Moreprecisely,motivatedbyEqn. 2weaim
reference for details. Clearly Wk,p is a vector space
toconstructaCFsuchthat(cid:107)f−ψ (cid:107) <(cid:107)f(cid:107) . This
N H H over R when addition and (scalar) multiplication are
connection with functional approximation offers the
definedpoint-wise. Forthespecialcasep=2weequip
possibility to leverage kernel methods for these prob-
Wk,2 with the inner product
lems, see e.g. [27, 31].
(cid:88)
(cid:104)f,g(cid:105) := I[DafDag]
ControlFunctionalErrorAnalysis. Considerpar- k
titioning x1:N into two sets u1:M and vM+1:N where a∈Nd,|a|≤k
0
1 < M < N and M/N → c ∈ (0,1) as N → ∞.
and denote this inner-product space Hk :=
The first set u1:M, possibly non-random, will be used
(Wk,2,(cid:104)·,·(cid:105) ). Defined in this way, Hk is a Hilbert
k
in a preliminary step to construct an approximation
space of functions whose (weak) derivatives exist up
f (·;u1:M) to f. Then the second set vM+1:N, possi-
M to order k. Moreover Hk can be made into a RKHS
bly random, is used to evaluate the ‘CF estimator’
with an appropriate choice of kernel (see below). Our
results below apply also to Sobolev spaces with non-
E[f;u1:M,vM+1:N] := Q[f −ψ (·;u1:M);vM+1:N]
N integer k, but this construction is more technical and
= Q[f −fM(·;u1:M);vM+1:N] we refer the reader to [27] for details.
+I[f (·;u1:M)]. (6)
M Approximation in Sobolev Spaces. Our assump-
tionsarenaturallystatedusingSobolevspaces: Given
We remark that if the points vn are ran- two Hilbert spaces H, H(cid:48), defined on the same ele-
dom and marginally distributed as U([0,1]d) then
ment set, with norms (cid:107)·(cid:107) , (cid:107)·(cid:107) , we say that H
H H(cid:48)
E[f;u1:M,vM+1:N] will be an unbiased estimator for and H(cid:48) are ‘norm-equivalent’, written H ≡H(cid:48), when-
I[f]. Error analysis for the CF estimator is based on
ever there exist positive constants c , c such that
1 2
the following:
c (cid:107)f(cid:107) ≤(cid:107)f(cid:107) ≤c (cid:107)f(cid:107) for all f ∈H.
1 H H(cid:48) 2 H
Theorem 1. Given f,f ∈H, we have
M Assumption 1: H ≡HαL where αL >d/2.
|E[f;u1:M,vM+1:N]−I[f]| Assumption 2: f ∈Hα where α≥αL.
≤e (vM+1:N)(cid:107)f −f (·;u1:M)(cid:107) . (7) Assumption 1 is a technical requirement to ensure the
H M H
spaceH (whereQMCisperformed)admitsconsistent
functional approximation. Assumption 2 ensures that
Proof. Since f,fM ∈ H we have that f −fM ∈ H. the test function f is ‘smooth enough’ for αL-QMC
The result then follows by applying the fundamental methodstoconvergeattheα -rate. Thisfollowsfrom
L
inequality from Eqn. 2 to the function f −fM and the fact that Sobolev spaces are nested, so that f ∈
using linearity of the integral operator I. Hα =⇒ f ∈HαL.
For consistent approximation of f it is necessary to
Thus the CF methodology produces an estimator baseourapproximationf inaspaceH offunctions
M ∗
E[f;u1:M,vM+1:N]thathasasymptoticallyzeroerror that are ‘at least as smooth’ as f:
relative to standard QMC estimators, providing that
itispossibletoconstructanapproximationf tof in Assumption 3: H∗ ≡HαU where αU ≥α.
M
such a way that (cid:107)f −f (·;u1:M)(cid:107) →0 as M →∞. It follows again from the nested property that f ∈
M H M
Thenextsectionsestablishconvergenceratesforfunc- HαL and thus the functional difference f −fM exists
tional approximation using kernel methods. in HαL. The Sobolev spaces H∗ can be characterised
as RKHS via an appropriate reproducing kernel K ,
Sobolev Spaces. To achieve consistent approxima- ∗
such as the well-known Mat´ern kernel.
tion (cid:107)f −f (cid:107) → 0 it is necessary to impose regu-
M H
larity conditions on H. Sobolev spaces are a general Finallyanapproximationf tof isconstructedbased
M
setting in which to formulate such regularity assump- on the points u1:M as follows:
tions; our main reference here is [27]. Firstly suppose
M
that k ∈ N0, k > d/2 and 1 ≤ p < ∞. For a multi- f (x;u1:M):= (cid:88)β K (x,un) (8)
index a∈Nd we write |a|=a +···+a . Define the M n ∗
0 1 d n=1
Control Functionals for Quasi-Monte Carlo Integration
where the weights β ∈ R are defined as the solution smoothnessisassumedsothat, asaruleofthumb, we
n
to the linear system of interpolation equations should try to select α as close as possible to α. Our
U
main result is stated below:
f (un;u1:M)=f(un), n=1,...,M. (9)
M
Corollary 1. When u1:M is quasi-uniform, CFs ac-
It is well-known that Eqn. 8 is the unique minimiser celerate αL-QMC by a factor O(N−(α−αL)/d).
of the H -norm under all functions in H that satisfy
∗ ∗
the linear system in Eqn. 9 [27]. In practice it may Remark: The improvement due to CFs appears to be
benecessarytoregularisethelinearsysteminorderto mainly limited to low-dimensional integrals (d small),
facilitate inversion, but we do not go into details here, but in fact CFs can in principle be extended to high-
see e.g. [27]. dimensional integrals under additional tractability as-
sumptions, as discussed in Sec. 5.
We note that I[f ] will not have a closed-form ex-
M
pression when the Mat´ern kernel is employed and for Remark: Optimising the bound in Theorem 2 enables
this technical reason we instead employ tensor prod- us to obtain the optimal scaling
ucts of polynomial kernels (these give rise to Sobolev
M α−α
spaces of mixed dominating smoothness - full details →c∗ = L,
are provided at the end of this section). N α
Theory: Deterministic Case. We begin by con- see the Supplement for full details.
sidering the case where the design points vM+1:N are
The overall convergence rate of the CF estimator de-
chosen deterministically. Define the ‘fill distance’
pends on how the design points vM+1:N are gener-
h(u1:M):= sup min(cid:107)x−un(cid:107), ated. For this there are many QMC methodologies
x∈[0,1]d n available, each leading to different convergence rates
for the worst case error e (vM+1:N); see [7] for a
the ‘separation radius’ HαL
recent survey of some of these approaches. Of partic-
1 ular interest in statistical applications is the case of
q(u1:M):= min(cid:107)uj −uk(cid:107)
random design points which we discuss below.
2 j(cid:54)=k
Theory: Randomised Case. Modern QMC meth-
and the ‘mesh ratio’ ρ(u1:M) := h(u1:M)/q(u1:M).
ods begin with a deterministic set/sequence of design
The set u1:M is called ‘quasi-uniform’ if ρ(u1:M) → 1
points (e.g. a Halton sequence or a Sobol sequence),
as M →∞.
then apply a random transformation leading to a low
Theorem 2. Under Assumptions 1-3 the CF estima-
discrepancy set with high probability. Below we con-
tor has error bounded by
sider three types of randomisation; shifting, folding
and scrambling.
|E[f;u1:M,vM+1:N]−I[f]|
≤CeHαL(vM+1:N)h(u1:M)α−αLρ(u1:M)αU−αL(cid:107)f(cid:107)Hα Shifting: In ‘random shift’ QMC the design points
vM+1:N are translated by a common uniform random
where C >0 is a constant that depends on α, α and vector ∆ ∈ [0,1]d, so that vn (cid:55)→ vn + ∆ for each
L
α but not on f, vM+1:N and u1:M. n = M + 1,...,N. For convenience we write this
U
‘shifted’ set as vM+1:N +∆. Applying Theorem 2 to
Proof. From [27] (Theorem 7.8) we have that the ker- vM+1:N +∆ and then marginalising over ∆ ∈ [0,1]d
nel estimator in Eqn. 8 is consistent for the non- produces a RMSE bound for the CF estimator:
parametric regression problem at a rate
Corollary 2. Under Assumptions 1-3 the random
shift CF estimator has error bounded by
(cid:107)f −f (·;u1:M)(cid:107)
M HαL
≤Ch(u1:M)α−αLρ(u1:M)αU−αL(cid:107)f(cid:107)Hα (cid:113)E|E[f;u1:M,vM+1:N +∆]−I[f]|2
wwhitehreECqn.de7pceonmdsploentelys tohneαp,rαooLf,.αU. Combining this ≤CesHhαL(vM+1:N)h(u1:M)α−αLρ(u1:M)αU−αL(cid:107)f(cid:107)Hα
where
For quasi-uniform u1:M, there is no asymptotic
(cid:90)
penaltyfromemployingakernelK thatimposes‘too
∗ (esh (vM+1:N))2 := e (vM+1:N +∆)2d∆
much smoothness’ on the approximation f , with HαL HαL
M [0,1]d
ρ→1. Inthiscaseh(u1:M)=O(M−1/d)and,sinceM
and N are proportional, h(u1:M) = O(N−1/d). How- and C > 0 is a constant that does not depend on f,
ever the rate constant C will increase when too much vM+1:N or u1:M.
Chris. J. Oates, Mark Girolami
For quasi-uniform u1:M, CFs accelerate random shift smoothness’ on f in order to be able to faithfully
M
αL-QMCbyafactorO(N−(α−αL)/d)(compareagainst approximate f (Assumption 3). Moreover, K∗ should
Sec. 5.2 of [7]). be tunable to achieve a pre-specified minimum level
of smoothness. Below we make an explicit connection
Folding: A shifted and ‘folded’ QMC rule takes the
between K and the order of the associated ‘native’
form ∗
Sobolevspacethatwillallowustosatisfythisrequire-
1 (cid:88)N ment. (ii) The functions K∗(·,y) can be integrated
Q (f;z1:N +∆):= f(b(zn+∆)) analytically, so that I[f ] is available in closed form.
b N M
n=1 This second requirement leads us to consider tensor
products of Sobolev spaces, as described below.
where b is the ‘baker’s transformation’, given by
b (t)=1−|2t −1|. Thistransformationreduceserror Toconstructanalyticallyintegrablefunctionalapprox-
i i
rates; for example, for f ∈ SH2([0,1]d) (defined be- imations we consider kernels that are given by poly-
low), folding and shifting a uniform lattice z1:N leads nomials. Wendland’s compactly supported functions
to a RMSE O(N−2+(cid:15)) that is smaller than the RMSE [34] are defined via the recursion
O(N−1+(cid:15)) for a shifted lattice (p. 59 of [7]). The CF
ϕ =Ik[ϕ ],
estimator here is d,k (cid:98)d/2(cid:99)+k+1
the base function ϕ (r) = (1 − r)(cid:96) with x :=
(cid:96) + +
Eb[f;u1:M,vM+1:N +∆] max{0,x}, and the integral operator
:=I[f (·;u1:M)]+Q [f −f (·;u1:M);vM+1:N +∆]. (cid:90) ∞
M b M
I[ϕ](r)= tϕ(t)dt
For convenience we denote the shifted and folded de- r
signpointsbyb(vM+1:N+∆). ApplyingTheorem2to (r ≥0), so that
b(vM+1:N+∆)andthenmarginalisingover∆∈[0,1]d (cid:26) (1−r)(cid:96)+kp (r), r ∈[0,1]
ϕ (r)= d,k
produces: d,k 0, r >1
Corollary 3. Under Assumptions 1-3 the shifted and
where (cid:96) = (cid:98)d/2(cid:99)+k+1 and p is a polynomial of
d,k
folded CF estimator has error bounded by
degree k (see e.g. p.87 of [9] for explicit formulae).
(cid:113) Then the kernel K∗(x,y) = ϕd,k((cid:107)x−y(cid:107)) has native
E|Eb[f;u1:M,vM+1:N +∆]−I[f]|2 space Hd/2+k+1/2 (where the restriction d > 3 is in
≤CesHhα,bL(vM+1:N)h(u1:M)α−αLρ(u1:M)αU−αL(cid:107)f(cid:107)Hα ppr.1in0c9ipolfe[9re])q.uiWreidthfotrhtishekesrpneecliawlecacasenkth=er0ef)o(rseeeguea.gr-.
where antee a minimum level of smoothness. By rescaling,
thekernel’ssupportcanbechangedfromtheunitball
(cid:90)
(esh,b (vM+1:N))2 := e (b(vM+1:N +∆))2d∆ (as above) to balls of smaller radius. This in turn en-
HαL [0,1]d HαL forcessparsityonthesystemofinterpolationequations
thatarethebasisoftheCFestimatorandreducesthe
and C > 0 is a constant independent of f, vM+1:N computational cost of inverting this linear system.
and u1:M.
Wendland’skernelcannotbeintegratedanalyticallyin
Again, for quasi-uniform u1:M, CFs accelerate shifted d≥2 dimensions, violating requirement (ii). However
andfoldedαL-QMCbyafactorO(N−(α−αL)/d)(com- we can exploit recent work by [29] that shows the d-
pare against Sec. 5.9 of [7]). dimensional tensor product space Hk([0,1]) ⊗ ··· ⊗
Hk([0,1]) is norm-equivalent to SHk = SHk([0,1]d),
Scrambling: An explicit α-QMC rule that applies for
theSobolevspacewithdominatingmixedsmoothness:
all integer values of α was recently discovered by [6].
Forsimplicityfocussingond=1,theserandomdesign SHk := {f :[0,1]d →R|Daf exists and
points achieve α-rates and, moreover, the RMSE is Daf ∈L ([0,1]d),∀a∈Nd with a ≤k}.
p 0 i
controlled by a norm of the form (cid:107)f(cid:107) . When α
Hα (The distinction with Hk([0,1]d) is that the multi-
is known and is an integer, one may achieve optimal
index a is now constrained component-wise, a ≤ k,
ratesandCFsprovidenorateimprovement. However, i
when α∈/ N, CFs can be used to transform these sub- rather than |a| ≤ k.) In particular SHk([0,1]d) ⊆
Hk([0,1]d) so that functions in SHk are at least as
optimal integrators into optimal integrators.
smooth as functions in Hk. We therefore propose to
Choice of Kernel: The QMC+CF methodology has employ the product kernel
some flexibility in terms of the choice of kernel K
∗ d
that is used to construct the approximation fM. Our K(k)(x,y)=(cid:89)ϕ (|x −y |) (10)
∗ 1,k i i
main requirements here are: (i) K imposes ‘enough
∗ i=1
Control Functionals for Quasi-Monte Carlo Integration
whose native space is SHk+1. The integral results. In the 6th, discontinuous case the QMC+CF
(cid:90) method does not out-perform QMC (at least in di-
K∗(k)(x,y)dx mension d > 1), as the functional approximation fM
[0,1]d is poor due to violation of our continuity assump-
of tensor products of Wendland functions in Eqn. 10 tion. In all cases the performance of QMC+CF ap-
cannowbeintegratedanalytically. Thisapproachpro- proachesthatofQMCasthedimensiondincreased. In
vides a convenient mechanism to control the degree of higher dimensions (d ≥ 5, not shown) the QMC+CF
smoothnessthatweimposeontheapproximationf . and QMC estimators demonstrated effectively identi-
M
cal performance, in line with theory.
4 Experimental Results The experiments were then repeated with rougher
(k = 0) and smoother (k = 2) regression kernels. Re-
Our methodology provides a variance reduction tech- sults in the Supplement (Figs. S3-8) demonstrated a
nique for QMC that is able to accelerate convergence slight improvement in the performance of QMC+CF
rates, yet is also practical. The first numerical study when k = 2, in line with theory, though generally es-
belowisa‘proof-of-principle’designedtovalidatethis timates were robust to the choice of regression kernel.
specific claim in the empirical setting. To further assess the generality of these conclusions,
further experiments were performed using a different
Simulation Study: For objective assessment we ex-
QMC rule (truncated Sobol sequence with scrambling
ploited the test package proposed by [11]. This pack-
due to [18]). Results in the Supplement showed that
age defines 6 function families, each of them charac-
thesameconclusionscanbedrawnineachcase. Taken
terized by some peculiarity, such as oscillation, dis-
together, these results demonstrate that CFs can ac-
continuity or corner peaks, with the property that
celerate QMC, at least in low-dimensional settings,
their exact integrals are available. The ‘discontinu-
and thus complete our ‘proof-of-principle’. MATLAB
ous’Genzfunctionprovidesanexamplewheresmooth-
code to reproduce these results is provided.
ness assumptions on the test function are violated.
We used the MATLAB implementation of [11] that Application to Robot Arm Data: Todemonstrate
is freely available at http://people.sc.fsu.edu/ the benefits of our methodology we consider the prob-
~jburkardt/m_src/testpack/testpack.html. lem of estimating the inverse dynamics of a seven
degrees-of-freedom robot arm. The task, as described
In the experiments below, we focus on the two QMC
in [25], is to map from a 21-dimensional input space
rules that are most widely used in practice. In the
(7 positions, 7 velocities, 7 accelerations) to the cor-
first experiment, the random QMC point set vM+1:N
responding 7 joint torques. Following [25] we present
was generated by truncating the Halton sequence,
resultsbelowonjustoneofthemappings,fromthe21
scrambling the digits of the resulting points using the
input variables to the first of the seven torques. The
reverse-radix algorithm [16] and applying a uniform
datasetconsistsof48,933input-outputpairs,ofwhich
randomshift. ThisQMCruleachievestheα =1rate
L
44,484 were used as a training set and the remaining
on the subsequence N = 2n when the test function
n
4,449wereusedasatestset. Theinputswerelinearly
has mixed partial derivatives of first order. To ensure
rescaled to have mean zero and unit variance on the
thattheseQMCruleswereimplementedfaithfully, we
training set. The outputs were centred to have mean
restricted attention to the case where M = N/2 so
zero on the training set.
that N −M was always a power of two. The train-
ingpointsu1:M weretakentobed-dimensionalsquare We consider a hierarchical model based on 21-
lattices in all experiments. dimensional Gaussian process (GP) regression. De-
note by Y ∈ R a measured response variable at
We considered the 6 Genz functions in d = 1,2,3,4 i
state z ∈ R21, assumed to satisfy Y = g(z ) +
dimensions. The performance of QMC with and with- i i i
(cid:15) where (cid:15) ∼ N(0,σ2) are independent for i =
out CFs was compared, in each case ensuring that the i i
1,...,n and σ > 0 will be assumed known. In or-
total number of evaluations of the integrand f was
der to use training data (y ,z )n to make predic-
equal for all methods. For CFs, the tensor-product i i i=1
tions regarding an unseen test point z , we place a
Wendland kernel with k = 1 was employed (i.e. ap- ∗
GP prior g ∼ GP(0,c(z,z(cid:48);θ)) where c(z,z(cid:48);θ) =
proximation with functions f ∈ H2, so α = 2).
M U θ exp(−1θ−2(cid:107)z−z(cid:48)(cid:107)2). Here θ =(θ ,θ ) are hyper-
Results are presented in Fig. 1. (For clarity we chose 1 2 2 2 1 2
parameters that control how training samples are
notpresentresultsforMC,sincethesewereinferiorto
used to predict the response at a new test point.
QMC methods in all cases considered.) For the first 5
A fully-Bayesian treatment aims to marginalise over
Genz functions it holds that f ∈ Hα with α = 2 and
thesehyper-parametersandweassignindependentpri-
theory (for the random case) guarantees an accelera-
ors θ ∼ Γ(α,β), θ ∼ Γ(γ,δ) in the shape/scale
tion of O(N−1/d); this is borne out in experimental 1 2
Chris. J. Oates, Mark Girolami
Oscillatory Test Function (Halton Set, k = 1) Product Peak Test Function (Halton Set, k = 1) Corner Peak Test Function (Halton Set, k = 1)
100 100 100
10−2
Error10−5 Error10−5 Error10−4
Square Square Square 10−6
Mean 10−10 Mean 10−10 Mean 10−8
10−10
10−15 10−15 10−12
2 3 4 ... 10 20 30 40... 100 2 3 4 ... 10 20 30 40... 100 2 3 4 ... 10 20 30 40... 100
Num. Function Evaluations (N) Num. Function Evaluations (N) Num. Function Evaluations (N)
Gaussian Test Function (Halton Set, k = 1) Continuous Test Function (Halton Set, k = 1) Discontinuous Test Function (Halton Set, k = 1)
100 100 105
100
Error10−5 Error10−5 Error
Square Square Square 10−5
Mean 10−10 Mean 10−10 Mean
10−10
10−15 10−15 10−15
2 3 4 ... 10 20 30 40... 100 2 3 4 ... 10 20 30 40... 100 2 3 4 ... 10 20 30 40... 100
Num. Function Evaluations (N) Num. Function Evaluations (N) Num. Function Evaluations (N)
Figure 1: Simulation study (Genz functions): Each panel represents one test function. Solid lines correspond
to standard QMC, dashed lines correspond to QMC+CF. represents dimension d = 1, (cid:52) represents d = 2,
(cid:3) represents d = 3 and ∗ represents d = 4. Experiments w(cid:35)ere replicated with 10 random seeds and error bars
denote standard error of the replicate mean. QMC points were generated from a shifted and scrambled Halton
sequence. A Wendland regression kernel was used with k =1.
parametrisation, which we write jointly as π(θ). Here we took a random subset of n=1,000 training points
σ =0.1, α=β =γ =δ =2. and a subset of regressors approximation with n(cid:48) =
100. The total computational time needed to obtain
To predict the value of the response Y corresponding
∗ these results was 268 core-hours.
toanunseenstatevectorz , ourestimatorwillbethe
∗
Bayesian posterior mean Foreachtestpointz thesamplingstandarddeviation
∗
of Yˆ was estimated from 10 independent realisations
(cid:90) ∗
Yˆ :=E[Y |y]= E[Y |y,θ]π(θ)dθ, (11) of the QMC procedures. For CF we used a randomly-
∗ ∗ ∗
shifted, scrambled Halton sequence (α = 1) and
L
Wendland kernels with k = 1 (α = 2), so that the-
where we implicitly condition on the covariates U
ory predicts an acceleration factor of O(N−1/2). The
z ,...,z ,z . Phrasing in terms of our earlier nota-
1 n ∗
estimator standard deviations were estimated for all
tion, the test function is
4,449 test points (with N = 28) and the full results
f(x)=E[Y∗|y,Π−1(x)]=C∗,n(Cn+σ2In×n)−1y are shown in Fig. 2. Note that each test point z∗ cor-
responds to a different test function f and thus these
where Π is the c.d.f for π, (Cn)i,j = c(zi,zj;θ) and resultsarequiteobjective,encompassingthousandsof
(C∗,n)1,j = c(z∗,zj;θ). Each evaluation of the in- different integration problems. For the vast major-
tegrand f(x) requires O(n3) operations due to the ity of integration problems, CF accelerated the stan-
matrix inversion and this entails a prohibitive level dard QMC estimator. Here the computational time
of computation. A partial solution is provided by a to construct a functional approximation (inverting a
‘subset of regressors’ approximation 16 × 16 matrix) was negligible (3%) in comparison
to the cost of evaluating the function f once. The
f(x)≈C (C C +σ2C )−1C y (12)
∗,n(cid:48) n(cid:48),n n,n(cid:48) n(cid:48) n(cid:48),n total additional computational time associated with
the QMC+CF methodology was 2% greater than for
where n(cid:48) < n denotes a subset of the full data; see
QMC, which is easily justified by the substantial vari-
Sec. 8.3.1 of [25] for full details. However even Eqn.
ance reductions (∼103%) that are realised in this ap-
12 still represents a substantial computational burden
plication. Supplementary results (Fig. S9) compare
in general. To facilitate the illustration below, which
QMC+CF to MC+CF (standard MC sampling).
investigates the sampling distribution of estimators,
Control Functionals for Quasi-Monte Carlo Integration
Test Function (1D Projection)
1 0.05
x) 0
f(
−1
0 0.5 1 0.04
x
1
Tractable Approximation
1 C] 0.03
f(x)N −010 0.5 1 Std[QM 0.02
x
1
(x)N 1Difference, for QMC 0.01
− f 0
x) −1 0
f( 0 0.5 1 0 0.01 0.02 0.03 0.04 0.05
x Std[QMC + CF]
1
Figure 2: Application to robot arm data. Left: Posterior predictive means were computed for the mechanical
torque experienced by one of the seven joints of the arm, for each of 4,449 joint configurations. Schematic
reproduced from [32]. Centre: Model hyper-parameters were integrated out; for this task we compared standard
QMC with the proposed QMC+CF approach (both implementations provided unbiased estimators). Right:
Examining the estimator sampling standard deviations, we see that, for all but a handful of the configurations,
QMC+CF was more accurate than QMC.
5 Discussion ple obtain dimension-independent rates by imposing a
(strong) assumption of polynomial tractability on the
QMC methods are becoming increasingly relevant in RKHS. This is achieved by generalising to weighted
modern statistics applications [12, 35] and it is surely Sobolev spaces, such that the integrand f ‘depends
a priority to target the rate constants governing the onlyweaklyonmostofthecomponentsofx’. Further
practical performance of these algorithms. CFs pro- details are provided in [7, 20] and form part of our
videoneroutetoachievethisgoal, providingsubstan- ongoing research.
tial variance reductions in many of the examples we
The methods that we describe are immediately appli-
considered. Indeed, CFsallowustouseasub-optimal
cable in a range of applications including marginalisa-
QMC rule (e.g. as built into existing software pack-
tion of hyper-parameters in classification [10], proba-
ages) and yet, with minimal additional coding, obtain
bilisticinferencefordifferentialequations[28,5],com-
a QMC+CF algorithm that attains optimal conver-
putation of model evidence [21] and approximation of
gence rates. The focus on unknown smoothness α dis-
the partition function in social network models [26].
tinguishes our work from previous literature on the
Finally we note that CFs generalise to other integra-
connectionbetweenintegrationandfunctionalapprox-
tion methods including Bayesian Quadrature [23, 4]
imation, e.g. [3, 13].
andrelatedkernel-basedquadraturerules[2],inwhich
Functional approximation, and hence our QMC+CF the worst case error is also controlled by an RKHS
methodology, has a computational cost associated norm (cid:107)f(cid:107) ; this will be the focus of our ongoing re-
H
with solution of a linear system. Whilst negligible in search.
our experiments, this cost could be reduced if neces-
saryusingstandardapproximationsand/orcompactly
supported kernels. On the other hand, we note that
QMC is often used when f is expensive to evaluate Acknowledgments
and in such situations it is likely that evaluation of
the integrand, rather than solution of a linear system,
The authors are grateful to Dan Simpson, Mathieu
will be the main computational bottleneck.
Gerber and Ben Collyer for helpful discussions. CJO
OurfocuswasonSobolevspaces, butitisknownthat was supported by EPSRC [EP/D002060/1] and the
a faster rate O(N−α+(cid:15)) is possible in the subspace ARC Centre of Excellence for Mathematics and Sta-
SHα([0,1]d), for any (cid:15) > 0, and explicit point sets tistical Frontiers. MG was supported by EPSRC
are available (for integer α) [6]. An immediate exten- [EP/J016934/1, EP/K034154/1], an EPSRC Estab-
sion is to establish optimal rates for CFs in this class lished Career Fellowship, the EU grant [EU/259348]
of functions. In a related direction, one can in princi- and a Royal Society Wolfson Research Merit Award.
Chris. J. Oates, Mark Girolami
References [14] F.J. Hickernell, C. Lemieux, and A.B. Owen.
Control variates for quasi-Monte Carlo. Statis-
[1] C. Aistleitner and J. Dick. Functions of
tical Science, 20(1):1–31, 2005.
bounded variation, signed measures, and a gen-
eral Koksma-Hlawka inequality. Acta Arith- [15] F.Husz´arandD.Duvenaud.Optimally-Weighted
metica, 167(2):143–171, 2015. Herding is Bayesian Quadrature. In Uncertainty
in Artificial Intelligence, pages 377–385, 2012.
[2] F. Bach. On the equivalence between quadrature
rules and random features. arXiv:1502.06800, [16] L. Kocis and W.J. Whiten. Computational in-
2015. vestigations of low-discrepancy sequences. ACM
TransactionsonMathematicalSoftware(TOMS),
[3] N.S. Bakhvalov. On the approximate calcula-
23(2):266–294, 1997.
tion of multiple integrals. Journal of Complexity,
31(4):502-516, 2015. [17] S. Lacoste-Julien, F. Lindsten, and F. Bach. Se-
quential Kernel Herding: Frank-Wolfe Optimiza-
[4] F.X. Briol, C.J. Oates, M. Girolami, M. Os-
tion for Particle Filtering. In International Con-
borne and D. Sejdinovic. Probabilistic Integra-
ference on Artificial Intelligence and Statistics,
tion: A Role for Statisticians in Numerical Anal-
volume 18, 2015.
ysis? arXiv:1512.00933, 2016.
[18] J.Matouˇsek.OntheL2-discrepancyforAnchored
[5] J. Cockayne, C.J. Oates, T. Sullivan and
Boxes. J. Complex., 14(4):527–556, 1998.
M. Girolami. Probabilistic Meshless Methods for
Bayesian Inverse Problems. In preparation. [19] E. Novak. Deterministic and stochastic error
bounds in numerical analysis. Springer-Verlag
[6] J. Dick. Higher order scrambled digital nets
Berlin, 1988.
achieve the optimal rate of the root mean square
errorforsmoothintegrands.TheAnnalsofStatis- [20] E. Novak and H. Wo´zniakowski. Tractability of
tics, 39(3):1372–1398, 2011. Multivariate Problems: Standard information for
functionals, volume 2. European Mathematical
[7] J. Dick, F.Y. Kuo, and I.H. Sloan. High-
Society, 2010.
dimensional integration: The quasi-Monte Carlo
way. Acta Numerica, 22:133–288, 2013. [21] C.J.Oates,T.PapamarkouandM.Girolami.The
Controlled Thermodynamic Integral for Bayesian
[8] J. Dick and F. Pillichshammer. Discrepancy the-
Model Evidence Evaluation. J. Am. Stat. Assoc.,
ory and quasi-Monte Carlo integration. Springer,
2016. To appear.
Berlin, 2010.
[22] C.J.Oates,M.Girolami,andN.Chopin. Control
[9] G.F.Fasshauer. Meshfree approximation methods
Functionals for Monte Carlo Integration. J. R.
with MATLAB. World Scientific Publishing Co.,
Statist. Soc. B, 2016. To appear.
Inc., 2007.
[23] A. O’Hagan. Bayes-Hermite quadrature. Journal
[10] M. Filippone and M. Girolami. Pseudo-Marginal
of Statistical Planning and Inference, 29(3):245–
Bayesian Inference for Gaussian Processes. Pat-
260, 1991.
tern Analysis and Machine Intelligence, IEEE
Transactions on, 36(11):2214–2226, 2014. [24] A.B. Owen. A constraint on extensible quadra-
ture rules. Numerische Mathematik, 2015. To
[11] A. Genz. Testing multidimensional integration
appear.
routines. In International Conference on Tools,
Methods and Languages for Scientific and En- [25] C.E. Rasmussen and C.K.I. Williams. Gaussian
gineering Computation, pages 81–94. Elsevier processesformachinelearning.MITPress,2(3):4,
North-Holland, Inc., 1984. 2006.
[12] M. Gerber and N. Chopin. Sequential Quasi- [26] G.Robins,P.Pattison,Y.Kalish,andD.Lusher.
Monte Carlo. Journal of the Royal Statisti- An introduction to exponential random graph
cal Society: Series B (Statistical Methodology), (p*) models for social networks. Social Networks,
77(3):509–579, 2015. 29(2):173–191, 2007.
[13] S. Heinrich. Random approximation in numeri- [27] R. Schaback and H. Wendland. Kernel tech-
cal analysis. Lecture Notes in Pure and Applied niques: from machine learning to meshless meth-
Mathematics, 150:123-171, 1994 ods. Acta Numerica, 15:543–639, 2006.
Control Functionals for Quasi-Monte Carlo Integration
[28] M. Schober, D.K. Duvenaud, and P. Hennig.
Probabilistic ODE solvers with Runge-Kutta
means. In Advances in Neural Information Pro-
cessing Systems, volume27, pages739–747.2014.
[29] W. Sickel and T. Ullrich. Tensor products of
Sobolev-Besovspacesandapplicationstoapprox-
imationfromthehyperboliccross. JournalofAp-
proximation Theory, 161(2):748–786, 2009.
[30] A.Smola,A.Gretton,L.Song,andB.Sch¨olkopf.
A Hilbert space embedding for distributions.
In Algorithmic Learning Theory, pages 13–31.
Springer, 2007.
[31] B.Tracey,D.Wolpert,andJ.J.Alonso. UsingSu-
pervised Learning to Improve Monte Carlo Inte-
gralEstimation. AIAAJournal,51(8):2015–2023,
2013.
[32] S. Vijayakumar and S. Schaal. Locally weighted
projection regression: An O(n) algorithm for in-
cremental real time learning in high dimensional
space. In International Conference on Machine
Learning, volume 16, pages 1079–1086, 2000.
[33] X.Wang.EnhancingQuasi-MonteCarloMethods
by Exploiting Additive Approximation for Prob-
lemsinFinance.SIAMJournalofScientificCom-
puting, 34(1):A283–A308, 2012.
[34] H.Wendland. Piecewisepolynomial,positivedef-
inite and compactly supported radial functions
of minimal degree. Advances in Computational
Mathematics, 4(1):389–396, 1995.
[35] J. Yang, V. Sindhwani, H. Avron, and M. Ma-
honey. Quasi-MonteCarlofeaturemapsforshift-
invariant kernels. In International Conference
on Machine Learning, volume 31, pages 485–493,
2014.