DiCE: The Infinitely Differentiable Monte Carlo Estimator

Jakob Foerster 1 Gregory Farquhar * 1 Maruan Al-Shedivat * 2 Tim Rocktäschel 1 Eric P. Xing 2
Shimon Whiteson 1

Abstract

arXiv:1802.05098v1 [cs.LG] 14 Feb 2018

The score function estimator is widely used for
estimating gradients of stochastic objectives in
Stochastic Computation Graphs (SCG), e.g., in
reinforcement learning and meta-learning. While
deriving the first order gradient estimators by differentiating a surrogate loss (SL) objective is computationally and conceptually simple, using the
same approach for higher order gradients is more
challenging. Firstly, analytically deriving and implementing such estimators is laborious and not
compliant with automatic differentiation. Secondly, repeatedly applying SL to construct new
objectives for each order gradient involves increasingly cumbersome graph manipulations. Lastly,
to match the first order gradient under differentiation, SL treats part of the cost as a fixed sample,
which we show leads to missing and wrong terms
for higher order gradient estimators. To address
all these shortcomings in a unified way, we introduce D I CE, which provides a single objective
that can be differentiated repeatedly, generating
correct gradient estimators of any order in SCGs.
Unlike SL, D I CE relies on automatic differentiation for performing the requisite graph manipulations. We verify the correctness of D I CE both
through a proof and through numerical evaluation of the D I CE gradient estimates. We also use
D I CE to propose and evaluate a novel approach
for multi-agent learning. Our code is available at
https://goo.gl/xkkGxN.

1. Introduction
The score function trick is used to produce Monte Carlo
estimates of gradients in settings with non-differentiable objectives, e.g., in meta-learning and reinforcement learning.
Estimating the first order gradients is computationally and
*

Order determined by a dice roll. 1 University of Oxford
2
Carnegie Mellon University. Correspondence to:
Jakob Foerster <jakob.foerster@cs.ox.ac.uk>.

conceptually simple. While the gradient estimators can be
directly defined, it is often more convenient to define an
objective whose derivative is the gradient estimator and let
the powerful automatic-differentiation (auto-diff) toolbox
as implemented in deep learning libraries do the work for
you. This is the method used by the surrogate loss (SL)
approach (Schulman et al., 2015a), which provides a recipe
for building a surrogate objective from a stochastic computation graph (SCG). When differentiated, the SL produces an
estimator for the first order gradient of the original objective.
However, estimating higher order gradients is more challenging. Such estimators are useful for a number of optimization
techniques, accelerating convergence in supervised settings
(Dennis & Moré, 1977). Furthermore, they are vital for
gradient-based meta-learning (Finn et al., 2017; Al-Shedivat
et al., 2017; Li et al., 2017), which differentiates an objective
after some number of first order learning steps. Higher order
gradient estimators have also proven useful in multi-agent
learning (Foerster et al., 2018), when one agent differentiates through the learning process of another agent.
Unfortunately, the first order gradient estimators mentioned
above are fundamentally ill-suited for calculating higher
order derivatives via auto-diff. Due to the dependency on
the sampling distribution, higher order gradient estimators
require repeated application of the score function trick. Simply differentiating the first order estimator again, as was for
example done by Finn et al. (2017), leads to missing terms.
To obtain higher order score function gradient estimators,
there are currently two unsatisfactory options. The first is
to analytically derive and implement the estimators. However, this is laborious, error prone, and does not comply
with the auto-diff paradigm. The second is to repeatedly
apply the SL approach to construct new objectives for each
further gradient estimate. However, constructing each of
these new objectives involves increasingly complex graph
manipulations, defeating the appeal of using a differentiable
surrogate loss.
Moreover, to match the first order gradient after a single
differentiation, the SL treats part of the cost as a fixed sample, severing the dependency on the parameters. We show
that this yields missing and incorrect terms in higher order
gradient estimators. We believe that these difficulties have

DiCE: The Infinitely Differentiable Monte Carlo Estimator

limited the usage and exploration of higher order methods
in reinforcement learning tasks and other application areas
that may be formulated as SCGs.
Therefore, we propose a novel technique, the Infinitely
Differentiable Monte-Carlo Estimator (D I CE), to address
all these shortcomings. D I CE constructs a single objective
that evaluates to an estimate of the original objective, but
can also be differentiated repeatedly to obtain correct gradient estimators of any order. Unlike the SL approach, D I CE
relies on auto-diff as implemented for instance in TensorFlow (Abadi et al., 2016) or PyTorch (Paszke et al., 2017)
to automatically perform the complex graph manipulations
required for these higher order gradient estimators.
D I CE uses a novel operator, M AGIC B OX ( ), which wraps
around the set of those stochastic nodes Wc that influence
each of the original losses in an SCG. Upon differentiation,
this operator generates the correct gradients associated with
the sampling distribution:
X
∇θ (Wc ) = (Wc )∇θ
log(p(w; θ)),
w ∈ Wc

while returning 1 when evaluated: (W)  1. The
M AGIC B OX-operator can easily be implemented in standard deep learning libraries as follows:

(W) = exp τ − ⊥(τ ) ,
X
τ=
log(p(w; θ)),
w∈W

where ⊥ is an operator that sets the gradient of the operand
to zero, so ∇x ⊥(x) = 0. In addition, we show how to use a
baseline for variance reduction in our formulation.
We verify the correctness of D I CE both through a proof
and through numerical evaluation of the D I CE gradient estimates. To demonstrate the utility of D I CE, we also propose
a novel approach for learning with opponent learning awareness (Foerster et al., 2018). We also open-source our code
in TensorFlow. We hope this powerful and convenient novel
objective will unlock further exploration and adoption of
higher order learning methods in meta-learning, reinforcement learning, and other applications of SCGs.

2. Background
Suppose x is a random variable, x ∼ p(x; θ), f is a function
of x and we want to compute ∇θ Ex [f (x)]. If the analytical gradients ∇θ f are unavailable or nonexistent, we can
employ the score function (SF) estimator:
∇θ Ex [f (x)] = Ex [f (x)∇θ log(p(x; θ))]

(2.1)

If instead x is a deterministic function of θ and another random variable z, the operators ∇θ and Ez commute, yielding

the pathwise derivative estimator or reparameterisation
trick. In this work, we focus on the SF estimator, which can
capture the interdependency of both the objective and the
sampling distribution on the parameters θ, and therefore requires careful handling for higher order gradient estimates.
2.1. Stochastic Computation Graphs
Gradient estimators for single random variables can be generalised using the formalism of a stochastic computation
graph (SCG, Schulman et al., 2015a). An SCG is a directed
acyclic graph with four types of nodes: input nodes, Θ;
deterministic nodes, D; cost nodes, C; and stochastic nodes,
S. Input nodes are set externally and can hold parameters
we seek to optimise. Deterministic nodes are functions of
their parent nodes, while stochastic nodes are distributions
conditioned on their parent nodes. The set of P
cost nodes C
are those associated with an objective L = E[ c∈C c].
Let v ≺ w denote that node v influences node w, i.e., there
exists a path in the graph from v to w. If every node along
the path is deterministic, v influences w deterministically
which is denoted by v ≺D w. See Figure 1 (top) for a
simple SCG with an input node θ, a stochastic node x and a
cost function f . Note that θ influences f deterministically
(θ ≺D f ) as well as stochastically via x (θ ≺ f ).
2.2. Surrogate Losses
In
P order to estimate gradients of a sum of cost nodes,
c∈C c, in an arbitrary SCG, Schulman et al. (2015a) introduce the notion of a surrogate loss (SL):
X
X
SL(Θ, S) :=
log p(w | DEPSw )Q̂w +
c(DEPSc )
w∈S

c∈C

Here DEPSw are the “dependencies” of w: the set of stochastic or input nodes that deterministically influence the node
w. Furthermore, Q̂w is the sum of sampled costs ĉ corresponding to the cost nodes influenced by w.
The SL produces a gradient estimator when differentiated
once (Schulman et al., 2015a, Corollary 1):
∇θ L = E[∇θ SL(Θ, S)].

(2.2)

The hat notation on Q̂w indicates that, inside the SL, these
costs are treated as fixed samples, thus severing the functional dependency on θ that was present in the original
stochastic computation graph. This ensures that the first
order gradients of the SL match the score function estimator,
which does not contain a term of the form log(p)∇θ Q.
Although Schulman et al. (2015a) focus on first order gradients, they argue that the SL gradient estimates themselves
can be treated as costs in an SCG and that the SL approach
can be applied repeatedly to construct higher order gradient

DiCE: The Infinitely Differentiable Monte Carlo Estimator
Stochastic Computation Graph

estimators. However, the use of sampled costs in the SL
leads to missing dependencies and wrong estimates when
calculating such higher order gradients, as we discuss in
Section 3.2.

x

θ

Surrogate Loss Approach

θ

x

f

D I CE

f

θ

x

f

(x)f

3. Higher Order Gradients
∇

In this section, we illustrate how to estimate higher order
gradients via repeated application of the score function (SF)
trick and show that repeated application of the surrogate
loss (SL) approach in stochastic computation graphs (SCGs)
fails to capture all of the relevant terms for higher order
gradient estimates.

∇θ L

log(p(x; θ))fˆ + f

∇2θ L

log(p(x; θ))ĝSL +
fˆ∇θ log(p(x; θ)) + ∇θ f

∇

gSL = fˆ∇θ log(p(x; θ)) + ∇θ f

∇θ ( (x)f )
∇

∇

ĝSL ∇θ log(p(x; θ)) +
fˆ∇2θ log(p(x; θ)) +
ˆ
∇f ∇θ log(p(x; θ)) + ∇2θ f

=0
∇nθ L

∇2θ ( (x)f )
∇
···

∇nθ ( (x)f )

3.1. Higher Order Gradient Estimators
We begin by revisiting the derivation of the score function
estimator for the gradient of the expectation L of f (x; θ)
over x ∼ p(x; θ):
∇θ L = ∇θ Ex [f (x; θ)]
X
= ∇θ
p(x; θ)f (x; θ)
x

=

X

=

X

Figure 1. Simple example illustrating the difference of the Surrogate Loss (SL) approach to D I CE. Stochastic nodes are depicted in
orange, costs in gray, surrogate losses in blue, D I CE in purple, and
gradient estimators in red. Note that for second-order gradients, SL
requires the construction of an intermediate stochastic computation
graph and due to taking a sample of the cost ĝSL , the dependency
on θ is lost, leading to an incorrect second-order gradient estimator.
Arrows from θ, x and f to gradient estimators omitted for clarity.


∇θ p(x; θ)f (x; θ)

x


f (x; θ)∇θ p(x; θ) + p(x; θ)∇θ f (x; θ)

bution. This leads to missing terms in the estimators of the
higher order gradients.

f (x; θ)p(x; θ)∇θ log(p(x; θ))

Consider the following example, where a single parameter
θ defines a sampling distribution p(x; θ) and the objective
is f (x, θ).

x

=

X
x


+ p(x; θ)∇θ f (x; θ)
= Ex [f (x; θ)∇θ log(p(x; θ)) + ∇θ f (x; θ)] (3.1)
= Ex [g(x; θ)].
The estimator g(x; θ) of the gradient of Ex [f (x; θ)] consists
of two distinct terms: (1) The term f (x; θ)∇θ log(p(x; θ))
originating from f (x; θ)∇θ p(x; θ) via the SF trick, and (2)
the term ∇θ f (x; θ), due to the direct dependence of f on θ.
The second term is often ignored because f is often only a
function of x but not of θ. However, even in that case, the
gradient estimator g depends on both x and θ. We might be
tempted to again apply the SL approach to ∇θ Ex [g(x; θ)] to
produce estimates of higher order gradients of L, but below
we demonstrate that this fails. In Section 4, we subsequently
introduce a practical algorithm for correctly producing such
higher order gradient estimators in SCGs.
3.2. Higher Order Surrogate Losses
While Schulman et al. (2015a) focus on the first order gradients, they state that a recursive application of SL can
generate higher order gradient estimators. However, as we
demonstrate in this section, because the SL approach treats
part of the objective as a sampled cost, the corresponding
terms lose a functional dependency on the sampling distri-

SL(L) = log p(x; θ)fˆ(x) + f (x; θ)
(∇θ L)SL = Ex [∇θ SL(L)]
= Ex [fˆ(x)∇θ log p(x; θ) + ∇θ f (x; θ)] (3.2)
= Ex [gSL (x; θ)].
The corresponding SCG is depicted at the top of Figure 1.
Comparing (3.1) and (3.2), note that the first term, fˆ(x) has
lost its functional dependency on θ, as indicated by the hat
notation and the lack of a θ argument. While these terms
evaluate to the same estimate of the first order gradient,
the lack of the functional dependency yields a discrepancy
between the exact derivation of the second order gradient
and a second application of SL.
SL(gSL (x; θ)) = log p(x; θ)ĝSL (x) + gSL (x; θ)
(∇2θ L)SL = Ex [∇θ SL(gSL )]
= Ex [ĝSL (x)∇θ log p(x; θ) + ∇θ gSL (x; θ)].
(3.3)
By contrast, the exact derivation of ∇2θ L results in the fol-

DiCE: The Infinitely Differentiable Monte Carlo Estimator

lowing expression:
∇2θ L = ∇θ Ex [g(x; θ)]
= Ex [g(x; θ)∇θ log p(x; θ) + ∇θ g(x; θ)].

(3.4)

Since gSL (x; θ) differs from g(x; θ) only in its functional
dependencies on θ, gSL and g are identical when evaluated.
However, due to the missing dependencies in gSL , the gradients w.r.t. θ, which appear in the higher order gradient
estimates in (3.3) and (3.4), differ:
∇θ g(x; θ) = ∇θ f (x; θ)∇θ log(p(x; θ))
+ f (x; θ)∇2θ log(p(x; θ))
+ ∇2θ f (x; θ),
∇θ gSL (x; θ) = fˆ(x)∇2θ log(p(x; θ))
+ ∇2θ f (x; θ).
We lose the term ∇θ f (x; θ)∇θ log(p(x; θ)) in the second
order SL gradient because ∇θ fˆ(x) = 0 (see left part of
Figure 1). This issue occurs immediately in the second
order gradients when f depends directly on θ. However, as
g(x; θ) always depends on θ, the SL approach always fails
to produce correct third or higher order gradient estimates
even if f depends only indirectly on θ.
3.3. Simple Failing Example
Here is a toy example to illustrate a possible failure case.
Let x ∼ Ber(θ) and f (x, θ) = x(1 − θ) + (1 − x)(1 + θ).
For this simple example we can exactly evaluate all terms:
L = θ(1 − θ) + (1 − θ)(1 + θ)
= −2θ2 + θ + 1
∇θ L = −4θ + 1
∇2θ L = −4

Evaluating the expectations for the SL-gradient estimators
analytically results in the following terms:
(∇θ L)SL = −4θ + 1
(∇2θ L)SL

= −2

Even with an infinite number of samples, the SL estimator
produces the wrong second order gradient. If, for example,
these wrong estimates were used in combination with the
Newton-Raphson method for optimising L, then θ would
never converge to the correct value. In contrast, this method
would converge in a single step using the correct gradients.
The failure mode seen in this toy example will appear whenever the objective includes a regularisation term that depends

on θ, and is also impacted by the stochastic samples. One
example in a practical algorithm is soft Q-learning for RL
(Schulman et al., 2017), which regularises the policy by
adding an entropy penalty to the rewards. This penalty encourages the agent to maintain an exploratory policy, reducing the probability of getting stuck in local optima. Clearly
the penalty depends on the policy parameters θ. However,
the policy entropy will also depends on the states visited,
which in turn depend on the stochastically sampled actions.
As a result, the entropy regularised RL objective in this algorithm will have the exact property leading to the failure
of the SL approach shown above. Unlike our toy analytic
example, the consequent errors will not just appear as a
rescaling of the proper higher order gradients, but will depend in a complex way on the parameters θ. Any second
order methods with such a regularised objective will therefore require an alternative strategy for generating gradient
estimators, even setting aside the awkwardness of repeatedly
generating new surrogate objectives.

4. Correct Gradient Estimators with DiCE
In this section, we propose the Infinitely Differentiable
Monte-Carlo Estimator (D I CE), a practical algorithm for
programatically generating correct gradients of any order
in arbitrary SCGs. The naive option is to recursively apply the update rules in (3.1) that map from f (x; θ) to the
estimator of its derivative g(x; θ). However, this approach
has two deficiencies: First, by defining gradients directly,
it fails to provide an objective that can be used in standard deep learning libraries. Second, these naive gradient
estimators violate the auto-diff paradigm for generating further estimators by repeated differentiation since in general
∇θ f (x; θ) 6= g(x; θ). Our approach addresses these issues,
as well as fixing the missing terms from the SL approach.
P
As before, L = E[ c∈C c] is the objective in an SCG. The
correct expression for a gradient estimator that preserves all
required dependencies for further differentiation is:
"
∇θ L = E

X
c∈C

c

X

∇θ log p(w | DEPSw )

w∈Wc

!#
+ ∇θ c(DEPSc )

,

(4.1)

where Wc = {w | w ∈ S, w ≺ c, θ ≺ w}, i.e. the set of
stochastic nodes that depend on θ and influence the cost c.
For brevity, from here on we suppress the DEPS notation,
assuming all probabilities and costs are conditioned on their
relevant ancestors in the SCG.
Note that (4.1) is the generalisation of (3.1) to arbitrary
SCGs. The proof is given by Schulman et al. (2015a, Lines
1-10, Appendix A). Crucially, in Line 11 the authors then

DiCE: The Infinitely Differentiable Monte Carlo Estimator

replace c by ĉ, severing the dependencies required for correct higher order gradient estimators. As described in Section 2.2, this was done so that the SL approach reproduces
the score function estimator after a single differentiation
and can thus be used as an objective for backpropagation in
a deep learning library.
To support correct higher order gradient estimators, we
propose D I CE, which relies heavily on a novel operator,
M AGIC B OX ( ). M AGIC B OX takes a set of stochastic
nodes W as input and has the following two properties by
design:

∇θ cn = cn+1 , which we prove as follows:
∇θ cn = ∇θ (cn (Wcn ))
= cn ∇θ (Wcn ) + (Wcn )∇θ cn


X
= cn (Wcn ) 
∇θ log(p(w; θ))
w∈Wcn

+
=

(Wcn )∇θ cn



X

(Wcn ) ∇θ cn + cn

∇θ log(p(w; θ))

w∈Wcn

(4.4)
1.

(W)  1,

2. ∇θ (W) =

=
(W)

P

w∈W

∇θ log(p(w; θ)).

Here,  indicates “evaluates to” in contrast to full equality,
=, which includes equality of all gradients. In the autodiff paradigm,  corresponds to a forward pass evaluation
of a term. Meanwhile, the behaviour under differentiation
in property (2) indicates the new graph nodes that will be
constructed to hold the gradients of that object. Note that
that (W) reproduces the dependency of the gradient on
the sampling distribution under differentiation through the
requirements above. Using , we can next define the D I CE
objective, L :
L =

X

(Wc )c.

(4.2)

c∈C

Below we prove that the D I CE objective indeed produces
correct arbitrary order gradient estimators under differentiation.
Theorem 1. E[∇nθ L ]  ∇nθ L, ∀n ∈ {0, 1, 2, . . . }.

(Wcn+1 )c

n+1

=c

n+1

.

(4.5)

To proceed from (4.4) to (4.5), we need two additional
steps. First, we require an expression for cn+1 . Substituting L = E[cn ] into (4.1) and comparing to (4.3), we find
the following map from cn to cn+1 :
X
cn+1 = ∇θ cn + cn
∇θ log p(w; θ).
(4.6)
w∈Wcn

The term inside the brackets in (4.4) is identical to cn+1 .
Secondly, note that (4.6) shows that cn+1 depends only
on cn and Wcn . Therefore, the stochastic nodes which
influence cn+1 are the same as those which influence cn . So
Wcn = Wcn+1 , and we arrive at (4.5).
To conclude the proof, recall that cn is the estimator for the
nth derivative of c, and that cn  cn . Summing over c ∈ C
then gives the desired result.
Implementation of D I CE. D I CE is easy to implement in
standard deep learning libraries:

(W) = exp τ − ⊥(τ ) ,
X
τ=
log(p(w; θ)),
w∈W

Proof. For each cost node c ∈ C, we define a sequence of
nodes, cn , n ∈ {0, 1, . . . } as follows:

Since ⊥(x)  x, clearly

c0 = c,
E[cn+1 ] = ∇θ E[cn ].

(4.3)

By induction it follows that E[cn ] = ∇nθ E[c] ∀n, i.e. that
cn is an estimator of the nth order derivative of the objective
E[c].
n

where ⊥ is an operator that sets the gradient of the operand
to zero, so ∇x ⊥(x) = 0. 1

∇θ

(W)  1. Furthermore:

(W) = ∇θ exp τ − ⊥(τ )

= exp τ − ⊥(τ ) ∇θ (τ − ⊥(τ ))
=
=

(W)(∇θ τ + 0)
X
(W)
∇θ log(p(w; θ)).
w∈W

n

We further define c = c (Wcn ). Since (x)  1,
clearly cn  cn . Therefore E[cn ]  E[cn ] = ∇nθ E[c],
i.e., cn is also a valid estimator of the nth order derivative
of the objective. Next, we show that cn can be generated
by differentiating c0 n times. This follows by induction, if

With this implementation of the -operator, it is now
straightforward to construct L as defined in (4.7). This
1
This operator exists in PyTorch as detach and in TensorFlow
as stop gradient.

DiCE: The Infinitely Differentiable Monte Carlo Estimator

procedure is demonstrated in Figure 2, which shows a reinforcement learning use case. In this example, the cost
nodes are rewards that depend
P on stochastic actions, and
the total objective
is
J
=
E[
rt ]. We construct a D I CE
P
objective J = t ({at0 , t0 ≤ t})rt . Now E[J ]  J
and E[∇nθ J ]  ∇nθ J, so J can both be used to estimate
the return and to produces estimators for any order gradients
under auto-diff, which can be used for higher order methods
such as TRPO (Schulman et al., 2015b).
Causality. In Theorem 1 of Schulman et al. (2015a), two
expressions for the the gradient estimator are provided:
1. The first expression handles causality by summing over
stochastic nodes, w, and multiplying ∇ log(p(w)) for
each stochastic node with a sum of the downstream
costs, Q̂w . This is forward looking causality.
2. In contrast, the second expression sums over costs, c,
and multiplies each cost with a sum over the gradients
of
P log-probabilities from upstream stochastic nodes,
w∈Wc ∇ log(p(w)). We can think of this as backward looking causality.
In both cases, integrating causality into the gradient estimator leads to reduction of variance compared to the naive
approach of multiplying the full sum over costs with the full
sum over grad-log-probabilities.
While the SL approach is based on the first expression,
D I CE uses the second formulation. As shown by Schulman
et al. (2015a), both expressions result in the same terms
for the gradient estimator. However, the second formulation leads to greatly reduced conceptual complexity when
calculating higher order terms, which we exploit in the definition of the D I CE objective. This is because each further
gradient estimator maintains the same backward looking
dependencies for each term in the original sum over costs,
i.e., Wcn = Wcn+1 . In contrast, the SL approach is centred
around the stochastic nodes, which each become associated with a growing number of downstream costs after each
differentiation. Consequently, we believe that our D I CE
objective is more intuitive, as it is conceptually centred
around the original objective and remains so under repeated
differentiation.
Variance Reduction. So far, the construction of the D I CE
objective addresses variance reduction only via implementing causality, since each cost term is associated with the
that captures all causal dependencies. However, we can
also include a baseline term in the definition of the D I CE
objective:
X
X
L =
(Wc )c +
(1 − ({w}))bw .
(4.7)
c∈C

w∈S

The baseline bw is a design choice and can be any function

θ
s1

s2

···

st

a1

a2

···

at

r1

r2

···

rt

(a1 )r1

(a1 , a2 )r2

∇n
∇nθ (a1 )r1

···

∇n
∇nθ (a1 , a2 )r2

(a1 , . . . , at )rt
∇n

· · · ∇nθ (a1 , . . . , at )rt

Figure 2. D I CE applied to a reinforcement learning problem. A
stochastic policy conditioned on st and θ produces actions, at ,
which lead to rewards rt and next states, st+1 . Associated with
each reward is a D I CE objective that takes as input the set of
all causal dependencies that are functions of θ, i.e., the actions.
Arrows from θ, ai and ri to gradient estimators omitted for clarity.

of nodes not influenced by w. As long as this condition
is met, the baseline will not change the expectation of the
gradient estimates, but can considerably reduce the variance
(including of higher order gradient estimators). A common
choice is the average cost.
Since (1− ({w}))  0, the addition of the baseline leaves
the evaluation of the estimator L of the original objective
unchanged, while reproducing the baseline term stated by
Schulman et al. (2015a) after differentiation.
Hessian-Vector Product. The Hessian-vector, v > H, is
useful for a number of algorithms, such as estimation of
eigenvectors and eigenvalues of H (Pearlmutter, 1994). Using D I CE, v > H can be implemented efficiently without
having to compute the full Hessian. Assuming v does not
depend on θ and using > to indicate the transpose:
v > H = v > ∇2 L
= v > (∇> ∇L )
= ∇> (v > ∇L ).
In particular, (v > ∇L ) is a scalar, making this implementation well suited for auto-diff.

5. Case Studies
While the main contribution of this paper is to provide a
novel general approach for any order gradient estimation in
SCGs, we also provide a proof-of-concept empirical evaluation for a set of case studies, carried out on the iterated

DiCE: The Infinitely Differentiable Monte Carlo Estimator

Estimate

2

(b) Hessian

1
0

0

1
1

2
0 2 4 6 8
Gradient dimension

1

0.9

0.8

0.8

0.6

0.7

0.4

0.6
0.4

0

50

Hessian dimension

100

Figure 3. For each of the two agents (1 top row, 2 bottom row)
in the iterated prisoner’s dilemma, shown is the flattened true
(red) and estimated (green) Gradient (left) and Hessian (right)
using the first and second derivative of D I CE and the exact value
function respectively. The correlation coefficients are 0.999 for the
gradients and 0.97 for the Hessian; the sample size is 100k.

prisoner’s dilemma (IPD). In IPD, two agents iteratively
play matrix games with two possible actions: (C)ooperate
and (D)efect. The possible outcomes of each game are DD,
DC, CD, CC with the corresponding first agent payoffs, -2,
0, -3, -1, respectively. This setting is useful because (1) it
has a nontrivial but analytically calculable value function,
allowing for verification of gradient estimates, and (2) differentiating through the learning steps of other agents in
multi-agent RL is a highly relevant application of higher
order policy gradient estimators in RL (Foerster et al., 2018).
Empirical Verification. We first verify that D I CE successfully recovers gradients and Hessians in stochastic computation graphs. To do so, we use D I CE to estimate gradients
and Hessians of the expected return for fixed policies in
IPD.
As shown in Figure 3, we find that indeed the D I CE estimator matches both the gradients (a) and the Hessians (b) for
both agents accurately. Furthermore, Figure 4 shows how
the estimate of the gradient improve as the value function
becomes more accurate during training, in (a). Also shown
is the quality of the gradient estimation as a function of
sample size with and without a baseline, in (b). Both plots
show that the baseline is a key component of D I CE for
accurate estimation of gradients.
D I CE for multi-agent RL. In learning with opponentlearning awareness (LOLA), Foerster et al. (2018) show
that agents which differentiate through the learning step
of their opponent converge to Nash equilibria with higher
social welfare in the IPD.
Since the standard policy gradient learning step for one
agent has no dependency on the parameters of the other

(b)

1.0

0.2

0.5

0

0

(a)

1.0

Correlation

Estimate

(a) Gradient

1

2

3

Standard deviation of value

0.0

0

10000

20000

Sample size

Figure 4. Shown in (a) is the correlation of the gradient estimator
(averaged across agents) as a function of the estimation error of the
baseline when using a sample size of 128 and in (b) as a function
of sample size when using a converged baseline (in blue) and no
baseline (in green). In both plots errors bars indicate the standard
deviation.

agent (which it treats as part of the environment), LOLA
relies on a Taylor expansion of the expected return in combination with an analytical derivation of the second order
gradients to be able to differentiate through the expected
return after the opponent’s learning step.
Here, we take a more direct approach, made possible by
D I CE. Let πθ1 be the policy of the LOLA agent and let πθ2
be the policy of its opponent and vice versa. Assuming that
the opponent learns using policy gradients, LOLA-D I CE
agents learn by directly optimising the following stochastic
objective w.r.t. θ1 :
 
L1 (θ1 , θ2 )LOLA = Eπθ1 ,πθ2 +∆θ2 (θ1 ,θ2 ) L1 , where
(5.1)
 
∆θ2 (θ1 , θ2 ) = α∇θ2 Eπθ1 ,πθ2 L2 ,
PT
Here α is a scalar step size and Li = t=0 γ t rti is the sum
of discounted returns for agent i.
To evaluate these terms directly, our variant of LOLA unrolls the learning process of the opponent, which is functionally similar to model-agnostic meta-learning (MAML,
Finn et al., 2017). In the MAML formulation, the gradient
update of the opponent, ∆θ2 , corresponds to the inner loop
(typically training objective) and the gradient update of the
agent itself to the outer loop (typically test objective).
Using the following D I CE-objective to estimate gradient
steps for agent i, we are able preserve all dependencies:
X n t0 ≤t o
Li (θ1 ,θ2 ) =
aj∈{1,2} γ t rti ,
(5.2)
t

n 0
o
≤t
where atj∈{1,2}
is the set of all actions taken by both
agents up to time t. We note that for computational reasons,
we cache the ∆θi of the inner loop when unrolling the outer
loop policies in order to avoid recalculating them at every
time step.

Avg. per step return

DiCE: The Infinitely Differentiable Monte Carlo Estimator

1.00
1.25
1.50
1.75
2.00

(a) Naive Learner
no lookahead

0

50

100 150 200

1.00
1.25
1.50
1.75
2.00

(b) LOLA (original)

0

50

100 150 200
Iterations

1.00
1.25
1.50
1.75
2.00

(d) LOLA (with DiCE)

Lookahead

1 step
2 step
3 step

0

50

100 150 200

Figure 5. Joint average per step returns for different training methods. (a) Agents naively optimize expected returns w.r.t. their policy
parameters only, without lookahead steps. (b) The original LOLA algorithm (Foerster et al., 2018) that uses gradient corrections. (c)
LOLA-D I CE with lookahead of up to 3 gradient steps. Shaded areas represent the 95% confidence intervals based on 5 runs. All agents
used batches of size 64, which is more than 60 times smaller than the size required in the original LOLA paper.

Importantly, using D I CE, differentiating through ∆θ2 produces the correct higher order gradients, which is vital
for LOLA to function. In contrast, simply differentiating
through the SL-based first order gradient estimator again,
as was done for MAML by Finn et al. (2017), results in
omitted gradient terms and a biased gradient estimator, as
pointed out by Al-Shedivat et al. (2017) and Stadie et al.
(2018).
Figure 5 shows a comparison between the LOLA-D I CE
agents and the original formulation of LOLA. In our experiments, we use a time horizon of 150 steps and a reduced
batch size of 64; the lookahead gradient step, α, is set to 1
and the learning rate is 0.3. Importantly, given the approximation used, the LOLA method was restricted to a single
step of opponent learning. In contrast, using D I CE we
can differentiate through an arbitrary number of opponent
learning steps.
The original LOLA implemented via second order gradient
corrections shows no stable learning, as it requires much
larger batch sizes (∼ 4000). In contrast, LOLA-D I CE
agents discover strategies of high social welfare, replicating
the results of the original LOLA paper in a way that is both
more direct, efficient and establishes a common formulation
between MAML and LOLA.

6. Related Work
Gradient estimation is well studied, although many methods
have been named and explored independently in different
fields, and the primary focus has been on first order gradients. Fu (2006) provides an overview of methods from the
point of view of simulation optimization.
The score function (SF) estimator, also referred to as the
likelihood ratio estimator or REINFORCE, has received
considerable attention in many fields. In reinforcement

learning, policy gradient methods (Williams, 1992) have
proven highly successful, especially when combined with
variance reduction techniques (Weaver & Tao, 2001; Grondman et al., 2012). The SF estimator has also been used in the
analysis of stochastic systems (Glynn, 1990), as well as for
variational inference (Wingate & Weber, 2013; Ranganath
et al., 2014).
Kingma & Welling (2013) and Rezende et al. (2014) discuss Monte-Carlo gradient estimates in the case where the
stochastic parts of a model are amenable to reparameterisation.
To easily make use of these estimates for optimizing neural
network models, automatic differentiation for backpropagation has been widely used (Baydin et al., 2015).
These approaches are formalized for arbitrary computation
graphs by Schulman et al. (2015a), but to our knowledge our
paper is the first to present a practical and correct approach
for generating higher order gradient estimators utilizing
auto-diff.
One rapidly growing application area for such higher order
gradient estimates is meta-learning for reinforcement learning. Finn et al. (2017) compute a loss after a number policy
gradient learning steps, differentiating through the learning
step to find parameters that can be quickly fine-tuned for
different tasks. Li et al. (2017) extend this work to also
meta-learn the fine-tuning step direction and magnitude.
Al-Shedivat et al. (2017) and Stadie et al. (2018) derive the
proper higher order gradient estimators for their work by
reapplying the score function trick. Foerster et al. (2018)
use a multi-agent version of the same higher order gradient estimators in combination with a Taylor expansion of
the expected return. None present a general strategy for
constructing higher order gradient estimators for arbitrary
graphs.

DiCE: The Infinitely Differentiable Monte Carlo Estimator

7. Conclusion
We presented D I CE, a general method for computing any
order gradient estimators for stochastic computation graphs.
D I CE resolves the deficiencies of current approaches for
computing higher order gradient estimators: analytical calculation is error-prone and incompatible with auto-diff,
while repeated application of the Surrogate Loss approach is
cumbersome and, as we show, leads to incorrect estimators
in many cases. We prove the correctness of D I CE estimators, introduce a simple practical implementation of D I CE
for use in deep learning frameworks, and validate its correctness and utility in a multi-agent reinforcement learning
problem. We believe D I CE will unlock further exploration
and adoption of higher order learning methods in metalearning, reinforcement learning, and other applications of
stochastic computation graphs.

Acknowledgements
This project has received funding from the European Research Council (ERC) under the European Unions Horizon
2020 research and innovation programme (grant agreement
number 637713). It was also supported by the OxfordGoogle DeepMind Graduate Scholarship, the UK EPSRC
CDT in Autonomous Intelligent Machines and Systems.
We would like to thank Misha Denil, Brendan Shillingford and Wendelin Boehmer for providing feedback on the
manuscript.

Finn, Chelsea, Abbeel, Pieter, and Levine, Sergey. Modelagnostic meta-learning for fast adaptation of deep networks. In Proceedings of the 34th International Conference on Machine Learning, ICML 2017, Sydney, NSW,
Australia, 6-11 August 2017, pp. 1126–1135, 2017.
Foerster, Jakob N., Chen, Richard Y., Al-Shedivat, Maruan,
Whiteson, Shimon, Abbeel, Pieter, and Mordatch, Igor.
Learning with opponent-learning awareness. In AAMAS,
2018.
Fu, Michael C. Gradient estimation. Handbooks in operations research and management science, 13:575–616,
2006.
Glynn, Peter W. Likelihood ratio gradient estimation for
stochastic systems. Communications of the ACM, 33(10):
75–84, 1990.
Grondman, Ivo, Busoniu, Lucian, Lopes, Gabriel AD, and
Babuska, Robert. A survey of actor-critic reinforcement
learning: Standard and natural policy gradients. IEEE
Transactions on Systems, Man, and Cybernetics, Part C
(Applications and Reviews), 42(6):1291–1307, 2012.
Kingma, Diederik P. and Welling, Max. Auto-encoding
variational bayes. CoRR, abs/1312.6114, 2013.
Li, Zhenguo, Zhou, Fengwei, Chen, Fei, and Li, Hang.
Meta-sgd: Learning to learn quickly for few shot learning.
CoRR, abs/1707.09835, 2017.

References
Abadi, Martı́n, Barham, Paul, Chen, Jianmin, Chen, Zhifeng,
Davis, Andy, Dean, Jeffrey, Devin, Matthieu, Ghemawat,
Sanjay, Irving, Geoffrey, Isard, Michael, Kudlur, Manjunath, Levenberg, Josh, Monga, Rajat, Moore, Sherry,
Murray, Derek Gordon, Steiner, Benoit, Tucker, Paul A.,
Vasudevan, Vijay, Warden, Pete, Wicke, Martin, Yu, Yuan,
and Zheng, Xiaoqiang. Tensorflow: A system for largescale machine learning. In 12th USENIX Symposium on
Operating Systems Design and Implementation, OSDI
2016, Savannah, GA, USA, November 2-4, 2016., pp.
265–283, 2016.
Al-Shedivat, Maruan, Bansal, Trapit, Burda, Yuri, Sutskever,
Ilya, Mordatch, Igor, and Abbeel, Pieter. Continuous
adaptation via meta-learning in nonstationary and competitive environments. CoRR, abs/1710.03641, 2017.
Baydin, Atilim Gunes, Pearlmutter, Barak A., and Radul,
Alexey Andreyevich. Automatic differentiation in machine learning: a survey. CoRR, abs/1502.05767, 2015.
Dennis, Jr, John E and Moré, Jorge J. Quasi-newton methods, motivation and theory. SIAM review, 19(1):46–89,
1977.

Paszke, Adam, Gross, Sam, Chintala, Soumith, Chanan,
Gregory, Yang, Edward, DeVito, Zachary, Lin, Zeming,
Desmaison, Alban, Antiga, Luca, and Lerer, Adam. Automatic differentiation in pytorch. 2017.
Pearlmutter, Barak A. Fast exact multiplication by the
hessian. Neural computation, 6(1):147–160, 1994.
Ranganath, Rajesh, Gerrish, Sean, and Blei, David M. Black
box variational inference. In Proceedings of the Seventeenth International Conference on Artificial Intelligence
and Statistics, AISTATS 2014, Reykjavik, Iceland, April
22-25, 2014, pp. 814–822, 2014.
Rezende, Danilo Jimenez, Mohamed, Shakir, and Wierstra,
Daan. Stochastic backpropagation and approximate inference in deep generative models. pp. 1278–1286, 2014.
Schulman, John, Heess, Nicolas, Weber, Theophane, and
Abbeel, Pieter. Gradient Estimation Using Stochastic
Computation Graphs. In Advances in Neural Information
Processing Systems 28: Annual Conference on Neural
Information Processing Systems 2015, December 7-12,
2015, Montreal, Quebec, Canada, pp. 3528–3536, 2015a.

DiCE: The Infinitely Differentiable Monte Carlo Estimator

Schulman, John, Levine, Sergey, Abbeel, Pieter, Jordan,
Michael, and Moritz, Philipp. Trust region policy optimization. In International Conference on Machine Learning, pp. 1889–1897, 2015b.
Schulman, John, Abbeel, Pieter, and Chen, Xi. Equivalence
between policy gradients and soft q-learning. CoRR,
abs/1704.06440, 2017.
Stadie, Bradly, Yang, Ge, Houthooft, Rein, Chen, Xi,
Duan, Yan, Wu, Yuhuai, Abbeel, Pieter, and Sutskever,
Ilya. Some considerations on learning to explore via
meta-reinforcement learning, 2018. URL https://
openreview.net/forum?id=Skk3Jm96W.
Weaver, Lex and Tao, Nigel. The optimal reward baseline
for gradient-based reinforcement learning. In Proceedings of the Seventeenth conference on Uncertainty in
artificial intelligence, pp. 538–545. Morgan Kaufmann
Publishers Inc., 2001.
Williams, Ronald J. Simple statistical gradient-following
algorithms for connectionist reinforcement learning. In
Reinforcement Learning, pp. 5–32. Springer, 1992.
Wingate, David and Weber, Theophane. Automated variational inference in probabilistic programming. CoRR,
abs/1301.1299, 2013.

