Year: 2020

Abstract:
We derive an unbiased estimator for expectations over discrete random variables based on sampling without replacement, which reduces variance as it avoids duplicate samples. We show that our estimator can be derived as the RaoBlackwellization of three different estimators. Combining our estimator with REINFORCE, we obtain a policy gradient estimator and we reduce its variance using a built-in control variate which is obtained without additional model evaluations. The resulting estimator is closely related to other gradient estimators. Experiments with a toy problem, a categorical Variational Auto-Encoder and a structured prediction problem show that our estimator is the only estimator that is consistently among the best estimators in both high and low entropy settings.

heading: 1 INTRODUCTION
text: Put replacement in your basement! We derive the unordered set estimator1: an unbiased (gradient) estimator for expectations over discrete random variables based on (unordered sets of) samples without replacement. In particular, we consider the problem of estimating (the gradient of) the expectation of f(x) where x has a discrete distribution p over the domain D, i.e.
Ex∼p(x)[f(x)] = ∑
x∈D p(x)f(x). (1)
This expectation comes up in reinforcement learning, discrete latent variable modelling (e.g. for compression), structured prediction (e.g. for translation), hard attention and many other tasks that use models with discrete operations in their computational graphs (see e.g. Jang et al. (2016)). In general, x has structure (such as a sequence), but we can treat it as a ‘flat’ distribution, omitting the bold notation, so x has a categorical distribution over D given by p(x), x ∈ D. Typically, the distribution has parameters θ, which are learnt through gradient descent. This requires estimating the gradient ∇θEx∼pθ(x)[f(x)], using a set of samples S. A gradient estimate e(S) is unbiased if
ES [e(S)] = ∇θEx∼pθ(x)[f(x)]. (2)
The samples S can be sampled independently or using alternatives such as stratified sampling which reduce variance to increase the speed of learning. In this paper, we derive an unbiased gradient estimator that reduces variance by avoiding duplicate samples, i.e. by sampling S without replacement. This is challenging as samples without replacement are dependent and have marginal distributions that are different from p(x). We further reduce the variance by deriving a built-in control variate, which maintains the unbiasedness and does not require additional samples.
Related work. Many algorithms for estimating gradients for discrete distributions have been proposed. A general and widely used estimator is REINFORCE (Williams, 1992). Biased gradients based on a continuous relaxations of the discrete distribution (known as Gumbel-Softmax or Concrete) were jointly introduced by Jang et al. (2016) and Maddison et al. (2016). These can be combined with the straight through estimator (Bengio et al., 2013) if the model requires discrete samples or be used to construct control variates for REINFORCE, as in REBAR (Tucker et al., 2017) or
1Code available at https://github.com/wouterkool/estimating-gradients-without-replacement.
RELAX (Grathwohl et al., 2018). Many other methods use control variates and other techniques to reduce the variance of REINFORCE (Paisley et al., 2012; Ranganath et al., 2014; Gregor et al., 2014; Mnih & Gregor, 2014; Gu et al., 2016; Mnih & Rezende, 2016).
Some works rely on explicit summation of the expectation, either for the marginal distribution (Titsias & Lázaro-Gredilla, 2015) or globally summing some categories while sampling from the remainder (Liang et al., 2018; Liu et al., 2019). Other approaches use a finite difference approximation to the gradient (Lorberbom et al., 2018; 2019). Yin et al. (2019) introduced ARSM, which uses multiple model evaluations where the number adapts automatically to the uncertainty.
In the structured prediction setting, there are many algorithms for optimizing a quantity under a sequence of discrete decisions, using (weak) supervision, multiple samples (or deterministic model evaluations), or a combination both (Ranzato et al., 2016; Shen et al., 2016; He et al., 2016; Norouzi et al., 2016; Bahdanau et al., 2017; Edunov et al., 2018; Leblond et al., 2018; Negrinho et al., 2018). Most of these algorithms are biased and rely on pretraining using maximum likelihood or gradually transitioning from supervised to reinforcement learning. Using Gumbel-Softmax based approaches in a sequential setting is difficult as the bias accumulates because of mixing errors (Gu et al., 2018).

heading: 2 PRELIMINARIES
text: Throughout this paper, we will denote with Bk an ordered sample without replacement of size k and with Sk an unordered sample (of size k) from the categorical distribution p.
Restricted distribution. When sampling without replacement, we remove the set C ⊂ D already sampled from the domain and we denote with pD\C the distribution restricted to the domain D \C:
pD\C(x) = p(x) 1− ∑ c∈C p(c) , x ∈ D \ C. (3)
Ordered sample without replacement Bk. Let Bk = (b1, ..., bk), bi ∈ D be an ordered sample without replacement, which is generated from the distribution p as follows: first, sample b1 ∼ p, then sample b2 ∼ pD\{b1}, b3 ∼ pD\{b1,b2}, etc. i.e. elements are sampled one by one without replacement. Using this procedure, Bk can be seen as a (partial) ranking according to the PlackettLuce model (Plackett, 1975; Luce, 1959) and the probability of obtaining the vector Bk is
p(Bk) = k∏ i=1 pD\B i−1 (bi) = k∏ i=1
p(bi) 1− ∑ j<i p(bj) . (4)
We can also restrict Bk to the domain D \ C, which means that bi 6∈ C for i = 1, ..., k:
pD\C(Bk) = k∏ i=1
pD\C(bi) 1− ∑ j<i pD\C(bj) = k∏ i=1
p(bi) 1− ∑ c∈C p(c)− ∑ j<i p(bj) . (5)
Unordered sample without replacement. Let Sk ⊆ D be an unordered sample without replacement from the distribution p, which can be generated simply by generating an ordered sample and discarding the order. We denote elements in the sample with s ∈ Sk (so without index) and we write B(Sk) as the set of all k! permutations (orderings)Bk that correspond to (could have generated) Sk. It follows that the probability for sampling Sk is given by:
p(Sk) = ∑
Bk∈B(Sk)
p(Bk) = ∑
Bk∈B(Sk) k∏ i=1
p(bi) 1− ∑ j<i p(bj) = ∏ s∈Sk p(s)  · ∑ Bk∈B(Sk) k∏ i=1
1 1− ∑ j<i p(bj) .
(6) The last step follows since Bk ∈ B(Sk) is an ordering of Sk, such that ∏k i=1 p(bi) = ∏ s∈S p(s). Naive computation of p(Sk) is O(k!), but in Appendix B we show how to compute it efficiently.
When sampling from the distribution restricted to D \ C, we sample Sk ⊆ D \ C with probability:
pD\C(Sk) = ∏ s∈Sk p(s)  · ∑ Bk∈B(Sk) k∏ i=1
1 1− ∑ c∈C p(c)− ∑ j<i p(bj) . (7)
The Gumbel-Top-k trick. As an alternative to sequential sampling, we can also sample Bk and Sk by taking the top k of Gumbel variables (Yellott, 1977; Vieira, 2014; Kim et al., 2016). Following notation from Kool et al. (2019c), we define the perturbed log-probability gφi = φi + gi, where φi = log p(i) and gi ∼ Gumbel(0). Then let b1 = arg maxi∈D gφi , b2 = arg maxi∈D\{b1} gφi , etc., soBk is the top k of the perturbed log-probabilities in decreasing order. The probability of obtaining Bk using this procedure is given by equation 4, so this provides an alternative sampling method which is effectively a (non-differentiable) reparameterization of sampling without replacement. For a differentiable reparameterization, see Grover et al. (2019).
It follows that taking the top k perturbed log-probabilities without order, we obtain the unordered sample set Sk. This way of sampling underlies the efficient computation of p(Sk) in Appendix B.

heading: 3 METHODOLOGY
text: In this section, we derive the unordered set policy gradient estimator: a low-variance, unbiased estimator of ∇θEpθ(x)[f(x)] based on an unordered sample without replacement Sk. First, we derive the generic (non-gradient) estimator for E[f(x)] as the Rao-Blackwellized version of a single sample Monte Carlo estimator (and two other estimators!). Then we combine this estimator with REINFORCE (Williams, 1992) and we show how to reduce its variance using a built-in baseline.

heading: 3.1 RAO-BLACKWELLIZATION OF THE SINGLE SAMPLE ESTIMATOR
text: A very crude but simple estimator for E[f(x)] based on the ordered sample Bk is to only use the first element b1, which by definition is a sample from the distribution p. We define this estimator as the single sample estimator, which is unbiased, since
EBk∼p(Bk)[f(b1)] = Eb1∼p(b1)[f(b1)] = Ex∼p(x)[f(x)]. (8)
Discarding all but one sample, the single sample estimator is inefficient, but we can use RaoBlackwellization (Casella & Robert, 1996) to signficantly improve it. To this end, we consider the distribution Bk|Sk, which is, knowing the unordered sample Sk, the conditional distribution over ordered samples Bk ∈ B(Sk) that could have generated Sk.2 Using Bk|Sk, we rewrite E[f(b1)] as
EBk∼p(Bk)[f(b1)] = ESk∼p(Sk) [ EBk∼p(Bk|Sk) [f(b1)] ] = ESk∼p(Sk) [ Eb1∼p(b1|Sk) [f(b1)] ] .
The Rao-Blackwellized version of the single sample estimator computes the inner conditional expectation exactly. Since Bk is an ordering of Sk, we have b1 ∈ Sk and we can compute this as
Eb1∼p(b1|Sk) [f(b1)] = ∑ s∈Sk P (b1 = s|Sk)f(s) (9)
where, in a slight abuse of notation, P (b1 = s|Sk) is the probability that the first sampled element b1 takes the value s, given that the complete set of k samples is Sk. Using Bayes’ Theorem we find
P (b1 = s|Sk) = p(Sk|b1 = s)P (b1 = s) p(Sk) = pD\{s}(Sk \ {s})p(s) p(Sk) . (10)
The step p(Sk|b1 = s) = pD\{s}(Sk \ {s}) comes from analyzing sequential sampling without replacement: given that the first element sampled is s, the remaining elements have a distribution restricted toD\{s}, so sampling Sk (including s) given the first element s is equivalent to sampling the remainder Sk \{s} from the restricted distribution, which has probability pD\{s}(Sk \{s}) (see equation 7).
2Note that Bk|Sk is not a Plackett-Luce distribution restricted to Sk!
The unordered set estimator. For notational convenience, we introduce the leave-one-out ratio.
Definition 1. The leave-one-out ratio of s w.r.t. the set S is given by R(Sk, s) = p D\{s}(Sk\{s})
p(Sk) .
Rewriting equation 10 as P (b1 = s|Sk) = p(s)R(Sk, s) shows that the probability of sampling s first, given Sk, is simply the unconditional probability multiplied by the leave-one-out ratio. We now define the unordered set estimator as the Rao-Blackwellized version of the single-sample estimator. Theorem 1. The unordered set estimator, given by
eUS(Sk) = ∑ s∈Sk p(s)R(Sk, s)f(s) (11)
is the Rao-Blackwellized version of the (unbiased!) single sample estimator.
Proof. Using P (b1 = s|Sk) = p(s)R(Sk, s) in equation 9 we have Eb1∼p(b1|Sk) [f(b1)] = ∑ s∈Sk P (b1 = s|Sk)f(s) = ∑ s∈Sk p(s)R(Sk, s)f(s). (12)
The implication of this theorem is that the unordered set estimator, in explicit form given by equation 11, is an unbiased estimator of E[f(x)] since it is the Rao-Blackwellized version of the unbiased single sample estimator. Also, as expected by taking multiple samples, it has variance equal or lower than the single sample estimator by the Rao-Blackwell Theorem (Lehmann & Scheffé, 1950).

heading: 3.2 RAO-BLACKWELLIZATION OF OTHER ESTIMATORS
text: The unordered set estimator is also the result of Rao-Blackwellizing two other unbiased estimators: the stochastic sum-and-sample estimator and the importance-weighted estimator.
The sum-and-sample estimator. We define as sum-and-sample estimator any estimator that relies on the identity that for any C ⊂ D
Ex∼p(x)[f(x)] = Ex∼pD\C(x) [∑ c∈C p(c)f(c) + ( 1− ∑ c∈C p(c) ) f(x) ] . (13)
For the derivation, see Appendix C.1 or Liang et al. (2018); Liu et al. (2019). In general, a sum-andsample estimator with a budget of k > 1 evaluations sums expectation terms for a set of categoriesC (s.t. |C| < k) explicitly (e.g. selected by their value f (Liang et al., 2018) or probability p (Liu et al., 2019)), and uses k− |C| (down-weighted) samples from D \C to estimate the remaining terms. As is noted by Liu et al. (2019), selecting C such that 1− ∑ c∈C p(c)
k−|C| is minimized guarantees to reduce variance compared to a standard minibatch of k samples (which is equivalent to setting C = ∅). See also Fearnhead & Clifford (2003) for a discussion on selecting C optimally. The ability to optimize C depends on whether p(c) can be computed efficiently a-priori (before sampling). This is difficult in high-dimensional settings, e.g. sequence models which compute the probability incrementally while ancestral sampling. An alternative is to select C stochastically (as equation 13 holds for any C), and we choose C = Bk−1 to define the stochastic sum-and-sample estimator:
eSSAS(Bk) = k−1∑ j=1 p(bj)f(bj) + 1− k−1∑ j=1 p(bj)  f(bk). (14) For simplicity, we consider the version that sums k − 1 terms here, but the following results also hold for a version that sums k−m terms and uses m samples (without replacement) (see Appendix C.3). Sampling without replacement, it holds that bk|Bk−1 ∼ pD\B k−1 , so the unbiasedness follows from equation 13 by separating the expectation over Bk into expectations over Bk−1 and bk|Bk−1:
EBk−1∼p(Bk−1) [ Ebk∼p(bk|Bk−1) [ eSSAS(Bk) ]] = EBk−1∼p(Bk−1) [E[f(x)]] = E[f(x)].
In general, a sum-and-sample estimator reduces variance if the probability mass is concentrated on the summed categories. As typically high probability categories are sampled first, the stochastic sum-and-sample estimator sums high probability categories, similar to the estimator by Liu et al. (2019) which we refer to as the deterministic sum-and-sample estimator. As we show in Appendix C.2, Rao-Blackwellizing the stochastic sum-and-sample estimator also results in the unordered set estimator. This even holds for a version that usesm samples and k−m summed terms (see Appendix C.3), which means that the unordered set estimator has equal or lower variance than the optimal (in terms of m) stochastic sum-and-sample estimator, but conveniently does not need to choose m.
The importance-weighted estimator. The importance-weighted estimator (Vieira, 2017) is
eIW(Sk, κ) = ∑ s∈Sk p(s) q(s, κ) f(s). (15)
This estimator is based on the idea of priority sampling (Duffield et al., 2007). It does not use the order of the sample, but assumes sampling using the Gumbel-Top-k trick and requires access to κ, the (k + 1)-th largest perturbed log-probability, which can be seen as the ‘threshold’ since gφs > κ ∀s ∈ Sk. q(s, a) = P (gφs > a) can be interpreted as the inclusion probability of s ∈ Sk (assuming a fixed threshold a instead of a fixed sample size k). For details and a proof of unbiasedness, see Vieira (2017) or Kool et al. (2019c). As the estimator has high variance, Kool et al. (2019c) resort to normalizing the importance weights, resulting in biased estimates. Instead, we use Rao-Blackwellization to eliminate stochasticity by κ. Again, the result is the unordered set estimator (see Appendix D.1), which thus has equal or lower variance.

heading: 3.3 THE UNORDERED SET POLICY GRADIENT ESTIMATOR
text: Writing pθ to indicate the dependency on the model parameters θ, we can combine the unordered set estimator with REINFORCE (Williams, 1992) to obtain the unordered set policy gradient estimator.
Corollary 1. The unordered set policy gradient estimator, given by eUSPG(Sk) = ∑ s∈Sk pθ(s)R(S k, s)∇θ log pθ(s)f(s) = ∑ s∈Sk ∇θpθ(s)R(Sk, s)f(s), (16)
is an unbiased estimate of the policy gradient.
Proof. Using REINFORCE (Williams, 1992) combined with the unordered set estimator we find:
∇θEpθ(x)[f(x)]=Epθ(x)[∇θ log pθ(x)f(x)]=ESk∼pθ(Sk) ∑ s∈Sk pθ(s)R(S k, s)∇θ log pθ(s)f(s) .
Variance reduction using a built-in control variate. The variance of REINFORCE can be reduced by subtracting a baseline from f . When taking multiple samples (with replacement), a simple and effective baseline is to take the mean of other (independent!) samples (Mnih & Rezende, 2016). Sampling without replacement, we can use the same idea to construct a baseline based on the other samples, but we have to correct for the fact that the samples are not independent.
Theorem 2. The unordered set policy gradient estimator with baseline, given by
eUSPGBL(Sk) = ∑ s∈Sk ∇θpθ(s)R(Sk, s) f(s)− ∑ s′∈Sk pθ(s ′)RD\{s}(Sk, s′)f(s′)  , (17) where
RD\{s}(Sk, s′) = p D\{s,s′} θ (S k \ {s, s′}) p D\{s} θ (S k \ {s}) (18)
is the second order leave-one-out ratio, is an unbiased estimate of the policy gradient.
Proof. See Appendix E.1.
This theorem shows how to include a built-in baseline based on dependent samples (without replacement), without introducing bias. By having a built-in baseline, the value f(s) for sample s is compared against an estimate of its expectation E[f(s)], based on the other samples. The difference is an estimate of the advantage (Sutton & Barto, 2018), which is positive if the sample s is ‘better’ than average, causing pθ(s) to be increased (reinforced) through the sign of the gradient, and vice versa. By sampling without replacement, the unordered set estimator forces the estimator to compare different alternatives, and reinforces the best among them.
Including the pathwise derivative. So far, we have only considered the scenario where f does not depend on θ. If f does depend on θ, for example in a VAE (Kingma & Welling, 2014; Rezende et al., 2014), then we use the notation fθ and we can write the gradient (Schulman et al., 2015) as
∇θEpθ(x)[fθ(x)] = Epθ(x)[∇θ log pθ(x)fθ(x) +∇θfθ(x)]. (19)
The additional second (‘pathwise’) term can be estimated (using the same samples) with the standard unordered set estimator. This results in the full unordered set policy gradient estimator:
eFUSPG(Sk) = ∑ s∈Sk ∇θpθ(s)R(Sk, s)fθ(s) + ∑ s∈Sk pθ(s)R(S k, s)∇θfθ(s)
= ∑ s∈Sk R(Sk, s)∇θ (pθ(s)fθ(s)) (20)
Equation 20 is straightforward to implement using an automatic differentiation library. We can also include the baseline (as in equation 17) but we must make sure to call STOP GRADIENT (DETACH in PyTorch) on the baseline (but not on fθ(s)!). Importantly, we should never track gradients through the leave-one-out ratioR(Sk, s) which means it can be efficiently computed in pure inference mode.
Scope & limitations. We can use the unordered set estimator for any discrete distribution from which we can sample without replacement, by treating it as a univariate categorical distribution over its domain. This includes sequence models, from which we can sample using Stochastic Beam Search (Kool et al., 2019c), as well as multivariate categorical distributions which can also be treated as sequence models (see Section 4.2). In the presence of continuous variables or a stochastic function f , we may separate this stochasticity from the stochasticity over the discrete distribution, as in Lorberbom et al. (2019). The computation of the leave-one-out ratios adds some overhead, although they can be computed efficiently, even for large k (see Appendix B). For a moderately sized model, the costs of model evaluation and backpropagation dominate the cost of computing the estimator.

heading: 3.4 RELATION TO OTHER MULTI-SAMPLE ESTIMATORS
text: Relation to Murthy’s estimator. We found out that the ‘vanilla’ unordered set estimator (equation 11) is actually a special case of the estimator by Murthy (1957), known in statistics literature for estimation of a population total Θ = ∑ i∈D yi. Using yi = p(i)f(i), we have Θ = E[f(i)], so Murthy’s estimator can be used to estimate expectations (see equation 11). Murthy derives the estimator by ‘unordering’ a convex combination of Raj (1956) estimators, which, using yi = p(i)f(i), are stochastic sum-and-sample estimators in our analogy.
Murthy (1957) also provides an unbiased estimator of the variance, which may be interesting for future applications. Since Murthy’s estimator can be used with arbitrary sampling distribution, it is straightforward to derive importance-sampling versions of our estimators. In particular, we can sample S without replacement using q(x) > 0, x ∈ D, and use equations 11, 16, 17 and 20, as long as we compute the leave-one-out ratio R(Sk, s) using q.
While part of our derivation coincides with Murthy (1957), we are not aware of previous work using this estimator to estimate expectations. Additionally, we discuss practical computation of p(S) (Appendix B), we show the relation to the importance-weighted estimator, and we provide the extension to estimating policy gradients, especially including a built-in baseline without adding bias.
Relation to the empirical risk estimator. The empirical risk loss (Edunov et al., 2018) estimates the expectation in equation 1 by summing only a subset S of the domain, using normalized probabilities p̂θ(s) =
pθ(s)∑ s′∈S pθ(s) . Using this loss, the (biased) estimate of the gradient is given by
eRISK(Sk) = ∑ s∈Sk ∇θ ( pθ(s)∑ s′∈Sk pθ(s ′) ) f(s). (21)
The risk estimator is similar to the unordered set policy gradient estimator, with two important differences: 1) the individual terms are normalized by the total probability mass rather than the leave-one-out ratio and 2) the gradient w.r.t. the normalization factor is taken into account. As a result, samples ‘compete’ for probability mass and only the best can be reinforced. This has the same effect as using a built-in baseline, which we prove in the following theorem. Theorem 3. By taking the gradient w.r.t. the normalization factor into account, the risk estimator has a built-in baseline, which means it can be written as
eRISK(Sk) = ∑ s∈Sk ∇θpθ(s) 1∑ s′′∈Sk pθ(s ′′) f(s)− ∑ s′∈Sk pθ(s ′) 1∑ s′′∈Sk pθ(s ′′) f(s′)  . (22) Proof. See Appendix F.1
This theorem highlights the similarity between the biased risk estimator and our unbiased estimator (equation 17), and suggests that their only difference is the weighting of terms. Unfortunately, the implementation by Edunov et al. (2018) has more sources of bias (e.g. length normalization), which are not compatible with our estimator. However, we believe that our analysis helps analyze the bias of the risk estimator and is a step towards developing unbiased estimators for structured prediction.
Relation to VIMCO. VIMCO (Mnih & Rezende, 2016) is an estimator that uses k samples (with replacement) to optimize an objective of the form log 1k ∑ i f(xi), which is a multi-sample stochastic lower bound in the context of variational inference. VIMCO reduces the variance by using a local baseline for each of the k samples, based on the other k − 1 samples. While we do not have a log term, as our goal is to optimize general E[f(x)], we adopt the idea of forming a baseline based on the other samples, and we define REINFORCE with replacement (with built-in baseline) as the estimator that computes the gradient estimate using samples with replacement Xk = (x1, ..., xk) as
eRFWR(Xk) = 1
k k∑ i=1 ∇θ log pθ(xi) f(xi)− 1 k − 1 ∑ j 6=i f(xj)  . (23) This estimator is unbiased, as Exi,xj [∇θ log pθ(xi)f(xj)] = 0 for i 6= j (see also Kool et al. (2019b)). We think of the unordered set estimator as the without-replacement version of this estimator, which weights terms by pθ(s)R(Sk, s) instead of 1k . This puts more weight on higher probability elements to compensate for sampling without replacement. If probabilities are small and (close to) uniform, there are (almost) no duplicate samples and the weights will be close to 1k , so the gradient estimate of the with- and without-replacement versions are similar.
Relation to ARSM. ARSM (Yin et al., 2019) also uses multiple evaluations (‘pseudo-samples’) of pθ and f . This can be seen as similar to sampling without replacement, and the estimator also has a built-in control variate. Compared to ARSM, our estimator allows direct control over the computational cost (through the sample size k) and has wider applicability, for example it also applies to multivariate categorical variables with different numbers of categories per dimension.
Relation to stratified/systematic sampling. Our estimator aims to reduce variance by changing the sampling distribution for multiple samples by sampling without replacement. There are alternatives, such as using stratified or systematic sampling (see, e.g. Douc & Cappé (2005)). Both partition the domain D into k strata and take a single sample from each stratum, where systematic sampling uses common random numbers for each stratum. In applications involving high-dimensional or structured domains, it is unclear how to partition the domain and how to sample from each partition. Additionally, as samples are not independent, it is non-trivial to include a built-in baseline, which we find is a key component that makes our estimator perform well.

heading: 4 EXPERIMENTS

heading: 4.1 BERNOULLI TOY EXPERIMENT
text: We use the code by Liu et al. (2019) to reproduce their Bernoulli toy experiment. Given a vector p = (0.6, 0.51, 0.48) the goal is to minimize the loss L(η) = Ex1,x2,x3∼Bern(σ(η)) [∑3 i=1(xi − pi)2 ] . Here x1, x2, x3 are i.i.d. from the Bernoulli(σ(η)) distribution, parameterized by a scalar η ∈ R, where σ(η) = (1 + exp(−η))−1 is the sigmoid function. We compare different estimators, with and without baseline (either ‘built-in’ or using additional samples, referred to as REINFORCE+ in Liu et al. (2019)). We report the (log-)variance of the scalar gradient ∂L∂η as a function of the number of model evaluations, which is twice as high when using a sampled baseline (for each term).
As can be seen in Figure 1, the unordered set estimator is the only estimator that has consistently the lowest (or comparable) variance in both the high (η = 0) and low entropy (η = −4) regimes and for different number of samples/model evaluations. This suggests that it combines the advantages of the other estimators. We also ran the actual optimization experiment, where with as few as k = 3 samples the trajectory was indistinguishable from using the exact gradient (see Liu et al. (2019)).

heading: 4.2 CATEGORICAL VARIATIONAL AUTO-ENCODER
text: We use the code from Yin et al. (2019) to train a categorical Variational Auto-Encoder (VAE) with 20 dimensional latent space, with 10 categories per dimension (details in Appendix G.1). To use our estimator, we treat this as a single factorized distribution with 1020 categories from which we can sample without replacement using Stochastic Beam Search (Kool et al., 2019c), sequentially sampling each dimension as if it were a sequence model. We also perform experiments with 102 latent space, which provides a lower entropy setting, to highlight the advantage of our estimator.
Measuring the variance. In Table 1, we report the variance of different gradient estimators with k = 4 samples, evaluated on a trained model. The unordered set estimator has the lowest variance in both the small and large domain (low and high entropy) setting, being on-par with the best of the (stochastic3) sum-and-sample estimator and REINFORCE with replacement4. This confirms the toy experiment, suggesting that the unordered set estimator provides the best of both estimators. In Appendix G.2 we repeat the same experiment at different stages of training, with similar results.
3We cannot use the deterministic version by Liu et al. (2019) since we cannot select the top k categories. 4We cannot compare against VIMCO (Mnih & Rezende, 2016) as it optimizes a different objective.
ELBO optimization. We use different estimators to optimize the ELBO (details in Appendix G.1). Additionally to the baselines by Yin et al. (2019) we compare against REINFORCE with replacement and the stochastic sum-and-sample estimator. In Figure 2 we observe that our estimator performs on par with REINFORCE with replacement (and built-in baseline, equation 23) and outperforms other estimators in at least one of the settings. There are a lot of other factors, e.g. exploration that may explain why we do not get a strictly better result despite the lower variance. We note some overfitting (see validation curves in Appendix G.2), but since our goal is to show improved optimization, and to keep results directly comparable to Yin et al. (2019), we consider regularization a separate issue outside the scope of this paper. These results are using MNIST binarized by a threshold of 0.5. In Appendix G.2 we report results using the standard binarized MNIST dataset from Salakhutdinov & Murray (2008).

heading: 4.3 STRUCTURED PREDICTION FOR THE TRAVELLING SALESMAN PROBLEM
text: To show the wide applicability of our estimator, we consider the structured prediction task of predicting routes (sequences) for the Travelling Salesman Problem (TSP) (Vinyals et al., 2015; Bello et al., 2016; Kool et al., 2019a). We use the code by Kool et al. (2019a)5 to reproduce their TSP experiment with 20 nodes. For details, see Appendix H.
We implement REINFORCE with replacement (and built-in baseline) as well as the stochastic sumand-sample estimator and our estimator, using Stochastic Beam Search (Kool et al., 2019c) for sampling. Also, we include results using the biased normalized importance-weighted policy gradient estimator with built-in baseline (derived in Kool et al. (2019b), see Appendix D.2). Additionally, we compare against REINFORCE with greedy rollout baseline (Rennie et al., 2017) used by Kool et al. (2019c) and a batch-average baseline. For reference, we also include the biased risk estimator, either ‘sampling’ using stochastic or deterministic beam search (as in Edunov et al. (2018)).
In Figure 3a, we compare training progress (measured on the validation set) as a function of the number of training steps, where we divide the batch size by k to keep the total number of samples equal. Our estimator outperforms REINFORCE with replacement, the stochastic sum-and-sample estimator and the strong greedy rollout baseline (which uses additional baseline model evaluations) and performs on-par with the biased risk estimator. In Figure 3b, we plot the same results against the number of instances, which shows that, compared to the single sample estimators, we can train with less data and less computational cost (as we only need to run the encoder once for each instance).
5https://github.com/wouterkool/attention-learn-to-route

heading: 5 DISCUSSION
text: We introduced the unordered set estimator, a low-variance, unbiased gradient estimator based on sampling without replacement, which can be used as an alternative to the popular biased GumbelSoftmax estimator (Jang et al., 2016; Maddison et al., 2016). Our estimator is the result of RaoBlackwellizing three existing estimators, which guarantees equal or lower variance, and is closely related to a number of other estimators. It has wide applicability, is parameter free (except for the sample size k) and has competitive performance to the best of alternatives in both high and low entropy regimes.
In our experiments, we found that REINFORCE with replacement, with multiple samples and a built-in baseline as inspired by VIMCO (Mnih & Rezende, 2016), is a simple yet strong estimator which has performance similar to our estimator in the high entropy setting. We are not aware of any recent work on gradient estimators for discrete distributions that has considered this estimator as baseline, while it may be often preferred given its simplicity. In future work, we want to investigate if we can apply our estimator to estimate gradients ‘locally’ (Titsias & Lázaro-Gredilla, 2015), as locally we have a smaller domain and expect more duplicate samples.

heading: ACKNOWLEDGMENTS
text: This research was funded by ORTEC. We would like to thank anonymous reviewers for their feedback that helped improve the paper.

heading: A NOTATION
text: Throughout this appendix we will use the following notation from Maddison et al. (2014):
eφ(g) = exp(−g + φ) Fφ(g) = exp(− exp(−g + φ)) fφ(g) = eφ(g)Fφ(g).
This means that Fφ(g) is the CDF and fφ(g) the PDF of the Gumbel(φ) distribution. Additionally we will use the identities by Maddison et al. (2014):
Fφ(g)Fγ(g) = Flog(exp(φ)+exp(γ))(g) (24)∫ b g=a eγ(g)Fφ(g)∂g = (Fφ(b)− Fφ(a)) exp(γ) exp(φ) . (25)
Also, we will use the following notation, definitions and identities (see Kool et al. (2019c)):
φi = log p(i) (26) φS = log ∑ i∈S p(i) = log ∑ i∈S expφi (27)
φD\S = log ∑ i∈D\S p(i) = log
( 1−
∑ i∈S p(i)
) = log(1− exp(φS)) (28)
Gφi ∼ Gumbel(φi) (29) GφS = max
i∈S Gφi ∼ Gumbel(φS) (30)
For a proof of equation 30, see Maddison et al. (2014).
B COMPUTATION OF p(Sk), pD\C(S \ C) AND R(Sk, s)
We can sample the set Sk from the Plackett-Luce distribution using the Gumbel-Top-k trick by drawing Gumbel variables Gφi ∼ Gumbel(φi) for each element and returning the indices of the k largest Gumbels. If we ignore the ordering, this means we will obtain the set Sk if mini∈Sk Gφi > maxi∈D\Sk Gφi . Omitting the superscript k for clarity, we can use the Gumbel-Max trick, i.e. that GφD\S = maxi 6∈S Gφi ∼ Gumbel(φD\S) (equation 30) and marginalize over GφD\S :
p(S) = P (min i∈S Gφi > GφD\S )
= P (Gφi > GφD\S , i ∈ S)
= ∫ ∞ gφD\S=−∞ fφD\S (gφD\S )P (Gφi > gφD\S , i ∈ S)∂gφD\S
= ∫ ∞ gφD\S=−∞ fφD\S (gφD\S ) ∏ i∈S ( 1− Fφi(gφD\S ) ) ∂gφD\S (31)
= ∫ 1 u=0 ∏ i∈S ( 1− Fφi ( F−1φD\S (u) )) ∂u (32)
Here we have used a change of variables u = FφD\S (gφD\S ). This expression can be efficiently numerically integrated (although another change of variables may be required for numerical stability depending on the values of φ).
Exact computation in O(2k). The integral in equation 31 can be computed exactly using the identity ∏
i∈S (ai − bi) = ∑ C⊆S (−1)|C| ∏ i∈C bi ∏ i∈S\C ai
which gives
p(S) = ∫ ∞ gφD\S=−∞ fφD\S (gφD\S ) ∏ i∈S ( 1− Fφi(gφD\S ) ) ∂gφD\S
= ∑ C⊆S (−1)|C| ∫ ∞ gφD\S=−∞ fφD\S (gφD\S ) ∏ i∈C Fφi(gφD\S ) ∏ i∈S\C 1∂gφD\S
= ∑ C⊆S (−1)|C| ∫ ∞ gφD\S=−∞ eφD\S (gφD\S )FφD\S (gφD\S )FφC (gφD\S )∂gφD\S
= ∑ C⊆S (−1)|C| ∫ ∞ gφD\S=−∞ eφD\S (gφD\S )Fφ(D\S)∪C (gφD\S )∂gφD\S
= ∑ C⊆S (−1)|C|(1− 0) exp(φD\S) exp(φ(D\S)∪C)
= ∑ C⊆S (−1)|C| 1− ∑ i∈S p(i) 1− ∑ i∈S\C p(i) . (33)
Computation of pD\C(S \ C). When using the Gumbel-Top-k trick over the restricted domain D \ C, we do not need to renormalize the log-probabilities φs, s ∈ D \ C since the Gumbel-Top-k trick applies to unnormalized log-probabilities. Also, assuming C ⊆ Sk, it holds that (D \C)\ (S \ C) = D \ S. This means that we can compute pD\C(S \ C) similar to equation 31:
pD\C(S \ C) = P ( min i∈S\C Gφi > Gφ(D\C)\(S\C))
= P ( min i∈S\C Gφi > GφD\S )
= ∫ ∞ gφD\S=−∞ fφD\S (gφD\S ) ∏ i∈S\C ( 1− Fφi(gφD\S ) ) ∂gφD\S . (34)
Computation of R(Sk, s). Note that, using equation 10, it holds that∑ s∈Sk pD\{s}(Sk \ {s})p(s) p(Sk) = ∑ s∈Sk P (b1 = s|Sk) = 1
from which it follows that p(Sk) = ∑ s∈Sk pD\{s}(Sk \ {s})p(s)
such that
R(Sk, s) = pD\{s}(Sk \ {s})
p(Sk) = pD\{s}(Sk \ {s})∑ s′∈Sk p D\{s′}(Sk \ {s′})p(s′) . (35)
This means that, to compute the leave-one-out ratio for all s ∈ Sk, we only need to compute pD\{s}(Sk \{s}) for s ∈ Sk. When using the numerical integration or summation inO(2k), we can reuse computation, whereas using the naive method, the cost is O(k · (k−1)!) = O(k!), making the total computational cost comparable to computing just p(Sk), and the same holds when computing the ‘second-order’ leave one out ratios for the built-in baseline (equation 17).
Details of numerical integration. For computation of the leave-one-out ratio (equation 35) for large k we can use the numerical integration, where we need to compute equation 34 with C = {s}. For this purpose, we rewrite the integral as
pD\C(S \ C) = ∫ ∞ gφD\S=−∞ fφD\S (gφD\S ) ∏ i∈S\C ( 1− Fφi(gφD\S ) ) ∂gφD\S
= ∫ 1 u=0 ∏ i∈S\C ( 1− Fφi ( F−1φD\S (u) )) ∂u
= ∫ 1 u=0 ∏ i∈S\C ( 1− uexp(φi−φD\S) ) ∂u
= exp(b) · ∫ 1 v=0 vexp(b)−1 ∏
i∈S\C
( 1− vexp(φi−φD\S+b) ) ∂v
= exp(a+ φD\S) · ∫ 1 v=0 vexp(a+φD\S)−1 ∏
i∈S\C
( 1− vexp(φi+a) ) ∂v.
Here we have used change of variables v = uexp(−b) and a = b − φD\S . This form allows to compute the integrands efficiently, as
∏ i∈S\C ( 1− vexp(φi+a) ) = ∏ i∈S ( 1− vexp(φi+a) )∏ i∈C ( 1− vexp(φi+a)
) where the numerator only needs to computed once, and, since C = {s} when computing equation 35, the denominator only consists of a single term.
The choice of a may depend on the setting, but we found that a = 5 is a good default option which leads to an integral that is generally smooth and can be accurately approximated using the trapezoid rule. We compute the integrands in logarithmic space and sum the terms using the stable LOGSUMEXP trick. In our code we provide an implementation which also computes all second-order leave-one-out ratios efficiently.

heading: C THE SUM-AND-SAMPLE ESTIMATOR

heading: C.1 UNBIASEDNESS OF THE SUM-AND-SAMPLE ESTIMATOR
text: We show that the sum-and-sample estimator is unbiased for any set C ⊂ D (see also Liang et al. (2018); Liu et al. (2019)):
Ex∼pD\C(x) [∑ c∈C p(c)f(c) + ( 1− ∑ x∈C p(c) ) f(x) ]
= ∑ c∈C p(c)f(c) +
( 1−
∑ c∈C p(c) ) Ex∼pD\C(x)[f(x)]
= ∑ c∈C p(c)f(c) +
( 1−
∑ c∈C p(c) ) ∑ x∈D\C
p(x) 1− ∑ c∈C p(c) f(x)
= ∑ c∈C p(c)f(c) + ∑ x∈D\C p(x)f(x)
= ∑ x∈D p(x)f(x)
= Ex∼p(x)[f(x)]

heading: C.2 RAO-BLACKWELLIZATION OF THE STOCHASTIC SUM-AND-SAMPLE ESTIMATOR
text: In this section we give the proof that Rao-Blackwellizing the stochastic sum-and-sample estimator results in the unordered set estimator.
Theorem 4. Rao-Blackwellizing the stochastic sum-and-sample estimator results in the unordered set estimator, i.e.
EBk∼p(Bk|Sk) k−1∑ j=1 p(bj)f(bj) + 1− k−1∑ j=1 p(bj)  f(bk)  = ∑ s∈Sk p(s)R(Sk, s)f(s). (36)
Proof. To give the proof, we first prove three Lemmas.
Lemma 1.
P (bk = s|Sk) = p(Sk \ {s}) p(Sk)
p(s) 1− ∑ s′∈Sk\{s} p(s ′) (37)
Proof. Similar to the derivation of P (b1 = s|Sk) (equation 10 in the main paper), we can write:
P (bk = s|Sk) = P (Sk ∩ bk = s)
p(Sk)
= p(Sk \ {s})pD\(Sk\{s})(s)
p(Sk)
= p(Sk \ {s}) p(Sk)
p(s) 1− ∑ s′∈Sk\{s} p(s ′) .
The step from the first to the second row comes from analyzing the event Sk∩bk = s using sequential sampling: to sample Sk (including s) with s being the k-th element means that we should first sample Sk \ {s} (in any order), and then sample s from the distribution restricted to D \ (Sk \ {s}).
Lemma 2.
p(S) + p(S \ {s}) 1−
∑ s′∈S p(s ′)
1− ∑ s′∈S\{s} p(s ′) = pD\{s}(S \ {s}) (38)
Dividing equation 33 by 1− ∑ s′∈S p(s ′) on both sides, we obtain
Proof.
p(S) 1− ∑ s′∈S p(s ′)
= ∑ C⊆S (−1)|C| 1 1− ∑ s′∈S\C p(s ′)
= ∑
C⊆S\{s}
( (−1)|C| 1
1− ∑ s′∈S\C p(s ′) + (−1)|C∪{s}| 1 1− ∑ s′∈S\(C∪{s}) p(s ′)
)
= ∑
C⊆S\{s}
(−1)|C| 1 1− ∑ s′∈S\C p(s ′) + ∑ C⊆S\{s} (−1)|C∪{s}| 1 1− ∑ s′∈S\(C∪{s}) p(s ′)
= ∑
C⊆S\{s}
(−1)|C| 1 1− p(s)− ∑ s′∈(S\{s})\C p(s ′) − ∑ C⊆S\{s} (−1)|C| 1 1− ∑ s′∈(S\{s})\C p(s ′)
= 1 1− p(s) ∑
C⊆S\{s}
(−1)|C| 1 1− ∑ s′∈(S\{s})\C p(s′) 1−p(s) − p(S \ {s}) 1− ∑ s′∈S\{s} p(s ′)
= 1 1− p(s) pD\{s}(S \ {s}) 1− ∑ s′∈S\{s} p(s′) 1−p(s) − p(S \ {s}) 1− ∑ s′∈S\{s} p(s ′)
= pD\{s}(S \ {s}) 1− p(s)− ∑ s′∈S\{s} p(s ′) − p(S \ {s}) 1− ∑ s′∈S\{s} p(s ′)
= pD\{s}(S \ {s}) 1− ∑ s′∈S p(s ′) − p(S \ {s}) 1− ∑ s′∈S\{s} p(s ′) .
Multiplying by 1− ∑ s′∈S p(s ′) and rearranging terms proves Lemma 2.
Lemma 3.
p(s) + 1− ∑ s′∈Sk p(s′) P (bk = s|Sk) = p(s)R(Sk, s) (39)
Proof. First using Lemma 1 and then Lemma 2 we find
p(s) + 1− ∑ s′∈Sk p(s′) P (bk = s|Sk) =p(s) +
1− ∑ s′∈Sk p(s′)  p(Sk \ {s}) p(Sk)
p(s) 1− ∑ s′∈Sk\{s} p(s ′)
= p(s)
p(Sk)
( p(Sk) + 1− ∑ s′∈Sk p(s ′)
1− ∑ s′∈Sk\{s} p(s ′) p(Sk \ {s})
)
= p(s)
p(Sk) pD\{s}(Sk \ {s})
=p(s)R(Sk, s).
Now we can complete the proof of Theorem 4 by adding p(bk)f(bk) − p(bk)f(bk) = 0 to the estimator, moving the terms independent of Bk outside the expectation and using Lemma 3:
EBk∼p(Bk|Sk) k−1∑ j=1 p(bj)f(bj) + 1− k−1∑ j=1 p(bj)  f(bk) 
=EBk∼p(Bk|Sk)  k∑ j=1 p(bj)f(bj) + 1− k∑ j=1 p(bj)  f(bk)  = ∑ s∈Sk p(s)f(s) + EBk∼p(Bk|Sk) 1− ∑ s′∈Sk p(s′)  f(bk)  = ∑ s∈Sk p(s)f(s) + ∑ s∈Sk 1− ∑ s′∈Sk p(s′)
P (bk = s|Sk)f(s) = ∑ s∈Sk p(s) + 1− ∑ s′∈Sk p(s′) P (bk = s|Sk)  f(s)
= ∑ s∈Sk p(s)R(Sk, s)f(s).

heading: C.3 THE STOCHASTIC SUM-AND-SAMPLE ESTIMATOR WITH MULTIPLE SAMPLES
text: As was discussed in Liu et al. (2019), one can trade off the number of summed terms and number of sampled terms to maximize the achieved variance reduction. As a generalization of Theorem 4 (the stochastic sum-and-sample estimator with k − 1 summed terms), we introduce here the stochastic sum-and-sample estimator that sums k −m terms and samples m > 1 terms without replacement. To estimate the sampled term, we use the unordered set estimator on them samples without replacement, on the domain restricted to D \ Bk−m. In general, we denote the unordered set estimator restricted to the domain D \ C by
eUS,D\C(Sk) = ∑
s∈Sk\C
p(s)RD\C(Sk, s)f(s) (40)
where RD\C(Sk, s) is the leave-one-out ratio restricted to the domain D \ C, similar to the second order leave-one-out ratio in equation 18:
RD\C(Sk, s) = p (D\C)\{s} θ ((S k \ C) \ {s}) p D\C θ (S k \ C) . (41)
While we can also constrain Sk ⊆ (D \C), this definition is consistent with equation 18 and allows simplified notation. Theorem 5. Rao-Blackwellizing the stochastic sum-and-sample estimator with m > 1 samples results in the unordered set estimator, i.e.
EBk∼p(Bk|Sk) k−m∑ j=1 p(bj)f(bj) + 1− k−m∑ j=1 p(bj)  eUS,D\Bk−m(Sk)  = ∑ s∈Sk p(s)R(Sk, s)f(s).
(42)
Proof. Recall that for the unordered set estimator, it holds that eUS(Sk) = Eb1∼p(b1|Sk) [f(b1)] = Ex∼p(x) [ f(x) ∣∣x ∈ Sk] (43)
which for the restricted equivalent (with restricted distribution pD\C) translates into eUS,D\C(Sk) = Ex∼pD\C(x) [ f(x) ∣∣x ∈ Sk] = Ex∼p(x) [f(x)∣∣x ∈ Sk, x 6∈ C] . (44) Now we consider the distribution bk−m+1|Sk, Bk−m: the distribution of the first element sampled (without replacement) after sampling Bk−m, given (conditionally on the event) that the set of k samples is Sk, so we have bk−m+1 ∈ Sk and bk−m+1 6∈ Bk−m. This means that its conditional expectation of f(bk−m+1) is the restricted unordered set estimator for C = Bk−m since
eUS,D\B k−m (Sk) = Ex∼p(x) [ f(x) ∣∣x ∈ Sk, x 6∈ Bk−m] = Ebk−m+1∼p(bk−m+1|Sk,Bk−m) [f(bk−m+1)] . (45)
Observing that the definition (equation 42) of the stochastic sum-and-sample estimator does not depend on the actual order of the m samples, and using equation 45, we can reduce the multisample estimator to the stochastic sum-and-sample estimator with k′ = k − m + 1, such that the result follows from equation 36.
EBk∼p(Bk|Sk) k−m∑ j=1 p(bj)f(bj) + 1− k−m∑ j=1 p(bj)  eUS,D\Bk−m(Sk) 
=EBk−m∼p(Bk−m|Sk) k−m∑ j=1 p(bj)f(bj) + 1− k−m∑ j=1 p(bj)  eUS,D\Bk−m(Sk) 
=EBk−m∼p(Bk−m|Sk) k−m∑ j=1 p(bj)f(bj) + 1− k−m∑ j=1 p(bj) Ebk−m+1∼p(bk−m+1|Sk,Bk−m) [f(bk−m+1)] 
=EBk−m+1∼p(Bk−m+1|Sk) k−m∑ j=1 p(bj)f(bj) + 1− k−m∑ j=1 p(bj)  f(bk−m+1) 
=ESk−m+1|Sk EBk−m+1∼p(Bk−m+1|Sk−m+1) k−m∑ j=1 p(bj)f(bj) + 1− k−m∑ j=1 p(bj)  f(bk−m+1) 
=ESk−m+1|Sk ∑ s∈Sk p(s)R(Sk, s)f(s)  = ∑ s∈Sk p(s)R(Sk, s)f(s). (46)

heading: D THE IMPORTANCE-WEIGHTED ESTIMATOR

heading: D.1 RAO-BLACKWELLIZATION OF THE IMPORTANCE-WEIGHTED ESTIMATOR
text: In this section we give the proof that Rao-Blackwellizing the importance-weighted estimator results in the unordered set estimator. Theorem 6. Rao-Blackwellizing the importance-weighted estimator results in the unordered set estimator, i.e.:
Eκ∼p(κ|Sk) ∑ s∈Sk
p(s)
1− Fφs(κ) f(s)  = ∑ s∈Sk p(s)R(Sk, s)f(s). (47)
Here we have slightly rewritten the definition of the importance-weighted estimator, using that q(s, a) = P (gφs > a) = 1 − Fφs(a), where Fφs is the CDF of the Gumbel distribution (see Appendix A).
Proof. We first prove the following Lemma:
Lemma 4. Eκ∼p(κ|Sk)
[ 1
1− Fφs(κ)
] = R(Sk, s) (48)
Proof. Conditioning on Sk, we know that the elements in Sk have the k largest perturbed logprobabilities, so κ, the (k + 1)-th largest perturbed log-probability is the largest perturbed logprobability in D \Sk, and satisfies κ = maxs∈D\Sk gφs = gφD\Sk ∼ Gumbel(φD\Sk). Computing p(κ|Sk) using Bayes’ Theorem, we have
p(κ|Sk) = p(S k|κ)p(κ) p(Sk) =
∏ s∈Sk(1− Fφs(κ))fφD\Sk (κ)
p(Sk) (49)
which allows us to compute (using equation 34 with C = {s} and gφD\S = κ) Eκ∼p(κ|Sk) [ 1
1− Fφs(κ) ] =
∫ ∞ κ=−∞ p(κ|Sk) 1 1− Fφs(κ) ∂κ
= ∫ ∞ κ=−∞ ∏ s∈Sk(1− Fφs(κ))fφD\Sk (κ) p(Sk)
1
1− Fφs(κ) ∂κ
= 1
p(Sk) ∫ ∞ κ=−∞ ∏ s∈Sk\{s} (1− Fφs(κ))fφD\Sk (κ)∂κ
= 1
p(Sk) pD\{s}(S \ {s})
=R(Sk, s).
Using Lemma 4 we find
Eκ∼p(κ|Sk) ∑ s∈Sk
p(s)
1− Fφs(κ) f(s)  = ∑ s∈Sk p(s)Eκ∼p(κ|Sk) [ 1 1− Fφs(κ) ] f(s)
= ∑ s∈Sk p(s)R(Sk, s)f(s).

heading: D.2 THE IMPORTANCE-WEIGHTED POLICY GRADIENT ESTIMATOR WITH BUILT-IN BASELINE
text: For self-containment we include this section, which is adapted from our unpublished workshop paper (Kool et al., 2019b). The importance-weighted policy gradient estimator combines REINFORCE (Williams, 1992) with the importance-weighted estimator (Duffield et al., 2007; Vieira, 2017) in equation 15 which results in an unbiased estimator of the policy gradient∇θEpθ(x)[fθ(x)]:
eIWPG(Sk, κ) = ∑ s∈Sk pθ(s) qθ,κ(s) ∇θ log pθ(s)f(s) = ∑ s∈Sk ∇θpθ(s) qθ,κ(s) f(s) (50)
Recall that κ is the (k + 1)-th largest perturbed log-probability (see Section 3.2). We compute a lower variance but biased variant by normalizing the importance weights using the normalization W (Sk) = ∑ s∈Sk pθ(s) qθ,κ(s) .
As we show in Kool et al. (2019b), we can include a ‘baseline’ B(Sk) = ∑ s∈Sk pθ(s) qθ,κ(s)
f(s) and correct for the bias (since it depends on the complete sample Sk) by weighting individual terms of the estimator by 1− pθ(s) + pθ(s)qθ,κ(s) :
eIWPGBL(Sk, κ) = ∑ s∈Sk ∇θpθ(s) qθ,κ(s) ( f(s) ( 1− pθ(s) + pθ(s) qθ,κ(s) ) −B(Sk) ) (51)
For the normalized version, we use the normalizationW (Sk) = ∑ s∈Sk pθ(s) qθ,κ(s) for the baseline, and Wi(S k) = W (Sk)− pθ(s)qθ,κ(s) + pθ(s) to normalize the individual terms:
∇θEy∼pθ(y) [f(y)] ≈ ∑ s∈Sk
1 Wi(Sk) · ∇θpθ(s) qθ,κ(s)
( f(s)− B(S k)
W (Sk)
) (52)
It seems odd to normalize the terms in the outer sum by 1 Wi(Sk) instead of 1 W (Sk) , but equation 52 can be rewritten into a form similar to equation 17, i.e. with a different baseline for each sample, but this form is more convenient for implementation (Kool et al., 2019b).

heading: E THE UNORDERED SET POLICY GRADIENT ESTIMATOR

heading: E.1 PROOF OF UNBIASEDNESS OF THE UNORDERED SET POLICY GRADIENT ESTIMATOR WITH BASELINE
text: To prove the unbiasedness of result we need to prove that the control variate has expectation 0: Lemma 5.
ESk∼pθ(Sk) ∑ s∈Sk ∇θpθ(s)R(Sk, s) ∑ s′∈Sk pθ(s ′)RD\{s}(Sk, s′)f(s′)  = 0. (53) Proof. Similar to equation 10, we apply Bayes’ Theorem conditionally on b1 = s to derive for s′ 6= s
P (b2 = s ′|Sk, b1 = s) = P (Sk|b2 = s′, b1 = s)P (b2 = s′|b1 = s′) P (Sk|b1 = s)
= p D\{s,s′} θ (S k \ {s, s′})pD\{s}θ (s′) p D\{s} θ (S k \ {s})
= pθ(s
′)
1− pθ(s) RD\{s}(Sk, s′). (54)
For s′ = s we have RD\{s}(Sk, s′) = 1 by definition, so using equation 54 we can show that∑ s′∈Sk pθ(s ′)RD\{s}(Sk, s′)f(s′)
= pθ(s)f(s) + ∑
s′∈Sk\{s}
pθ(s ′)RD\{s}(Sk, s′)f(s′)
= pθ(s)f(s) + (1− pθ(s)) ∑
s′∈Sk\{s}
pθ(s ′)
1− pθ(s) RD\{s}(Sk, s′)f(s′)
= pθ(s)f(s) + (1− pθ(s)) ∑
s′∈Sk\{s}
P (b2 = s ′|Sk, b1 = s)f(s′)
= pθ(s)f(s) + (1− pθ(s))Eb2∼pθ(b2|Sk,b1=s) [f(b2)] = Eb2∼pθ(b2|Sk,b1=s) [pθ(b1)f(b1) + (1− pθ(b1))f(b2)] .
Now we can show that the control variate is actually the result of Rao-Blackwellization:
ESk∼pθ(Sk) ∑ s∈Sk ∇θpθ(s)R(Sk, s) ∑ s′∈Sk pθ(s ′)RD\{s}(Sk, s′)f(s′)  = ESk∼pθ(Sk) ∑ s∈Sk pθ(s)R(S k, s)∇θ log pθ(s) ∑ s′∈Sk pθ(s ′)RD\{s}(Sk, s′)f(s′)
 = ESk∼pθ(Sk) ∑ s∈Sk P (b1 = s|Sk)∇θ log pθ(s) ∑ s′∈Sk pθ(s ′)RD\{s}(Sk, s′)f(s′)
 = ESk∼pθ(Sk) Eb1∼pθ(b1|Sk) ∇θ log pθ(b1) ∑
s′∈Sk pθ(s
′)RD\{b1}(Sk, s′)f(s′)  = ESk∼pθ(Sk) [ Eb1∼pθ(b1|Sk) [ ∇θ log pθ(b1)Eb2∼pθ(b2|Sk,b1) [pθ(b1)f(b1) + (1− pθ(b1))f(b2)]
]] = ESk∼pθ(Sk) [ EBk∼pθ(Bk|Sk) [∇θ log pθ(b1) (pθ(b1)f(b1) + (1− pθ(b1))f(b2))]
] = EBk∼pθ(Bk) [∇θ log pθ(b1) (pθ(b1)f(b1) + (1− pθ(b1))f(b2))]
This expression depends only on b1 and b2 and we recognize the stochastic sum-and-sample estimator for k = 2 used as ‘baseline’. As a special case of equation 13 for C = {b1}, we have
Eb2∼pθ(b2|b1) [(pθ(b1)f(b1) + (1− pθ(b1))f(b2))] = Ei∼pθ(i) [f(i)] . (55)
Using this, and the fact that Eb1∼pθ(b1) [∇θ log pθ(b1)] = ∇θEb1∼pθ(b1) [1] = ∇θ1 = 0 we find
ESk∼pθ(Sk) ∑ s∈Sk ∇θpθ(s)R(Sk, s) ∑ s′∈Sk pθ(s ′)RD\{s}(Sk, s′)f(s′)  = EBk∼pθ(Bk) [∇θ log pθ(b1) (pθ(b1)f(b1) + (1− pθ(b1))f(b2))] = Eb1∼pθ(b1) [ ∇θ log pθ(b1)Eb2∼pθ(b2|b1) [(pθ(b1)f(b1) + (1− pθ(b1))f(b2))]
] = Eb1∼pθ(b1) [ ∇θ log pθ(b1)Ex∼pθ(x) [f(x)]
] = Eb1∼pθ(b1) [∇θ log pθ(b1)]Ex∼pθ(x) [f(x)] = 0 · Ex∼pθ(x) [f(x)] = 0

heading: F THE RISK ESTIMATOR

heading: F.1 PROOF OF BUILT-IN BASELINE
text: We show that the RISK estimator, taking gradients through the normalization factor actually has a built-in baseline. We first use the log-derivative trick to rewrite the gradient of the ratio as the ratio times the logarithm of the gradient, and then swap the summation variables in the double sum that arises:
eRISK(S) = ∑ s∈S ∇θ ( pθ(s)∑ s′∈S pθ(s ′) ) f(s)
= ∑ s∈S pθ(s)∑ s′∈S pθ(s ′) ∇θ log ( pθ(s)∑ s′∈S pθ(s ′) ) f(s)
= ∑ s∈S pθ(s)∑ s′∈S pθ(s ′)
( ∇θ log pθ(s)−∇θ log
∑ s′∈S pθ(s ′)
) f(s)
= ∑ s∈S pθ(s)∑ s′∈S pθ(s ′) ( ∇θpθ(s) pθ(s) − ∑ s′∈S ∇θpθ(s′)∑ s′∈S pθ(s ′) ) f(s)
= ∑ s∈S ∇θpθ(s)f(s)∑ s′∈S pθ(s ′) − ∑ s,s′∈S pθ(s)∇θpθ(s′)f(s)(∑ s′∈S pθ(s ′) )2
= ∑ s∈S ∇θpθ(s)f(s)∑ s′∈S pθ(s ′) − ∑ s,s′∈S pθ(s ′)∇θpθ(s)f(s′)(∑ s′∈S pθ(s ′) )2
= ∑ s∈S ∇θpθ(s)∑ s′∈S pθ(s ′) ( f(s)− ∑ s′∈S pθ(s ′)f(s′)∑ s′∈S pθ(s ′) )
= ∑ s∈S ∇θpθ(s)∑ s′′∈S pθ(s ′′)
( f(s)−
∑ s′∈S pθ(s ′)∑ s′′∈S pθ(s ′′) f(s′)
) .

heading: G CATEGORICAL VARIATIONAL AUTO-ENCODER

heading: G.1 EXPERIMENTAL DETAILS
text: We use the code6 by Yin et al. (2019) to reproduce their categorical VAE experiment, of which we include details here for self-containment. The dataset is MNIST, statically binarized by thresholding at 0.5 (although we include results using the standard binarized dataset by Salakhutdinov & Murray (2008); Larochelle & Murray (2011) in Section G.2). The latent representation z is K = 20 dimensional with C = 10 categories per dimension with a uniform prior p(zk = c) = 1/C, k = 1, ...,K. The encoder is parameterized by φ as qφ(z|x) = ∏ k qφ(zk|x) and has two fully connected hidden layers with 512 and 256 hidden nodes respectively, with LeakyReLU (α = 0.1) activations. The decoder, parameterized by θ, is given by pθ(x|z) = ∏ i pθ(xi|z), where xi ∈ {0, 1} are the pixel values, and has fully connected hidden layers with 256 and 512 nodes and LeakyReLU activation.
ELBO optimization. The evidence lower bound (ELBO) that we optimize is given by L(φ,θ) = Ez∼qφ(z|x) [ln pθ(x|z) + ln p(z)− ln qφ(z|x)] (56)
= Ez∼qφ(z|x) [ln pθ(x|z)]−KL(qφ (z|x)||p(z)) . (57)
For the decoder parameters θ, since qφ(z|x) does not depend on θ, it follows that ∇θL(φ,θ) = Ez∼qφ(z|x) [∇θ ln pθ(x|z)] . (58)
For the encoder parameters φ, we can write∇φL(φ,θ) using equation 57 and equation 19 as ∇φL(φ,θ) = Ez∼qφ(z|x) [∇φ ln qφ(z|x) ln pθ(x|z)]−∇φKL(qφ (z|x)||p(z)) . (59)
6https://github.com/ARM-gradient/ARSM
This assumes we can compute the KL divergence analytically. Alternatively, we can use a sample estimate for the KL divergence, and use equation 56 with equation 19 to obtain
∇φL(φ,θ) = Ez∼qφ(z|x) [∇φ ln qφ(z|x)(ln pθ(x|z) + ln p(z)− ln qφ(z|x)) +∇φ ln qφ(z|x)] (60)
= Ez∼qφ(z|x) [∇φ ln qφ(z|x)(ln pθ(x|z)− ln qφ(z|x))] . (61)
Here we have left out the term Ez∼qφ(z|x) [∇φ ln qφ(z|x)] = 0, similar to Roeder et al. (2017), and, assuming a uniform (i.e. constant) prior ln p(z), the term Ez∼qφ(z|x) [∇φ ln qφ(z|x) ln p(z)] = 0. With a built-in baseline, this second term cancels out automatically, even if it is implemented. Despite the similarity of the equation 56 and equation 57, their gradient estimates (equation 60 and equation 59) are structurally dissimilar and care should be taken to implement the REINFORCE estimator (or related estimators such as ARSM and the unordered set estimator) correctly using automatic differentiation software. Using Gumbel-Softmax and RELAX, we take gradients ‘directly’ through the objective in equation 57.
We optimize the ELBO using the analytic KL for 1000 epochs using the Adam (Kingma & Ba, 2015) optimizer. We use a learning rate of 10−3 for all estimators except Gumbel-Softmax and RELAX, which use a learning rate of 10−4 as we found they diverged with a higher learning rate. For ARSM, as an exception we use the sample KL, and a learning rate of 3 · 10−4, as suggested by the authors. All reported ELBO values are computed using the analytic KL. Our code is publicly available7.

heading: G.2 ADDITIONAL RESULTS
text: Gradient variance during training. We also evaluate gradient variance of different estimators during different stages of training. We measure the variance of different estimators with k = 4 samples during training with REINFORCE with replacement, such that all estimators are computed for the same model parameters. The results during training, given in Figure 4, are similar to the results for the trained model in Table 1, except for at the beginning of training, although the rankings of different estimator are mostly the same.
Negative ELBO on validation set. Figure 5 shows the -ELBO evaluated during training on the validation set. For the large latent space, we see validation error quickly increase (after reaching a minimum) which is likely because of overfitting (due to improved optimization), a phenomenon observed before (Tucker et al., 2017; Grathwohl et al., 2018). Note that before the overfitting starts, both REINFORCE without replacement and the unordered set estimator achieve a validation error similar to the other estimators, such that in a practical setting, one can use early stopping.
Results using standard binarized MNIST dataset. Instead of using the MNIST dataset binarized by thresholding values at 0.5 (as in the code and paper by Yin et al. (2019)) we also experiment with the standard (fixed) binarized dataset by Salakhutdinov & Murray (2008); Larochelle & Murray (2011), for which we plot train and validation curves for two runs on the small and large domain in Figure 6. This gives more realistic (higher) -ELBO scores, although we still observe the effect of overfitting. As this is a bit more unstable setting, one of the runs using REINFORCE with replacement diverged, but in general the relative performance of estimators is similar to using the dataset with 0.5 threshold.
7https://github.com/wouterkool/estimating-gradients-without-replacement

heading: H TRAVELLING SALESMAN PROBLEM
text: The Travelling Salesman Problem (TSP) is a discrete optimization problem that consists of finding the order in which to visit a set of locations, given as x, y coordinates, to minimize the total length of the tour, starting and ending at the same location. As a tour can be considered a sequence of locations, this problem can be set up as a sequence modelling problem, that can be either addressed using supervised (Vinyals et al., 2015) or reinforcement learning (Bello et al., 2016; Kool et al., 2019a).
Kool et al. (2019a) introduced the Attention Model, which is an encoder-decoder model which considers a TSP instances as a fully connected graph. The encoder computes embeddings for all nodes (locations) and the decoder produces a tour, which is sequence of nodes, selecting one note at the time using an attention mechanism, and uses this autoregressively as input to select the next node. In Kool et al. (2019a), this model is trained using REINFORCE, with a greedy rollout used as baseline to reduce variance.
We use the code by Kool et al. (2019a) to train the exact same Attention Model (for details we refer to Kool et al. (2019a)), and minimize the expected length of a tour predicted by the model, using different gradient estimators. We did not do any hyperparameter optimization and used the exact same training details, using the Adam optimizer (Kingma & Ba, 2015) with a learning rate of 10−4 (no decay) for 100 epochs for all estimators. For the baselines, we used the same batch size of 512, but for estimators that use k = 4 samples, we used a batch size of 5124 = 128 to compensate for the additional samples (this makes multi-sample methods actually faster since the encoder still needs to be evaluated only once).

Review summaries:
In this paper, an unbiased estimator for expectations over discrete random variables is developed based on a sampling-without-replacement strategy. The proposed estimator is shown to be a Rao-Blackwellization of three existing unbiased estimators with guaranteed reduction in estimation variance. The connections of the method to other gradient estimators are discussed. Experimental results on several toy and real-data DL/RL problems are reported to demonstrate the applicability of the proposed estimators in the practice of machine learning. 

This paper introduces an gradient estimator for loss functions that are expectations over discrete random variables. The basic idea is that an estimator over a discrete distribution can be Rao-Blackwellized by conditioning on the event that the discrete realization was produced by being the first sample drawn from an unordered set of samples drawn with replacement. Much of the paper is spent showing how this Rao-Blackwellized estimator can be computed in practice and how it compares to other known estimators.