U-Statistics for Importance-Weighted Variational Inference
Authors: Javier Burroni
Kenta Takatsu
Justin Domke
Daniel Sheldon
Venue: Accepted by TMLR ()

Abstract: We propose the use of U-statistics to reduce variance for gradient estimation in importance-weighted variational inference.
The key observation is that, given a base gradient estimator that requires $m > 1$ samples and a total of $n > m$ samples to be used for estimation, lower variance is achieved by averaging the base estimator on overlapping batches of size $m$ than  disjoint batches, as currently done.
We use classical U-statistic theory to analyze the variance reduction, and propose novel approximations with theoretical guarantees to ensure computational efficiency.
We find empirically that U-statistic variance reduction can lead to modest to significant improvements in inference performance on a range of models, with little computational cost.


Published in Transactions on Machine Learning Research (02/2023)
U-Statistics for Importance-Weightee Variational Inference
Javier Burroni
jburvoni&cs. www.://iisssssssss.oo.om eduu
University of Massachusetts Amherst
Kenta Takatsu
kakalsu@@andreu.cmee edu
Carnegie Mellon University
Justin Domke University of Mossachusetts Amherst
domke wmass. edu
Daniel Sheldon
sheldon wmass. edu
University of Massachusetts Amherst
Reviewed on OpenReview: https: openreview. net/ forum? id= oXmuAP18 Vw
Abstract
We propose the of U-statistics to reduce variance for gradient estimation in importance- weighted variational inference. The key observation is that, given base gradient estimator that requires > samples and total of n > samples to be used for estimation, lower variance is achieved by averaging the base estimator on overlapping batches of size m than disjoint batches, as currently done. We use classical U-statistic theory to analyze the variance reduction, and propose novel approximations with theoretical guarantees to ensure computational efficiency. We find empirically that U-statistic variance reduction can lead to modest to significant improvements in inference performance on range of models, with little computational cost.
1 Introduction
An important recent development in variational inference (VT) the use of ideas from Monte Carlo sampling obtain tighter variational bounds (Burda et al., 2016; Maddison al., 2017; Le al., 2018; Naesseth al., 2018; Domke & Sheldon, 2019). Burda al. (2016) first introduced the importonce- weighted autoencoder (IWAE), deep generative model that uses the importance-weighted evidence lower bound (IW-ELBO) as its variational objective. The IW-ELBO uses samples from a proposal distribution to bound the lig-likelihood more tightly than the conventional evidence lower bound (ELBO), which uses only 1 sample. Later, the IW-ELBO was also connected to obtaining better approximate posterior distributions for pure inference applications of VI (Cremer et al., 2017; Domke & Sheldon, 2018), or "IWVI". Similar connections were made for other variational bounds (Naesseth al., 2018; Domke & Sheldon, 2019).
The IW-ELBO is attractive because, under certain assumptions [see Burda al. (2016); Domke & Sheldon (2018)), it gives tunable knob to make VI more accurate with more computation. The most obvious downside is the increased computational cost (up to factor of m) to form single estimate of the bound and its gradients. A more subtle tradeoff is that the signal-to- noise ratio of some gradient estimators degrades with m (Rainforth et al., 2018), which makes stochastic optimization ion of the bound harder and might hurt overall inference performance.
To take advantage of the tighter bound while controlling variance, one can average over independent replicates of base gradient estimator (Rainforth al., 2018). This idea is often used in practice and requires total of
samples from the proposal distribution.
Work done while UMass.

---

Published in Transactions on Machine Learning Research (02/2023)
Our main contribution is the observation that, whenever using r 1 replicates, it is possible to reduce variance with little computational overhead using ideas from the theory of U-statistics. Specifically, instead of rumning the base estimator on r independent batches of m samples from the proposal distribution and averaging the result, using the same n = rm samples we can run the estimator on à > r overlapping batches of m samples and average the result. In practice, the extra computation from using more batches is a small fraction of the time for model computations that are already required to be done for cach of the n samples. Specifically:
We describe how to take an m-sample base estimator for the IW-ELBO or its gradient and reduce variance compared to averaging over r replicates by forming a complete U-statistic, which averages the base estimator applied to every distinct batch of size m. This estimator has the lowest variance possible among estimators that average the base estimator over different batches, but it is usually not tractable in practice due the very large number of distinct batches.
We then show how to achieve most of the variance reduction with much less computation by using incomplete U-statistics, which average over a smaller number of overlapping batches. We introduce a novel way of selecting batches and prove that it attains a (1 1/E) fraction of the possible variance reduction with fr batches.
As an alternative to incomplete U-statistics, we introduce novel and fast approximations for IW-ELBO complete U-statisties. The extra computational step compared to the standard estimator is a single sort of the input samples, which is very fast. We prove accuracy bounds and show the approximations perform very well, especially in earlier iterations of stochastic optimization.
We remonstrate on a diverse set of inference problems that U-statistic-based variance reduction for the IW-ELBO either does not change, or leads to modest to significant gains in black-box VI performance, with no substantive downsides. We recommend always applying these techniques for black-box IWVI with r > 1.
We empirically show that U-statistic-basee estimators also reduce variance during IWAE training and lead to models with higher training objective values when used with either the standard gradient estimator the doubly-reparamettriae gradient (DReG) estimator (Tucker et al., 2018).
2 Importance-Weighee. Variational Inference
Assume target distribution p(2, æ) where E dx is observed and E R2 is latent. VI uses the following evidence lower bound (ELBO), given approximating distribution 4s with parameters ø E Rd,, to approximate p(x( (Saul et al., 1996; Blei et al., 2017):
p( 1 x)
L L= E
In
Z Gp-
(2(
The inequality follows from Jensen's inequality and the fact that E p(Z.z) go(2)
p()), that is, the importance
weight p(Z.z) is an unbiased estimate of p(x).
(o(Z)
Burda et al. (2016) first showed that a tighter bound can be obtained by using the average of m importance
weights within the logarithm. The importance-weighhed ELBO (IW-ELBO) is
Cm = E In
m
p(Z,, x²)
go (ii) < Inp(F),
lid Z₁ go.
(1)
This bound again follows from Jensen's incquality and the fact that 1 Em M(Z,.))
(*(2)) which is the sample
average of m unbiased estimates, remains unbiased for p(z). Moreover, we expect Jensen's inequality to provide a tighter bound bocause the distribution of this sample average more concentrated around p(x)
than the distribution of one estimate. Indeed, L. 2 Lm for m m' and C. + In p(x()) as m - oo (Burda et al., 2016).
2

---

Published in Transactions on Machine Learning Research (02/2023)
In inpirtanne--eiglttdd VI (IWVI), the IW-ELBO Lm is maximized with respect to the variational param- eters & to obtain the tightest possible lower bound to In p(x), which simultaneously finds an approximating distribution that is close in KL divergence to p(z ()) (Domke & Sheldon, 2018). In practice, the IW-ELBO and its gradients are estimated by sampling within a stochastic optimization routine. It convenient to define the log-weight random variables V2 = In 22, x) - In (Z())) for Z ~ 9s and rewrite the IW-ELBO as
h((1:m) In 1
fm R(h(V1mm))
(2)
Then, an unbiased IW-ELBO estimate with replicates, using n = rm i.idd. log-weights (V))) is
mn = Σ
(3)
En,m, we use the subscript n to denote the total number of input samples used for estimation and m for the number of arguments of h, which determines the IW-ELBO objective be optimized.
For gradient estimation, unbiased estimate for the TW-ELBO gradient VaLm is:
On. = - - 1( Zj,m), 5=1
ud Zja ~ Gd
(4)
where g(21m)) is any one of several unbiased "base" gradient estimators that operates on a batch of m samples from 9a, including the reparameterization gradient estimator (Kingma & Welling, 2013; Rezende et al., 2014), the doubay-reparrmeeeriied gradient (DReG) estimator (Tucker et al., 2018), or the score function estimator (Fu, 2006; Kleijnen & Rubinstein, 1996).
2.1 IWVI Tradeoffs: Bias, Variance, and Computation
Past research has shown that by using a tighter variational bound, IWVI can improve both learning and inference performance, but also introduce tradeoffs such as those pointed out by Rainforth et al. (2018). fact, there are several knobs to consider when using IWVI that control its bias, variance, and amount of computation. These tradcoffs our setting, with the goal of understanding when and how IWVI can be helpful and providing self-contained evidence that the setting where U-statistics are beneficial can and does arise in practice.
be complex so it is helpful to review the key elements as they relate to
Consider the task of maximizing an IW-ELBO objective L. to obtain the tightest final bound on the log- likelihood. This requires estimating L. and its gradient with respect to the variational parameters in each iteration of stochastic optimization procedure. Assume there a fixed budget of independent samples per iteration, where, for convenience, be adjusted to control the estimation bias and variance at the cost of increased computation. Specifically:
rm for an integer 2 1, above. The parameters m and r can
For a fixed m, by setting r' >> r,, we can reduce the variance of the estimator in Equation (4) by increasing the computational cost. to r'm rm samples per iteration.
For fixed r, by setting m' > m can reduce the bias of the objective - that is, the gap in the bound Lm < Inp(x)) - by increasing the computational cost to rm' > rm samples per iteration.
However, Rainforth et al. (2018) observed that increasing m may also have the negative effect of worsening the signal-to-noise (SNR) ratio of gradient estimation (but also that this could be counterbalande by increasing r). Later, Tucker et al. (2018) showed that, for the DReG gradient estimator, increasing can increase SNR; see also the paper by (Finke & Thiery, 2019) for a detailed discussion of these issues.
Overall, while the effoct of increasing the number of replicates r to roduce variance is quite clear, the effects of increasing m are sufficiently complex that is difficult to predict advance when it will be beneficial.
3

---

Published in Transactions on Machine Learning Research (02/2023)
1e-5
However, an important premise of our work is that the optimal setting of m is often strictly between and n, since this is the setting where U-statistics can be used reduce variance. To understand this, we can first reason from the perspective of user that is willing to spend more computation to get better model. Assuming the variational bound is not already tight, this user can increase m as much as desired to tighten the bound, and then, increase r as needed to control the gradient variance. This argument predicts that, for a sufficiently large com- putational budget and complex enough model (so that the bound is not already tight with 1), value > 1 will often be optimal.
3.75
3.50
3.25
3.00
2.75
16
2.50
128
256
2.00
8- 32 128
Figure 1: Distribution of the distance (error) between the distribution's covariance and that of the approx- imating distribution as a function of m for different numbers of sampled points n, after training an approx- imating distribution using the standard IW-ELBO es- timator. As we increase n, the optimal m also in- creases, but slow rate. [See Section 6 for details.]
From the perspective of a user with fixed computa- tional budget, in which the number of optimization iterations is also being fixed, we could instead ask: "for a fixed n, what are the optimal choices of m and r == n/m"? This question can be addressed em- pirically. Rainforth et al. (2018) reported in their Figure 6 that the extreme values, i.e, m 1 or m = n. were never the best values. We found empirically that for some models, this result also holds for black box VI, ie., the optimal choice of m is strictly greater than 1 and less than n, as shown in Figure 1. See also Figure 5, which shows that similar observations apply when using the DReG estimator. In our analysis of real Stan and UCI models, with n = 16, around half of them achieved the best performance for an intermediate value of m, depending on the approximating distribution and base gradient estimator seee Table 8 and 9 in Appendix G). And we further conjecture that the fraction of real-world models with this property will increase as increases.
For the rest of this work we focus on methods that can reduce variance for the case when 1 < m < n.
3 U-Statistic Estimators
We now introduce estimators for the IW-ELBO and its gradients based on U-statistics, and apply the theory of U-statistics to relate their variances. The theory of U-statistics was developed a seminal work by Hoeffding (1948) and extends the theory of unbiased estimation introduced by Halmos (1946). For detailed background, see the original works or the books by Lee (1990) and van der Vaart (2000).
The standard estimators in Eqs. (3) and (4) average the basc estimators h(vim) and g(31:)) on disjoint batches of the input samples. The key insight of U-statistics is that variance can be reduced by averaging the base estimators on a larger number of overlapping sets of samples.
We will consider general IW-ELBO estimators of the form
[s(V1.)) = Σ hhuee...... Dom.)
1
(5)
S
SES
where is any non-empty collection of size-m subsets of the indices [m] == fI,....,n,, .n....,] and 8; is the ith smallest index in the set s E S. Since E ((V.m))) = Cm, is clear (by symmetry and linearity) that T. (Vi.)) Lm, that is, the estimatoor unbiased. For now, we will call this a "U-statistic with kernel h", as it is clear the same construction can be generalized by replacing h by any other symmetric function of variables' or "kernel", while preserving the expected value. Later, will distinguish between different types of U-statistics based on the collection S.
Recall that symmetric function function invariant under all permutations of arguments.
4

---

Published in Transactions on Machine Learning Research (02/2023)
We can form U-statistics for gradient estimators by using base gradient estimators as kernels. Let g(31m)) be any symmetric base estimator such that ggg( (1mm)))) VaLm- The corresponding U-statistic is
(7
1
Σ g (( ..... 7mmm
(6)
|5]
SES
and satisfies (9 = VoLm.
3.1 Variance Comparison
How much variance reduction is possible for IWVI by using U-statistics? In this section, we first define the standard TW-ELBO estimator and complete U-statistic IW-ELBO estimator, and then relate their variances. For concreteness, restrict our attention to IW-ELBO objective estimators, but analagous results hold for gradients by using basc gradient cstimator as the kernel of the U-statistic.
We first express the standard IW-ELBO estimator 6n,m the terminology of Eq. (5):
Estimator The standard IW-ELBO estimator fn,m of Eq. (3) is the U-statistic Ls formed by taking S to be partition of [n]] into disjoint sets, i.e., S = ((I....,,), {m {m+1,,... 2m)..... {(r-((r--- 1)m +,...,, +l...,,m. rm)]
Estimator 2. The complete U-statistic IW-ELBO estimator C is the U-statistic Ls with S = (Inl) the set of all distinct subsets of [n] with exactly m elements.
We will show that the variance of the L is never more than that Ên.m: and is strictly less under certain conditions (that occur in practice), using classical bounds on U-statistic variance due to Hoeffding (1948).
Since ńU is an average of terms, one for each s E (Ixl), its variance depends on the covariances between pairs of terms for index sets s and s', which in turn depend on how many indices are shared by s and s'. This motivates the following definition:
Definition 3.1. Let Vi...., Vem be ii.d. log- weight,, For S c S m, take s, s, E ((2ml) with |s sss = c. Using h from Eq. (2), define
Se CovCov(((V,, Vam).. h(V Ver.
which depends only on and not the porticular and s'.
In words, this the covariance between two IW-ELBO estimates, each using one batch of i.i.d. log-whights, and where the two batches share log-weights in common. For example, when m have
So == 0, C1 = =ovo[lnn(e + 1),,,, In(ee Vh + 1)%,], and, C2 = Var [ar[in(3e + k))
Then, due to Hoeffding's classical result, Proposition 3.2. With S1, and Sm defined as above, the standard IW-ELBO estimator Pn.m (Estimator 1) and complete U-stalistic estimator (Estimator 2) with n = rm and r E salisfy
m² G1 VarlCh, " Sm Var Var[Lnnm)]
Moreover, for a fixed m. the quantity n Varl lends to ils lower bound m²C1 as increases.
Proof. The inequalities and asymptotic statement follow directly from Theorem 5.2 of Hoeffding (1948). The equality follows from the definition (m.
Hoeflding proved that mς1 S Cm- We observe in practice that there is gap bet the two variances that leads to practical gains for the complete U-statistic estimator in real VI problems.
A classical result of Halmos (1946) also shows that complete U-statistics are optimal in a certain sense: we describe how this result. applies estimator L0 f in Appendix B.
5

---

Published in Transactions on Machine Learning Research (02/2023)
Finally, we conchude this discussion by stating the main analogue of Proposition 3.2 for gradient estimation. The result, also following from Theorem 5.2 of Hoeffding (1948), states that the complete U-statistic gradient estimator has total variance and expected squared norm no larger than that of the standard estimator:
Proposition 3.3. Let 9n.m and 9 be the standard and complete- U-statistic gradient estimators formed using symmetric base gradient estimator g(21m)) that is unbiased for V VaLm and the same inder sels as fn.m and C respectively. Then tr( ar(Var(Sn,m)) S tr tr(Var((n,..) and E 11 11 S 1E/|s1155 181|
We provide proof in Appendix B.1.
3.2 Computational Complexity
There are two main factors to consider for the computational complexity of an IW-ELBO estimator:
1) The cost to compute leightights Vi = Inpp(ZZ, 2²) - In (Z,) for i E [n], and
2) the cost to compute the estimator given the bog-weights.
A problem with the complete U-statistics S and 5% is that they use ((())) (w) distinct subsets of indices in Step 2), which is expensive. It should be noted that these log- log-weight manipulations are very simple, while, for many probabilistic models, computing each log-weight is expensive, so, for modest m and n, the computation may still be dominated by Step 1). However, for large enough m and n, Step 2) is impractical.
4 Incomplete U-Statistic Estimators
In practice, we can achieve most of the variance reduction of the complete U-statistic with only modest computational cost. by averaging over only h < (*) subsets of indices selected in some way. Such an estimator is called an incomplete U-statistic. Incomplete U-statistics were introduced and studied by Blom (1976).
A general incomplete U-statistic for the IW-ELBO has the form in Eq. (5) where S (Im)) is collection of size-m subsets of [n] that does not include every possible subset. We will also allow S to be multi-set, that the same subset may appear more than once. Note that the standard IW-ELBO estimator Cn,m is itself an incomplete U-statistic, where the k- = index sets are disjoint. We can improve on this by selecting t: > r sets.
Estimator 3 (Random subsets). The random-subedd incomplete- Caamppete-U-statiinn U-statistic estimator for the IW-ELBO is the estimator Ls, where Sy is a set of le subsets (s:) drawn uniformly at random (with replacement) from (la)).
We next introduce a novel incomplete U-statistic, which is both very simple and enjoys strong theoretical propertios.
Estimator 4 (Permuted block). The permuted block estimator is computed by repeating the standard IW-ELBO estimator é times with randomly permuted bog-weights and averaging the results. Formally, the permuded-bee incomplete- U-statistie estimator for the IW-ELBO is the estimator Lsg with the collection st defined as follows. Let π denote a permutation of [n]. Define Sş as the collection obtained by permuting indices according to and then dividing them into disjoint sets of size m. That is,
S.
............... [(m)) (m)),, {m(m= + D)...,,((m)) ..... (=(((= 1)m + 1), ... (rm)))
+1)....,((mm)]}
Now, let St = Unen Sn where II is a collection of & random permutations and ( denotes union as a multiset. The total number of sets in St is k = ré.
Both incompleto- U-statisic estimators can achieve variance reduction in practice for a large enough number of sets but the permuted block estimator has an advantage: its variance with k subsets is never more than that of the random subset estimator with k subsets, and never more than the variance of the standard

---

Published in Transactions on Machine Learning Research (02/2023)
IW-ELBO estimator (and usually smaller). On the other hand, the variance of the random subset estimator is more than that of the standard estimator unless k 2 ko for some threshold ko > r.
Proposition 4.1. Given m and rm, the variances of these estimators satisfy the following partial ordering:
standard
(a) Varl permuted black
Var Var[
Var Var(Sn,m m,
VarlSsg] Lsg
(7)
Var [Save]
complete
random subset
Moreover, if the number of permutations f and Varl Var then (b) is strict; ifr= 1, then (c) is strict. (Note that the permuted and random subset estimators both use k = rł subsets.)
Proof. By Def. 3.1, if and s are uniformly drawn from (lll) and ((1r)) we have
B[n(V,,. ... Vs,, (V
Vvr
E[G = Σ
EICH, 12 Varl Cy
(8)
x2
R²
(,
(()))(
Let ........ . We be the random permutations. Observe that for s, E Sx. distinct, i.e., two distinct sets within the ith block, s and s' are disjoint and then h(Ve..... V...) ) is independent of h(Ve,..... Vein). Hence, all dependencies between different sets are due to relations between permutations, i.e., each of the fr terms will have a dependency with the (f (1-1)r terms not in the same permutation. Therefore, it follows from (8) that the total variance of Lsg is
VarlCsş] 5m + (1- 7) Var(Cy,m)
(9)
i.e., a convex combination of fm Var[Cn,m] and Var(Cy,ml Hence, using Proposition 3.2, (a) and (b) holds.
By a similar argument, the total variance of Ls.. is
VarlCs,d 1 Sm (1- 1e) Var(Cy
Then, (c) holds because
VarlCsg] Varl = t -1) Var [r S 0.
A remarkable property of the permuted-block estimator is that we can choose the number of permutations to guarantee what fraction of the variance reduction of the complete estimator we want to achieve. Say we would like to achieve 90% of the variance reduction; then it suffices to set & = 10. The following Proposition formalizes this result.
Proposition 4.2. Given m and n = rm, for & E N the permuted--laed block estimator achieves a (1 (1-1/E) fraction of the variance reduction provided by the complete U-slatistic TW-ELBO estimator, i.e.,
Var[Cn.m]
Var[Cs;]
= 7)( [Var[ Var LV (H........
standard
permuted block
standard
complete
Proof. This follows directly from Eq. (9).
The conclusions of Propositions 4.1 and 4.2 do not depend on the kernel. This means they provide strong guarantees for our novel and simple permuted-block incomplete U-statistic with any kernel, which may be of general interest, and also imply the following result for gradients:
Proposition 4.3. The conclusion of Proposition 4.2 holds with Var(E] replaced by either [|[ or tr([Var())), for each pair (L, 9) of objective estimator and gradient estimator that use the same collection S of index sets, and for any base gradient estimator g(vim))-

---

Published in Transactions on Machine Learning Research (02/2023)
5 Efficient Lower Bounds
In the last section, we approximated the complete U-statistic by averaging over < (m) subsets. For example, by Proposition 4.2, we could achieve 90% of the variance reduction with 10x more batches than the standard estimator, and the extra running-time cost is often very small in practice. An even faster alternative is to approximate the kernel in such way that iterating over subsets. In this section, we introduce such approximation for the IW-ELBO objective, where the extra running- time cost is a single sort of the n log neg-weights, which is extremely fast. Furthermore, Proposition 5.2 below will show that it is always a lower bound to L and has bounded approximation error, so expectation lower bounds f. and In p(z); thus, can be used as a surrogate objective within VT that behaves well under maximization. We then introduce "socond-order" lower bound, which has provably lower error. Unlike the last two sections, these approximations do not have analogues for arbitrary gradient estimators such DReG score function estimators. For optimization, gradients of the surrogate objective.
compute the complete U-statistic without
reparameterization
Estimator The approrimate complete U-statistic IW-ELBO estimator is
n,m (Vim))
11
Σ (m))
max max(Vs............... (V......, Von) In m.
xE (
This estimator uses the approximation In Em= e ze max.ei..... Um for log-sum-cxp. The following Propo- sition shows that we can compute LA exactly without going over the (m) subsets but instead taking only OnnIn)) time. The intuition is that cach of the log- weights will be maximum clement of some number of size-m subsets, and each such term in the summation for L will be the same. Moreover, we can reason in advance how many times each log- weight will be a maximum.
Proposition 5.1. For any Vran R it holds that
- n
AA n.m (Vi.)) =
6rin Inm,
where b; = -(==i), if E [n (m 1) (and O otherwisc), and [-]: :: [n]] - [m] is permutation s.t., the sequence of wei-weights .. M[n] is non-increasing.
Proof. For (In!), let v (Vap...., V...) We can see that max vs Uyy where the smallest index in s. Thus,
Σ max V = Σ
i=1
where is the number of sets [x]] with minimum index equal to i. The conclusion follows because there are i indices larger than but can take 1 of them only when i E [n (m 1)]. j
To further understand both the computational simplification and the quality of this approximation, con- sider this real example of computing the (non-approximate) complete U-statistic IW-ELBO estimator 2y Suppose that the sampled log-weights are
v (-6034.091, -4351. 43351.335, -4157.236, 55419.201)
Given the (2) sets, we can evaluate the kernel h(v, vy) In(e" + +e*) In on each of them generate the following table:

---

Published in Transactions on Machine Learning Research (02/2023)
(0,, 0/)
h(0,, 04)
6034.091, 4331.3335
4352.028
6034.091, -4157.236)
4157.930
6034.091, -5419.201)
5419.895
4351.335. -4157.236)
4157.930
4351.335,
-5419.201)
4352.028
4157.236, -5419.201
4157.930
Mean
4432.956
At three decimal points of precision, we see that k() = max(e; 01) + In 2 and therefore -4157.930, -4352.028, and -5419.895 each appear (?) times, (7) times, and once, respectively.
-000
-1150
300
objectiver
objestiv
-1200
400
My
standard permutted
-50
-1.00
2000
1000
6000
8000
16930
500
1000
1500
2000
0
1800
2000
3000
2000
5000
iterations
iteratione
iterations
Figure 2: Median envelope of the objective using the permuted-block and standard IW-ELBO estimators for the mushrooms (left), mesquite (center) and electric-one -rrd (right) models. In all cases used = 16 and m 8. For reference there is line segment of length similar to the average objective difference, respectively, 202.08, 0.97, and -391..
5.1 Accuracy and Properties of the Approximation
It is straightforward to derive both upper and lower bounds of the complete U-statistic IW-ELBO estimator At from this approximation.
Proposition 5.2. For any set of log-weights v1.n R", it holds that
(v Lu L. (via) Inm.
(10)
Moreover, the first inequality is strict unless On the other hand, the second inequality is an equality when all weightss equal.
Proof. This is direct application of well-known incqualitics for log-sum-cxp. Let h(v1, ...,, Um) In EM c
and f(i....., um) = max(vi,...... J. Then, for all U1am m
f(umm) h((1m)) S f(11m) + In m.
(11)
To this, write 6 max ( (U,,..., )- Then,
1
< 1 Σ 2 < c
(12)
Eq. (11) follows from applying In to (12). Eq. (10) then follows from (11) and the definitions of LA and AU b
One comment about the approximation quality is in order: in the limit as the variance of the log- weights de- creases, the second incquality in the bounds above becomes tight, and the approximation error of ÉA (vim) approaches its maximum m. This can be seen during optimization when maximizing the IW- IW-ELBO, which tends to reduce log-weight variance (ef. Figure 4).
9

---

Published in Transactions on Machine Learning Research (02/2023)
5.2 Recond-Order Approximation
Based on our understanding of the approximation properties of n,mo we can add correction term to obtain second-order approximation.
Estimator 6. For 2 S S n, the second-order approzimate complete- U-statiiie 1W-ELBO estimator is
n 1--
L (Vi:n) = 2n,m (Vim) + (m)
Σ
8, In(1 Vi)),
(13)
where A = V/--1 - Vg and 6 (mm-2)
This can still be computed in (nn) time and gives a tighter approximation than L m Proposition 5.3. For all vim E R"
f (Vi:n) < L. (V1m)) S L (11=a)-
Moreover, the second inequality is equality evactly when m = 2.
Proof. The first inequality follows directly because the terms in the summation of (13) are positive reals.
For the second inequality, take s E (Iml) and let be the smallest index in s. If one of the (mm-2) sets on which i is the smallest index and + E s, then
1 Le (1 + e -((41-10)) = 6 + c
1
Σ
SE8
m Ei=1e's, ex We finish by applying logarithm both inequalities and If + 1 Sss, s, we know that Le s
1 5m
m the definition of LA and Ly C
In contrast to LAm. the second-order approximation is not U-statistic. However, it is tighter lower-bound of 20 m
Note 5.4. To use the approximations as an objective, we need them to be differentiable, If the distribution of W is absolutely continuous, then the approximations are almost surely differentiable because sort is almost surely differentiable, with Jacobian given by the permutation matrix it represents [cf. Blondel et al. (2020)).
6 Experiments
In this section, we empirically analyze the methods proposed in this paper. We do so in three parts: we first study the gradient variance, VI performance, and running time for IWVI in the "black-box" setting?; we then focus on case where the posterior has closed-form solution, using random Dirichlet distributions; and finally, study the performance of the estimators for Importance- Weightt Autoencoders.
For black-box IWVI, we experiment with two kinds of models: Bayesian logistic regression with 5 different UCI datasets (Dua& & Graff, 2017) using both diagonal and full covariance Gaussian variational distributions, 3 and a suite of 12 statistical models from the Stan example models (Stan Development Team, 2021: Carpenter et al., 2017), with both diagonal (all models) and full covariance Gaussian (10 models") approximating
2That is, VT that uses only black-hox access to In 1n2(2, x) and its gradients.
3That p(y () II Bernoulliyy:: (y; logistic(@ (gg.)tic(TTx.)) for fixed x; Rd and p(0) Ni0; oI) and V In p(0.)))In g()) for 0 ~ (0) with either g(0) A p, diag(w)) q(0) = /V μ,, LLT); optimize over (μ, w) (n, L), with constrained to be positive (via exponential transformation) and L constrainod to be lower triangular with positive diagonal (via softplus transformation). Parameters were randomly initialized prior transformations from fid standard Gaussians.
"The irt-multilevel model diverged for all configurations using a full covariance Gaussian.
10

---

Published in Transactions on Machine Learning Research (02/2023)
distributions. We provide additional information regarding the models in Appendix C. For each model, the variational parameters were optimized using stochastic gradient descent with fixed learning rate for 15 different logarithmically spaced learning rates. We used n = 16 samples per iteration except for the running time analysis, and experimented with m E (2,4,8) Since this is stochastic optimization problem, we ran every combination of model, learning rale, n, and m, using different. random seeds to assess typical performance. We used the reparameterization gradient estimator as the base gradient estimator, and also provide in Appendix D and G (very similar) results for the doublyyrreprrameterieed (DReG) gradient estimator.
0.8
callo of of lntal variaao
matio of abjoetive's vuriante
0.0
0.7
0.8
0.6
0.6
0.3
0.0
0.4
complete-U
aporox
appres. 2nd
mardon
permoited
complete-U
approc
appoor.
randon subs. peraotee
(a) Ratio gradient's total variance.
(b) Ratio of objective's variances.
Figure 3: Ratios of the brace the variance (i.e., the total variance) of different proposed gradient estimators to that of the standard gradient estimator, and objective's variances (b) for the mushrooms dataset (d 96). All ratios are below 1, which indicates variance reduction. The cstimators can be ordered by variance: the complete- U-stattt- estimator and second-order approximation are lowest, followed by the permuted- block and first-order approximation, and finally the random subsets estimator. Since f = 20, we expect the permuted-block estimator to achieve 1 - 20 = 0.95 of the variance reduction; the estimated variance reduction is 91.24% and 95.72% for the objective.
Gradient Variance We first confirm empirically that U-statistics reduce the variance of gradients within IWV1. For each random seed, we performed IWVI using the complete U-statistic 2 for 10,000 iterations. Every 200 iterations, we computed the gradients, given the values the parameters at that time, for each of the alternative gradient estimators: the standard est imtaoor,, the complete U-statistic estimator, its approximations, the permuted-bloot block estimator with 20, and the random subsets estimator with 20 (a number of sets equal to the permuted version). In all cases we used 16 and m 8. For cach gradient estimator č. we estimate the total variance tr( tr(Var(@) using 200 independent gradient samples.
Figure 3 (a) shows the total variance of each estimator as a fraction of that of the standard estimator (that is, the ratio tr(Var[g))/ tr tr(Var((n.m))) for Bayesian logistic regression with the mushrooms dataset. The ratios are between 60% and 70% for all methods, with the random subsets estimator showing the highest variance and the complete U-statistic the lowest. This confirms is possible to reduce gradient variance with U-statistics. Moreover, the estimators can be ordered by their gradients' total variance. The complete U-statistic estimator and the 2nd order approximation have the smallest variance, the permuted- block estimator has slightly higher variance, and the random subsets estimator has the highest variance (but still less than that of the standard estimator). Recall that, according to Prop. 4.2, & = 20 implics that the permuted estimator achieves 95% of the variance reduction provided by the complete-U-staatiste U-statistic IW-ELBO. In this case, we estimated the variance reduction of the permuted- block estimator to be 91.24% that of the complete- Complete-U-statistic U-statistic estimator. We also show the ratio of the objective's variances in Figure 3-(b). Most estimators have a ratio of around 80%, but the permuted-blodd estimator achieves 95. 72% variance reduction provided by the complete U-statistic estimator.
11

---

Published in Transactions on Machine Learning Research (02/2023)
Table 1: For Bayesian logistic regression (a) and Stan models (b), difference in nats of the average objective (higher is better) when trained using the permuted estimator vs. the standard IW-ELBO estimator. The variational distribution is Gaussian distribution using either a full rank covariance matrix (first column) diagonal one (sccond column). The entry is "_" when the model diverged for all configurations, and NaN when it diverged for the specific configuration (other configurations are found in the Appendix).
(a) Bayesian logistic regression models.
(b) Stan models.
permuted standard IW-ELBO
permuted standard IW-ELBO
Dataset
8
Dataset
8
Full Covariance
Diagonal
Full Covariance
Diagonal
ala australian
112.45
4.48
congress
19.80
7.33
3.36
1.38
election88
1133.70
6.94
ionosphere
16.58
0.06
election88Exp
NaN
32.76
mushrooms
202.56
8.69
electric
80.46
4.32
sonar
50.62
0.19
electric-one-pred
-3.45
-3.91
hepatitis
NaN
0.65
hiv-chr irt
283.19
15.84
(6077.03
1.00
irt-multilevel
62.32
mesquite
1.41
2.00
radon
268.98
14.83
wells
-0.03
-0.11
VI Performance Ultimately, goal is to provide a more efficient optimization method. To measure typ- ical stochastic optimization performance, we first took the maximum objective value across learning rates in each iteration to construct the optimization envelope for each method and random seed [ef. Geffner & Domke (2018)]. The purpose the envelope is to eliminate the learning rale as a nuisance parameter since stochastic optimization methods are very sensitive to learning rate, and one common benefit of variance reduction is to allow a larger learning rate. Then, for each method we used the median envelope across the random soods measure of its lypical optimization behavior metric for each method computed the average objective value (of the median envelope) across iterations up to 10,000 iterations, 5 excluding the first 50 iterations, which were highly noisy and sensitive to initialization. This is useful summary metric to measure the tendency of one method to "stay ahead" of another (see the examples Figure 2). Agrawal et al. (2020) found a similar metrie effective for learning rate selection.
iterations. Examples can be seen in Figure 2. As final
Table 2: Times for 1000 iterations of optimiza- tion with different estimators on the mushrooms dataset with 24, 12, averaged over 100 trials.
Table 1 shows the average objective difference between the permuted-block and standard IW-ELBO estimators for 8, with positive numbers indicating better per- formance for permuled--lock. We focus on permuted- block here because it consistently achieves an excellent tradcoff between variance reduction and running time. Appendix we present similar results for two addi- tional methods- 2nd order approximation and the permuted-block estimator with DReG as the base gra- dient estimator- and for different values of m; in Ap- pendix G we show the median envelopes themselves for many combinations of models, methods, and m. The ex- amples in Figure were selected to show cases where the difference is big (left), small (center), and negative (right); contextualize our summary metric, we also added ref- erence vertical bar showing an iteration where the difference between the two envelopes approximately
Time (s)
Method
Mean
Std
Lu112 standard IW-ELBO
5.47
0.04
Cha.12 complete U
1573.27
2.12
Csz024 random subsets
6.49
0.09
$20 permuted block
6.45
0.09
C approx.
5.25
0.02
2 approx. 2nd order
5.54
0.04
c For datasets, such sonar, observed early convergence by visual inspection and computed metric only up to that point. See Figures Appendix G.
12

---

Published in Transactions on Machine Learning Research (02/2023)
equal to the average objective difference. These results make it clear that the permuted-block estima- for improves the convergence of stochastic optimization for VT across range of models and settings. In electric-one-pret permuted-blooe was consistently worse, but we verified that it still had lower-variance gradients; we speculate this is an unstable model where higher variance gradients help escape local optima.
Running Time Table 2 shows the times required to complete 1000 iterations of optimization with different estimators for Bayesian logistic regression with the mushrooms data sel, averaged over 100 trials. Here we used n = 24 and m = 12, which makes a challenging setting for the complete U-statistic estimator, because there are (22 2,704,156 sets. As expected, the complete U-statistic orders of magnitude slower. The approximations are faster than the standard estimator because the smallest m 1 log- oo-weightss do not contribute to the objective, and thus their gradients are not needed. The permuted-block estimator incurs extra cost of less than ms per iteration compared to the standard IW-ELBO estimator for this model (a 18% increase). However, the increased time only depends on m, n, and f, and not on the model. Even for a very complex model, we would expect the extra time for these settings to be on the order of order of 1 ms per iteration, and be negligible compared to other costs. For example, for the irt, the standard estimator took 16.62s (0.11), while the permuted-block estimator took 17.54s (0.12), ie., a 5% increase.
Incomplete U-Statistics and Approximations Proviously, we analyzed the methods by com- paring them to the standard IW-ELBO estimator. In this part we will use the complete U- statistics as baseline: given a realization of log-weights "h....,"n. we measure the differ- ence between the objective value assigned by the complete U-statistic and the alternatives. For this experiment, we will use mushrooms. In Figure 4 we plot the difference measured in nats as a function of the iteration step. From that plot (especially the insct), it is clear that the approximations are underostimators. It is also interesting to see the approximations and the incomplete U-statistic being complementary: as the optimization progresses, the error of the approx- imations increases, but the error made by the in- complete U-statistics decreases. We expected this result because the variance of the log- weights de- creases with the optimization. (The soper-bound of Eq. (10) is achieved when all iú are equal; but this is exactly the case when all the incomplete U-statistics coincide.)
16,
8 and the Bayesian logistic regression dataset
60
siaudard IW
0.50
50
Approx:
0.25
apprexx 2nd
40
Zandlom subests
0.00
30
permatted
-0.25
nola
20
-0.50
-10
Dirichlet Experiments We conducted experi- ments with random Dirichlet distributions as de- scribed in (Domke & Sheldon, 2018). The goal was twofold. First, this is setting where exact inference evaluate IWVI with different istiors on the accuracy of posterior inference di- rectly, instead of using the IW-ELBO as a proxy. Secondly, this is a simple setting to demonstrate that the optimal value of m is often strictly between and n, which the regime in which our variance reduction methods are useful (all but the approxi- mations coincide when E {1,m)). We again used SGD with 15 different learning rates and selected, for each configuration, the learning rate that achieved the best mcan objective after 10k itcrations. Wo optimized each configuration using, for this experiment, 100 different random seeds. We estimated the accu- racy of the approximation by computing the distance (error) between the distribution's covariance and the estimated covariance of the learned approximation. Figure shows the error as a function of m for different values of n when using the standard IW-ELBO estimator for random Dirichlet with 50 parameters. The figure shows that. the optimal m increases with n, but slowly. Figure shows similar results for other esti-
0
2000
4000
6000
5000
10000
itenatonns
Figure 4: Difference between estimated value using any of the methods and that of the complete U- statistic, in nats, for the mushrooms dataset. (25th and 75th percentiles shown with dashed lines.) As op- timization progresses, the error of the incomplete U- statistics decreases, but the error of the approximation increases. The inset shows the permuted and both ap- proximations in region that is 0.5 nats of the target value.
is possible,
13

---

Published in Transactions on Machine Learning Research (02/2023)
mators: permuted, DReG, and permuted-DReG. In all cases, confirm that, for this model, the optimal m lies strictly between and n. We provide in Appendix E additional details.
16-5
permuted
1e-5
DReG
16-5
permuted-DReG
2,5
2,5
2.3
2.3
arro
2.1
2.1
128
256
2.5
2.0
1.9
2
2
Figure 5: Distance between the covariance of random Dirichlet distribution with 50 parameters and the covariance of its approximation as function of m for different values of after training using the permuted (left), DReG (center) or permuted-DReG (right) estimators.
6.1
Autoencoders
To evaluate the performance of the proposed methods on IWAEs, we trained IWAEs on different datasets: MNIST, KMNIST, FMNIST, and Omniglot. We compare the standard IW-ELBO estimator and DReG estimators to their permuted versions, i..., the permuted and permuted-DReG estimators. We also evaluate the sccond- order approximation to the complete-U-staaistie U-statistic estimator. We trained each combination of dalaset, method, and value of m using five different random seeds, and the optimization was run for 100 epochs using Adam (Kingma & Ba, 2015).
In Figure 6, we present the final lesting objective for different values of m (using n = 50 all cases) for the KMNIST dataset, and show results for the rest of the datasets in Figure in the Appendix F along with further details on the experiments. The figure shows that the permuted versions consistently improved the base versions, i.c., the permuted estimator improves over the standard-IW estimator in the same way as the permuted-DRReG estimator improves over the DReG est.mmaoor Additionally, we can see that the second-order approximation outperforms the permuted estimator for small values of m. However, as m increases, the permuted estimator takes the lead, which expected since the approximation crror grows with m.
Jrmnist
knnnist
10
kmnnist
201.5
X
IW abjective
standard
-20255
X
X
appeor. 2nd
permuted
200.0
X
DReG
permuted DRtG
-00.5
-20..
method
method
method
Figure 6: Objective's distribution for KMNIST with n 50 and different combinations of methods and m.
We also compared the total wall-clock time required to complete the optimization with different estimators in Figure 9 in the Appendix. can be seen that there is not significant time increase for using our proposed methods.
14

---

Published in Transactions on Machine Learning Research (02/2023)
7 Related and Future Work
Gradient variance reduction is an active topic in VI because of its impact on stochastic optimization. Our complete- and incomplete-U-statiste methods are complementary to other variance reduction techniques: they are compatible with different base estimators, including the Doubly Reparameterized Gradient Estima- tor (DReG) Tucker et al. (2018) and the generalization of Bauer & Mnih (2021). Another broad approach variance reduction the use of control variates (Miller al., 2017; Mnih Gregor, 2014; Ranganath al., 2014; Geffner & Domke, 2018; 2020). In the case of IWVT, the control variates of Mnih & Rezende (2016) and Liévin al. (2020), which designed for the score function estimator, could work base estimator from which U-statistic can be built. We leave its empirical evaluation for future work.
weighted estimators are also being used for the Reweighted Wake-Sleep (RWS) procedure (Born- schein & Bengio, 2015; Le et al., 2020) and its variations (Dieng & Paisley, 2019; Kim et al., 2020). Given the connection between the gradient estimators RWS and that of the IW-ELBO [see Kim et al. (2020)), these estimators could be potentially improved by using the ideas of the complete- and incomplete---statittii U-statistic methods.
The numerical approximations of Section 5 follow a different principle of approximating the objective; it is an open question if such an approximation can be used in conjunction with other variance reduction methods. Interestingly, the first-order approximation expresses the objective as a convex combination of the ordered log- weights (minus a constant), which has a form similar to the objective presented in Wang et al. (2018), albeit with different coefficients.
It would be an interesting future line of work to extend the order of Proposition 4.1 a partial order of random variables in the sense of Mattei & Frellsen (2022).
Nowozin (2018) introduced Jackknife-VT (JVT), which uses complete-U statistics to reduce bias instead of variance. In Appendix A we briefly discuss possible applications of our methods to JVI.
8 Conclusion
We introduced novel methods based on U-statistics to reduce gradient and objective variance for importance- weighted variational inference, and found empirically that the methods improve black-box VT performance and IWAEs training. We recommend using the permuted-block estimator in any situation with r > I replicates: it never increases variance, and can be tuned based on computational budget to achieve any desired fraction of the possible variance reduction. In practice, a 95% fraction of possible variance reduction very low cost. The approximations of Section 5 are extremely fast and provide substantial variance reduction, but are not universally better than the standard estimator because they introduce some bias that can hurt performance, especially in easier models near the end of optimization.
can be achieved
Acknowledgments
This material is based upon work supported by the National Science Foundation under Grant Nos. 1749854 and 1908577. JB would like to thank Tomás Geffner and Miguel Fuentes for their helpful discussions.
15

---

Published in Transactions on Machine Learning Research (02/2023)
References
Abhinav Agrawal, Justin Domke, and Daniel Sheldon. Advances in black-box VI: Normalizing flows, im- portance weighting, and optimization. In Advances in Neural Information Processing Systems (NeurIPS), pp. 1-8, 2020.
Matthias Bauer and Andriy Mnih. Generalized doubly reparameterized gradient estimators. In International Conference on Machine Learning, pp. 738-74777 747. PMLR, 2021.
Eli Bingham, Jonathan P. Chen, Martin Jankowiak, Fritz Obermeyer, Neeraj Pradhan, Theofanis Karaletsos, Rohit Singh, Paul Szerlip, Paul Horsfall, and Noah D. Goodman. Pyro: Deep Universal Probabilistic Programming. Journal of Machine Learning Research, 2018.
David M Blei, Alp Kucukelbir, and Jon D McAuliffe.
Variational inference: A review for statisticians. Journal of the Americon statistical Association, 112((18):859 877, 2017.
Gunnar Blom. Some propertics of incomplete U-statistics. Biometrika, 63(3)5573-555 580, 1976.
Mathicu Blondel, Olivier Teboul, Quentin Berthet, and Josip Djolonga. Fast differentiable sorting and ranking. In Hal Daumé III and Aarti Singh (eds.), Proceedings of the 37th International Conference Machine Learning, volume 119 of Proceedings of Machine Learning Research, pp. 950-959. PMLR, 13-18 Jul 2020. URL https://prroceednngs mlr wtt.:/in.......eee//////////////
html.
Jörg Bornschein and Yoshua Bengio. Reweighted wake-sleep. In ICLR, 2015.
Yuri Burda, Roger Grosse, and Ruslan Salakhutdinov. Importance Weighted Autoencoders. In ICLR, 2016.
Bob Carpenter, Andrew Gelman, Matthew D Hoffman, Daniel Lee, Ben Goodrich, Michael Betancourt, Marcus Brubaker, Jiqiang Guo, Peter Li, and Allen Riddell. Stan: A probabilistic programming language. Journal of statistical software, 76(1), 2017.
Chih-Chune Chang and Chih-Jen Lin. Libsvm: a library for support vector machines. ACM transactions on intelligent systems and technology (TIST), 2(3):1-27, 2011.
Tarin Clanuwat, Mikel Bober-Irizar, Asanobu Kitamoto, Alex Lamb, Kazuaki Yamamoto, and David Ha. Deep learning for classical japanese literalure, 2018.
Chris Cremer, Quaid Morris, and David Duvenand. Reinterpreting inportance-weighted autoencoders. arXiv preprint Xiv: (eee:104400010 2017.
Adji B Dieng and John Paisley. Reweighted expectation maximization. arXiv preprint arXiw: 05850. 2019.
Justin Domke and Daniel Sheldon. Importance weighting and variational inference. In Proceedings of the 32nd International Conference on Neurra Information Processing Systems, pp. 4475-4484, 2018.
Justin Domke and Daniel R. Sheldon. Divide and couple: Using Monte Carlo variational objectives for posterior approximation. In Advances in Neural Information Processing Systems 32: Annual Conference on Neural Information Processing Systems 2019, NeurIPS 2019, December 8-14, 2019, Vancouser, BC, Canada, pp. 338-347, 2019.
Dheeru Dua and Casey Graff. UCI machine learning repository, 2017. URL http:///rrchive.iie uci edu. ml.
Axel Finke and Alexandre Thiery. On importance-weighted arXiv:1907.10477 2019.
autoencoders.
arXie preprint
Michael C Fu. Gradient estimation. Handbooks in operations research and management science, 13:575-616 2006.
16

---

Published in Transactions on Machine Learning Research (02/2023)
Tomas Geffner and Justin Domke. Using large ensembles of control variates for variational inference. In Proceedings of the 32nd International Conference on Neural Informalion Processing Systems, pp. 9982- 9992, 2018.
Tomas Geffner and Justin Domke. Approximation basod variance reduction for reparameterization gradients. Advances in Newral Information Processing Systems, 33, 2020.
Andrew Gelman and Jennifer Hill. Data analysis using regression and multtbeeve//hierrcchiaal models. Cam- bridge university press, 2006.
Paul R Halmos. The theory of unbiased estimation. The Annals of Mathematical Statistics, 17(1):34-43, 1946.
Wassily Hoeffding. A class of statistics with asymptotically normal distribution. Annals of Mathematical Statistics, 19:273-322 325, 1948.
Dongha Kim, Jaesung Hwang, and Yongdai Kim. On casting importance weighted autoencoder an em algorithm to learn deep generative models. In Silvia Chiappa and Roberto Calandra (eds.), Proceedings of the Twenly Third International Conference on Arlificial Intelligence and Stalistics, volume 108 of Proceedings of Machine Learning Research, pp. 2153 2163. PMLR, 26 28 Aug 2020. URL https:// proceedings.mk. mlr. wat.aags.mmmm....e/eeeee0/0e kiml...html.
Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. In ICLR (Poster), 2015.
Diederik P Kingma and Max Welling. Auto-cncoding variational bayes arXiv preprint arXiv:1312.611 2013.
Jack PC Kleijnen and Reuven Y Rubinstein. Optimization and sensitivity analysis of computer simulation models by the score function method. European Journal of Operational Research, 88(3):413 1171112--777-- 1996.
Brenden M Lake, Ruslan Salakhutdinov, and Joshua B Tenenbaum. Human-level concept learning through probabilistic program induction. Science, 350(6266):1332 1338, 2015.
Tuan Anh Le, Maximilian Igl, Tom Rainforth, Tom Jin, and Frank Wood. Auto-encoding sequential Monte Carlo. In 6th International Conference on Learning Representations, ICLR 2018, Vancowver BC, Canada, April 30 May 3, 2018, Conference Track Proceedings, 2018.
Tuan Anh Le, Adam R. Kosiorek, N. Siddharth, Yee Whye Teh, and Frank Wood. Revisiting reweighted wake- sleep for models with stochastic control flow. In Ryan P. Adams and Vibhav Gogate (eds.), Proceedings of The 35th Uncertainty Artificial Intelligence Conference, volume 115 of Proceedings of Machine Learning Research, pp. 1039 1049. PMLR, 22 25 Jul 2020. URL https://ppti //ppce/e/n/occco.......o mlr. ppwas/v115/1e2000 html.
Yann LeCun, Corinna Cortes, and CJ Burges. Mnist handwritten digit. database. Available: http://yann, http///nnnnneeun..o www.eoooomooiooe /ee///dd///m///ee /misst, 2, 2010.
ATT Labs [Online].
Laurence Lock Lee. U-statistics. Theory and Practice. CRC Press, 1990. ISBN 978135140580
Valentin Liévin, Andrea Dittadi, Anders Christensen, and Ole Winther. Optimal variance control of the score- function gradient estimator for seoootancc-woighted bounds. Advances in Neural Information Processing Systems, 33:16591- 16602, 2020.
David J Lunn, Andrew Thomas, Nicky Best, and David Spiegelhalter. Winbugs a bayesian modelling framework: concepts, structure, and extensibility. Statistics and computing, 10(4):325-335 2000.
Chris J. Maddison, Dieterich Lawson George Tucker, Nicolas Heess, Mohammad Norouzi, Andriy Mnih, Arnaud Doucet, and Yee Whye Tch. Filtering variational objectives. In Advances in Neural Information Processing Systems 30: Annual Conference on Neural Information Processing Systems 2017, December 4-9, 2017, Long Beach, CA, USA, pp. 6573 05733,653,441 2017.
17

---

Published in Transactions on Machine Learning Research (02/2023)
Nathan Mantel. 232. note: Assumption estimators using U statistics and relationship to the Jackkmife method. Biometrics, pp. 567-571, 1967.
Pierre-Alexandre Mattei and Jes Frellsen. Uphill roads to variational tightness: Monotonicity and monte carlo objectives. arXiv preprint Xrx::2201. 10989, 2022.
AC Miller, NJ Foti, A D Amour, and Ryan P Adams. Reducing reparameterization gradient variance. Advances in Newral Information Processing Systems, 2017.
Andriy Mnih and Karol Gregor. Neural variational inference and learning in belief networks. In International Conference Machine Learning, pp- 1791-1799. PMT.R, 2014.
Andriy Mnih and Damilo Rezende. Variational inference for monte carlo objectives. In International Con- ference on Machine Learning, pp. 2188 2196. PMI.R, 2016.
Christian A. Naesseth, Scolt. W. Linderman, Rajesh Ranganath, and David M. Blei. Variational sequential Monte Carlo. In International Conference on Artificial Intelligence and Statisties, AISTATS 2018, 9- 11 April 2018, Playa Blanca, Lanzarote, Canary Islands, Spain, volume 84 of Proceedings of Machine Learning Research, pp. 868-9777 PMLR, 2018.
Sebastian Nowozin. Debiasing evidence approximations: On weighttd autoencoders and Jack- knife variational inference. In International Conference on Learning Representations, 2018
Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, Jamos Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala. Pytorch: An imperative style, high-performance deep learning library. In H. Wallach, H. Larochelle, A. Beygelzimer, F. d' d'Alcée-Buc, E. Fox, and R. Garnett (eds.), Advances in Neural Informa- tion Processing Systems 32, pp. 8024 8035. Curran Associates, Inc., 2019. URL
neurips. hiags-llbbrryyyre pdf.
John C Platt. Sequential minimal optimization: fast algorithm for training support vector machines. Tn Advances in Kernel Methods-Support Vector Learning, 1999.
Tom Rainforth, Adam Kosiorek, Tuan Anh Le, Chris Maddison, Maximilian 1gl, Frank Wood, and Yee Whye Tch. Tighter variational bounds are not necossarily better. In International Conference on Machine Learning, pp. 4277-4285. PMLR, 2018.
Rajesh Ranganath, Sean Gerrish, and David Blei. Black box variational inference. In Artificial intelligence and statistics, pp. $14- PMLR, 2014.
Danilo Jimenez Rezende, Shakir Mohamed, and Daan Wierstra. Stochastic backpropagation and approximate inference in deep generative models. In International conference on machine learning, pp. 1278 1286. PMLR, 2014.
Lawrence K Saul, Tommi Jaakkola, and Michael I Jordan. Mean field theory for sigmoid belief networks. Journal of Artificial Intelligence Research, 4:61-76, 1996.
Stan Development Team. example-modee models.
Stan Example models, 2021. URL https://gtthub..oomsttaar dev/
George Tucker, Dieterich Lawson, Shixiang Gu, and Chris J Maddison. Doubly reparameterized gradient estimators for Monte Carlo objectives. In International Conference on Learning Representations, 2018.
Aad W van der Vaart. Asymptotic statistics, volume 3. Cambridge University Press, 2000.
Dilin Wang, Hao Liu, and Qiang Liu. Variational inference with (ail-adaptive (ariaaptiivvddiieegeeeeee f-divergence. Advances in Newral Information Processing Systems, 31, 2018.
Han Xiao, Kashif Rasul, and Roland Vollgraf. Fashion-mman novel image dataset for benchmarking machine learning algorithms, 2017.
18

---

Published in Transactions on Machine Learning Research (02/2023)
A Experiments with Jackknife
The relation ween the jackknife estimator and complete U-statistics was made explicit early on by Mantel (1967). Recently, Nowozin (2018) used the jackknife estimator as way to diminish the bias in IW-VI, proposing jackknife VI (JVT). Using the notation of Section 3, the jackknife estimator is
(V = Σ c(n, rn,, -
(14)
j=0
where nn, is the complete U-statistic IW-ELBO estimator, and the r.jr are the Sharot coefficients (ef. Nowozin (2018)].
In the original version (14), it evaluates collection of complete U-statistics with ranging from to r. However, there need to constrain in that way, i.e, we can instead compute the following estimator
(Vin)) = Σ c(m,r,j)4 n- (Vin), for r m S n,
because the bias is function of m. This means that once is fixed, pick the number of independent samples 2 to reduce the variance of the estimation.
For our experiment, we optimized a variational approximation to the posterior of the mushrooms dataset as in Section 6. We used the complete-U. U- statistic IW-ELBO estimator for optimization (n' 16 and m 8), and we choose the configuration with the highest final bound.
methhd
-100
Compdarie-U
apyenr. 2nd
permonsed (2 (00)
-180
We evaluated the trained model using the Jackknife estimator with n = 24, r = and m = 8. For the inner estimator we used the complete-U-stttiisii IW-ELBO estimator, a variation of the permuted- block IW-ELBO estimator" with = 20 and with t= 100, and the second order approximation. Fig- ure 7 shows that, when using the permuted estima- tor with & 20, the increased variance gets trans- lated into an increased variance in the final estima- tion. However, it can be reduced by increasing the number of permutations to l = 100. In the follow- ing table, we show the time taken compute the Jackknife estimator without accounting for the time
U-statistic
200
-210
complete-U
approx. 2nd
permatte (f - 23) permuted (t- 100)
Figure 7: Distributions of the objective using the Jack- knife estimator Li on an approximation of the pos- terior of the mushrooms dataset, using different esti- mators.
of building the index s
Method
Mean time (ms)
Std
complete-U U
23.98
2.28
approx. 2nd
0.71
0.05
permuted (8 (6=30 20)
0.69
0.02
permuted (f = 100)
0.86
0.03
In this case, we observed that the alternatives are approximately 30 times faster than using the complete- U-statistic.
6Since not integer multiple of m reduced the total number of scts Em T We pre-computed the index set for the complete-U-snee U-statistic, which this requires 735471 346104 1081575 sets, taking 1.34 seconds.
19

---

Published in Transactions on Machine Learning Research (02/2023)
B Additional Theoretical results
In this section, we apply result of Halmos (1946) to the estimation of the IW-ELBO. Subject to certain conditions, the estimator C has the smallest variance of any unbiased estimator of the IW-ELBO. The technical conditions are needed to define the class of "unbiased estimators" as ones that are unbiased for all log-weight distributions in non-trivial class.
Proposition B.1. Let Ey[-] and Varel.1 denote expectation and variance with respect to log- weights Vi...., Vn drown independently from distribulion F, and let (F)F) Er in ex)] be the IW- ELBO with log-weight distribution F. Let Ff denote the set of distributions supported on finite subset of R. Suppose 4 is any estimator such that [p(D(Vim)] (m) for all F E F. Then,
Varr (Visn)] S Varr
whenever the lotter quantity is defined, for any distribution F on the real numbers (up to conditions of measurability and integrability).
Proof. The result direct application Theorem Halmos (1946).
For IW-ELBO estimation, the conditions are rather mild: we expect IW-ELBO estimator to work for generic log-wer weight distributions. For gradient estimation, we take the conclusion lightly, because gradient estimators often use specific properties the underlying distributions, such as having a reparameterization.
B.1 Additional Proofs
In this section, provide a proof of Proposition 3.3. We first need define quantity similar to Def- inition 3.1. Recall, from the statement of the Proposition, that y: Rdz Rdo, and let gi denote its ith component.
Definition B.2. Let Zi,..., . Zem i.i.d. drawn from qø and 1 sis dg. For 0 S cs m, take s, ( (12ml) with Is sns' = Using g from Proposition 3.3, define
(0) = Cov (Z.,... Zam) 9e (z Zun.
which depends only and not the particular and s'.
We can now proceed with the proof.
Proof of Proposition 3.3. For S i s do, using Eq. (4) and (6), it follows from Theorem 5.2 of Hoeffding (1948) that
Vari < Sm = Var
From the definition of the covariance matrix, we get
de Va
du
tr tr(Var(@u,,.)
bi <
Varl(Gn.m)l tr(Varl@n..)).
1=1
Using again that G and &n.m are unbiased, that is, ElGu = E E(Gn,m), then
do (Var(( + E( (m
E g
n.m 112
1113)
(Var((G.,.)))) + (E((Gn,m))))) 119 12.
20

---

Published in Transactions on Machine Learning Research (02/2023)
C Dataset description
We provide a brief description of the datasets and models used for the experiments. The models used for Bayesian logistic regressions were taken from the UCI Machine Learning Repository Dua & Graff (2017). The rest of the models are part of the Stan Example models Stan Development Team (2021); Carpenter et al. (2017).
For the dataset used for Bayesian logistic regression, whenever there was a categorical variable with k categories, dummified it by creating k- dummies variables. Additionally, for the aia dataset, continuous variables were discretized into quintiles following the work of Platt (1999). However, since unable to find the file describing the actual process used for the discretization, some discrepancies remained.
Table 3: Description of datasets/models.
Num. of variables
Num. of records
Name
Comments
First 1605 instances of the Adult Data Set, following LIBSVM Chang & Lin (2011), + discretized continous and dummified.
ala australian
105
1605
35
690
From UCI + dummified.
ionosphere
35
351
From UCT
mushrooms
96
8124
From UCI + dummified.
sonar
61
208
From UCI
congress
343
Gelman & Hill (2006) Ch. 7
election88
95
2015
Gelman & Hill (2006) Ch. 19
election88Exp
96
2015
Gelman & Hill (2006) Ch. 19
electric
100
192
Gelman & Hill (2006) Ch. 23
clectric-ono-pred
192
Golman & Hill (2006) Ch. 23
hepatitis
218
288
WinBUGS Lunn et al. (2000) examples
hiv-chr
173
369
Gelman & Hill (2006) Ch.
irt
501
30105
Gelman & Hill (2006) Ch. 14
irt-multilevel
604
30015
Gelman & Hill (2006) Ch. 14
mesquite
3
46
Gelman & Hill (2006) Ch. 4
radon
88
919
radon-chr from Gelman & Hill (2006) Ch. 19
wells
3020
Gelman & Hill (2006) Ch. 7
MNIST
784
60000 + 10000
LeCun et al. (2010) Mashion-MNIST Xiao et al. (2017)
FMNIST
784
60000
10000
KMNIST
784
60000
10000
Kuzushiji-MNIST Clanuwat al. (2018)
Omniglot
784
24343
8070
Lake et al. (2015) from Burda et al. (2016)
D Pairwise comparison
In this section we present the mean difference of the medians of the envelopes as described in Section 6. We compare the methods that used the reparameterized gradients as based gradient estimator, i.e., the permuted- block estimator and the 2nd order approximation, to the standard IW-ELBO estimator. Additionally, we compare the standard IW using DReG as based gradient estimator with a version of the permuted-bles that uses the DReG base gradient estimator, namely, the permuted DReG.
Interestingly, in the settings presented in Table 7, only the proposed methods, i.e., the complete-U U statistic with its two approximations, the boonuted-blockk, and the random subsets, converged at some point. All the other methods diverged, which explains why we cannot. compute the difference.
21

---

Published in Transactions on Machine Learning Research (02/2023)
Table 4: Bayesian logistic regression models using Gaussian approximation with a covariance matrix of full rank. Difference in nats of the average objective (higher values are better).
permuted standard 1W
approx. 2nd
standard 1W
permuted DReG DReG
Datasct
m
2
4
8
4
8
2
4
8
ala
45.36
100.56
112.42
47.53
105.30
122.34
27.30
111.74
119.51
australian
1.31
2.61
3.36
1.07
2.37
3.22
1.45
1.87
3.94
tonosphere
3.89
13.17
16.58
4.11
13.55
17.55
4.34
15.74
17.91
mushrooms
64.46
145.85
202.56
67.28
153.58
214.31
93.45
186.01
179.02
30.15
61.09
50.62
32.94
63.34
54.14
27.99
69.54
90.86
Table 5: Bayesian logistic regression models using a diagonal Gaussian approximation. Difference in nats of the average objective (higher values are better).
permuted standard IW
approx. 2nd standard IW
permuted DReG: DReG
Dataset
m
2
8
4
ala
1.54
4.01
4.48
1.49
4.08
4.45
1.40
12.67
12.86
australian
0.02
1.00
1.38
0.05
0.96
1.28
-0.07
0.06
1.43
ionosphere
-0.10
-0.10
0.06
-0.12
-0.21
0.00
-0.12
-0.08
0.31
mushrooms
1.88
2.76
8.69
1.74
3.30
9.16
1.94
4.69
8.50
sonar
0.03
-0.15
0.19
-0.02
-0.18
0.15
0.03
-0.28
0.21
Table 6: Stan models using diagonal Gaussian approximation. Difference in nats of the average objective (higher values are better).
permuted standard IW
approx. 2nd
standard IW
permuted DReG DReG
Dataset
2
4
8
2
4
8
2
4
8
congress
2.50
4.61
7.33
2.89
4.76
7.63
2.37
4.68
7.02
election88
0.12
2.66
6.94
0.12
2.84
7.06
0.10
2.67
6.83
election88lixp
0.82
98.52
32.76
4.73
117.78
55.27
-1.89
100.09
32.53
electric
0.26
1.53
4.32
0.16
1.54
4.52
0.24
1.56
4.63
electric-one-pred
0.66
-0.77
-3.91
0.74
-0.76
-4.38
0.69
-0.77
-3.93
hepatitis
0.90
-0.06
0.65
2.06
156.53
1.86
-0.30
0.92
0.69
hiv-chr
0.16
2.03
15.84
0.34
2.12
21.74
-0.08
1.45
12.91
irt
0.19
0.80
1.00
0.15
0.72
0.93
0.11
0.61
1.40
irt-multilevel
35.69
43.79
62.32
29.74
48.20
53.64
34.66
50.26
55.22
mesquite
0.20
0.58
2.00
-0.06
0.28
1.74
-0.29
0.39
1.99
radon
7.88
5.79
14.83
7.85
8.91
65.49
8.16
7.56
60.92
wells
-0.02
0.01
-0.11
-0.20
-0.30
-0.35
-0.02
-0.04
-0.14
22

---

Published in Transactions on Machine Learning Research (02/2023)
Table 7: Stan models using a full covariance Gaussian approximation. Difference in nats of the average objective (higher values are better).
permuted standard IW
арргоx. 2nd standard IW
permuted DReG DReG
Dataset
772
8
2
4
2
congress
11.62
12.02
19.80
12.33
12.57
20.46
13.55
13.11
20.96
clection88
NaN
1785
1133
NaN
2494
2170
NaN
1776
1116
election88Exp
NaN
NaN
NaN
NaN
NaN
NaN
NaN
NaN
NaN
electric
NaN
-38.02
80.46
NaN
-77.53
89.06
NaN
-43.16
34.91
electric-one-pred
-1811
-4.73
-3.45
-3.18
-4.77
-4.37
-.779
-4.72
-3.46
hepatitis
NaN
NaN
NaN
NaN
NaN
NaN
NaN
NaN
NaN
hiv-chr
NaN
NaN
283.19
NaN
NaN
325.79
NaN
NaN
NaN
irt
17793
20064
16077
19399
22000
17686
NaN
NaN
NaN
mesquite
2.57
1.20
1.41
2.43
0.95
1.19
2.67
0.53
0.74
radon wells
NaN
1150
268.98
NaN
1316
303.83
NaN
11675
269.26
0.02
0.07
-0.03
-0.29
-0.31
-0.29
0.04
0.02
-0.02
E
Random Dirichlet experiment
We follow Domke & Sheldon (2018) for the Random Dirichlet experiment. For a randomly-sampled Dirichlet Distribution with 50 parameters, we approximate using a (50 - 1)-dimensional Gaussian distribution parameterized with full rank covariance matrix, with its domain constraine to the simplex using Py PoToc''s distributions (Paszke al., 2019).
We optimize each configuration using 100 different random seeds. We select the learning rate that achieved the highest objective among all learning rates that converged for all seeds. For each seed, compute the Frobenius norm between the empirical covariance of the approximating distribution and that of the theoretical distribution (the error). The distribution of this error is shown in Figure 1 and 5. We had to exclude eight outliers with errors greater than 10 and up to 0.5. Interestingly, those outliers used either the DReG or permuted-DReG estimators.
F
VAE details.
For the variational autoencoders, we used, for all datasets, the architecture used by Burda et al. (2016). We trained each configuration for a fixed number of epochs (100) using Adam (Kingma & Ba, 2015) with a learning rate of 10-4, In all cases, we used a batch size of 500, and latent variable of dimension 50, while taking 50 samples. Datascts
taken from PyTorch, except for the Omniglot, for which used the construction provided by Burda el al. (2016). We evaluated using the slandard IW-ELBO estimator, regardless of the estimator used for the optimization.
To get consistent wall-clock time measurements, we trained only using CPU on dedicated servers, with disabled hyper-threading and single task per core. Additionally, we used set set_flush fos_fffss___orrm_nn denormal to avoid creating denormal numbers because some estimators create many of such numbers (especially DReG-like estimators), having a substantial negative impact on performance. Our implementation of DReG is based Pyro's (Bingham et al,, 2018) nottyee-integrated implementation. We are not aware of PyTorch imple- mentation without the extra time penalty.
In the following plots, we provide the objective distribution for all Mataset/methodd /mm configurations and the distribution of the wall-clock time.
23

---

Published in Transactions on Machine Learning Research (02/2023)
knunist
lemnist 25
201.5
IW abjective
-202.0
41
202.5
x
x
203.5
204.0
-244.25
-211.50
IW abjcesive
244.75
X
-245.00
245.25
X
245.50
slamdard IW
permuted
omniglot
omniglot
omniglot
appror. 2nd
DReG
148
Pemmatot-DReC
IW abjeatinn
X
-152
x
154
156
X
106.50
106.755
TW abjectiv
107.25
X
107.75
X
108.25
method
method
method
Figure 8: Distribution of the objective for different combinations datasets, methods and m, using 50 samples.
24

---

Published in Transactions on Machine Learning Research (02/2023)
40000
42500
40000
total time (6)
57500
35000
32500
90000
x
27500
x
X
Emnist
450000
42500
40000
tosal tibe (6)
37500
32500
27100
slamdard IW
permated
ornniglot
omnigiot
cmaniglor
approt. 2nd
DReG Pemmatot-DReC
18000
16000
tatal time [8]
14000
13000
12000
11000
X
45000
40000
totel tiae (8
37500
35000
32500
27500
method
method
method
Figure 9: Distribution of the time taken run epochs for different combinations of datasets, methods and m, using n 50 samples.
25

---

Published in Transactions on Machine Learning Research (02/2023)
G Figures of median envelope
For some of the methods presented in the paper, we compute the median envelope during training as described in Section 6.
Zeagao
-40
-40
-193
258
130
230
150 104E
-140
253
-1290
1
-12.00
-1560
-12.00
30
11.00
11/00
-10
-20
-30
1200
11000
130
-12993
-150
114.0
brveskllond 11:01
113.07
-1250
-1220
-1120
-1250
-1000
-15500
-11294
-15790
200 2010
380 1800
-00
180
200
300 3000
180
1800
1300
200
1009
1000
20 13
104
-20
-2X
-30
-220
3001
2010
180
7.08
1800
200
-35
-30
-20
114 130
30 18
101 X00
30 14 108
3.1
180
110
-120
-1000
-130
200
3000
1e000-ieeee pad-
4b000-0a pad
-300
-100
-190
-230
500 130 10K 2019
50 300 196 3000
53 106 1500
200
30 300 13N 3100
brgettt 2
-20
3.4
1il 130
-20
30
1.0
(00
-30
540 20
3.3 1006
-2393
-216
-5311
-210
-2312
-210
-334
-311
-254
-244
-2393
250
130
-2.0
250
236
750
3000
-254
750
104
-2100
223 300
2zudio
Figure 10: Median envelope for models when using diagonal Gaussian as approximating distribution for the estimators complete-U, permuted and the standard 1W.
26

---

Published in Transactions on Machine Learning Research (02/2023)
-12
-13
-225
-150
350
750
100
-15
550
300
750 3400
130
-1.00
-10.0
-1963
-200
-1013
-1511
-10
-2000
-1000
1.30
-1903
-220
2400
-300
-15
-45
Oattr MDg
11000
100
290 30 13 308
28 DM 13 306
26 38 150 340
20 20 79 D20
500
200
1508
3000
200
1988
050
1808
-100
400
-120
-130
-12)
34
358 18 18 208
26 150 240
251
-210
-211
-340
-00
-211
-330
-012
-212
330
30
138 303
-230
2.0
138
310
-211
23 100
750
340
-25
Figure 11: Median envelope for models when using a diagonal Gaussian as approximating distribution using the complete-U DReG, permuted DReG and the standard DReG gradient estimators.
27

---

Published in Transactions on Machine Learning Research (02/2023)
-1000
201
100
30 13 100
70 130
-34
-30
-242
Figure 12: Median envelope for models when using a Gaussian distribution with full-rank covariance as approximating distribution for the estimators complete-U, U, permuted and the standard IW.
28

---

Published in Transactions on Machine Learning Research (02/2023)
2003
100
150
400p 10
-120
+22
Figure 13: Median envelope for models when using a Gaussian distribution with full-rank covariance as approximating distribution using the complete-U DReG, permuted DRcG and the standard DRcG gradient estimators.
29

---

Published in Transactions on Machine Learning Research (02/2023)
Table 8: Median objective averaged over the last 200 iterations when using the estimators complete- U, permuted and the standard IW. It can be seen that for at least 10 models out of 17, using either the di- agonal Gaussian the full rank covariance Gaussian approximation, the best objective is achieved with intermediate value of m, and it is at least nat larger than the objective with m 16. These mod- els are: congress, election88, election88Exp, electric, electric-one-pred hepatitis, hiv-chr, art-multieevel, mushrooms and radon. Optimizations using m = 1 are not shown.
Diagonal Gaussian
Full Rank Covariance Caussian
model
method
m
als
standard IW
-652.8
-649.77
-648.0
-617.6
-639.0
-660.7
-738.1
-772.1
complete-U
-65266
-648.9
-646.7
-647.6
-637.8
-639
-663.9
-772.1
permuted
652.6
-649.0
-647.0
-647.7
-63.88
-639.2
-664.1
-771.2
australian
standard IW
264.4
-261.2
259.11
-288.0
256.8
256.88
-256.8
-257.2
complete-U
-27.2
permuted
264.7
-261.4
259.0
-27.99
-56.88
-566.8
-256.8
-257.1
congress
standard
417.1
419.5
419.5
417.0
416.9
417.2
412.4
403.9
complete-U
118.8
420.3
120.3
117.6
419.4
120.3
119.2
403.9
permuted
418.5
420.2
420.4
417.1
419.4
420.2
418.9
402.9
election88
standard IW
1529.5
-1523.3
-525.0
1535..6
NaN
59644.
-4383.2
NaN
complete-U
1529.77
1521.5
1519.88
1535.6
N&N
4943.1
2046.2
NaN
permuted
5443.7
3300099
NaN
election88Exp
standard TW
-175..7
1570.8
1502.7
1461.2
NaN
NAN
NaN
NaN
complete-U
-1700.0
11496.5
1467.0
-161.2
N&N
NAN
-3748.2
N&N
permuted
17666.6
-1512.3
1482.7
116199
N&N
NGN
NaN
N&N
electric
standard IW
830.5
-27.77
826.1
$30.9
N&N
1421.1
11668.2
1207.2
complete-U
-30.6
-827.0
828.0
-30.9
NAN 1459.2
1090.3
1207.2
permuted
-80..6
-827.1
-2335
$30.9
N&N
14136
1098.1
1203.4
eleetric-one-pred
standard
1153.0
complete-U
1148.3
1145.2
1146.8
1144.6
1150.7
11444
1140.3
1141.2
permuted
-1148.3
-114.66
1146.5
-1111
-115.0
11115
1140.0
1111.2
hepatitis
staudard IW
-56133
-774.9
773.4
-774.4
N&N
NaN
NaN
1693.T
complete-1U
-56.22
-564.3
773.1
-774.4
N&N
1664.4
1592.8
1693.7
permuted
-561.2
-773.8
-75522
-774.4
N&N
1779.6
1716.7
1682.1
biv-chr
standard IW
-606.1
-60144
607.5
-603.8
N&N
NaN
1879.99
1945.0
complete-1/
-00.88
3395.8
1945.0
permuted
-006.2
-601.0
-003.1
-603.7
N&N
NaN
1486,9
1960.9
ionosphere
staudard
-1331
-12.22
-127.2
125.6
-125.3
126.6
132.5
142.3
complete-U
133.2
-29.66
-127.3
-125.6
-124.9
125.2
127.2
142.3
permuted
-13333
-129.6
-177.3
-125.7
-12499
-25.2
127.3
142.2
irt.
standard IW
15887.5
11887.1
11886.8
15886.7
-36563
-64934
68147
N&N
complete-U
118877
-15877.0
11888..8
15886.7
-33230
37383
50316
N&N
permuted
15887.44
15887.0
15886.8
15886.6
-36200
37547
54763
irt-multilevel
standard IW
15201.77
-15191.1
-15111.3
15196.8
N&N
NaN
NaN
N&N
complete-U
15198.7
-15164.0
15185.8
15196.8
NaN
NaN
NaN
NaN
permuted
15200.3
-15173.0
-15186.2
15197.0
NaN
NaN
NaN
N&N
mesquite
standard IW
-9.9
-2..7
-9.66
-9.3
-9.8
-9.77
29.6
-9.22
complete-U
-2.88
-2..6
29.8
permuted
-9.9
-2..7
-9.6
-9.2
-9.8
-9.77
29.6
-29.2
standard
339.0
complete-U
-210.6
201.3
200.8
-215.5
180..2
-18.77
185.6
339.0
permuted
-20088
204.5
201.4
215.44
180.2
-18.77
117.66
337.3
radon
standard IW
-1210.5
-1210.4
-213..3
-2100.4
NaN
2422.9
1595..8
1636.5
complete-U
-1210.5
-2200.2
-1210.2
-2100.4
N&N
-1448.0
1445.9
1636.5
permuted
N&N
sonar
staudard IW
-136.2
-126.3
1213
-117.9
-138.0
-1549
200.99
-25.77
complete-U
-117.9
-120.9
-156.3
-266.7
permuted
-136.5
-127.3
-121.5
-117.9
-116.7
-121.5
-158.2
-228.5
wells
standard IW
-2022.1
204119
2041.7
-2041.2
-2041.9
2041.8
2041.7
-2041.1
completa-U
-2042.2
-2422.0
2041.77
-2041.2
-2041.9
2041.9
-201188
-2041.1
permuted
-2012.2
-011.9
-201.77
-2011.2
-20119
-21188
-00118
-2011.1
30

---

Published in Transactions on Machine Learning Research (02/2023)
Table 9: Median objective averaged over the last 200 iterations when using the complete-U U DReG, permuted DReG and the standard DReG gradient estimators. It can be seen that for at least 8 models out of 17, using either the diagonal Gaussian the full rank covariance Gaussian approximation, the best objective is achieved with an intermediate value of m, and it is at least nat larger than the objective with m 16. These models are: congress, election88, election88Exp. electric, electric-one-prre irt-multileved mushrooms and radon. Optimizations using m are not shown.
Diagonal Gaussian
Full Rank Covariance Gaussian
model
method
als
DReG
-652.7
-649.9
-648.0
-617.0
-659.7
-170.3
-936.6
-2205.1
comp.-DReG
-65255
-648.6
-646.5
-647.0
-635.7
-667
-894.2
-1205.
perm.-DReG
-652.5
-648.7
-646.4
-647.1
-655.3
-725.2
-874.7
12209.8
australian
DReG
-66.44
261.22
-259.0
-57.88
-256.7
-26.77
-256.7
-256.9
comp.DRReG
-26.99
perm--DReG
-664.7
261.55
258.99
257.8
-26.77
-266.7
-256.6
-256.9
congress
417.9
419.7
419.9
418.2
418.5
418.9
comp.-DReG
119.6
420.5
120.7
118.2
420.5
120.8
119.8
401.5
perm.-DiteG
419.4
420.5
420.7
417.9
420.4
420.7
419.8
404.5
elertion88
DReG
-1529.2
-1222.3
1524.22
1534.9
NaN
5964.4
-434922
NaN
comp.-DReG
-1529.1
1520077
11518.4
1534.99
N&N
4950.7
2079.0
NaN
perm.-DRee
NaN 5439.8
3041.2
NaN
election88Exp
DReG
-155..8
15711.9
1502.0
1461.9
NaN
NAN
NaN
NaN
comp-DReeG
-733.22
-1495.9
1168.8
1461.9
N&N
NAN
366400
N&N
perm.-DReG
-1766.8
-1512.3
1483.3
11660.
N&N
NGN
3947.5
N&N
electric
DReG
-830.0
.827.2
-24.8
826.2
N&N
1417.4
4
1291.3
1314.0
comp.-DReG
-830.2
-826.1
822.0
-26.2
NAN
1459.2
1219.0
1314.0
perm.-DReG
-829.9
-826.3
-822.6
-22.3
N&N
1427.22
1239.1
1326.1
elottric-one-pred
1153.0
1145.8
comp--DReG
-1148.5
1146.0
1146.8
1144.6
1150.7
1144.1
1140.3
1141.2
perm--DReG
-1118.3
-1116.7
1146.5
1141.
-115.0
S1TT.
1140.00
-1111.
hepatitis
DReG
-561.3
-77633
775.1
-774.0
N&N
NGN
NaN
NaN
comp.-DReG
-561.1
772.0
-77.6
-774.0
N&N
NAN
NaN
NaN
perm.-DReG
-561.3
-773.5
-774.8
-774.0
N&N
NAN
NaN
N&N
biv-chr
DReG
-606.2
-601.1
605.4
-602.7
N&N
N&N
NaN
N&N
comp--DReG
perm.-DReG
-606.
-603.6
-002.9
-602.7
N&N
NaN
NaN
N&N
ionosphere
DReG
-133.1
129.33
-127.1
125.6
1243
125.8
130.7
142.1
comp.-DRee
133.22
1129.6
-17.22
-125.6
124.2
124.2
124.7
142.1
perm.-DReG
-133.3
-129.6
127.22
-125.6
124.2
-24.3
125.7
142.0
irt.
DReG
-1887.33
15886.9
16886.6
15886.3
N&N
NAN
NaN
N&N
comp.-DReG
-5887.3
-15886.9
11886.5
15886.3
N&N
NAN
NaN
N&N
perm--DReG
15887.3
15886.9
-1886..5
15886.3
NAN
irt-multilevel
DReG
15226.
-15199.1
-15195.8
1122114
N&N
NaN
N&N
N&N
comp.-DReG
15226699
-15188.6
11518800
1522444
N&N
NAN
NaN
NaN
perm.-DReG
15214.00
-15191.6
-15188.4
15222.2
NaN
NAN
NaN
N&N
mesquite
DReG
-2.99
-2.88
-29.6
-9944
29.8
29.7
-9.77
-9.3
comp-DReG
-2.99
-2.88
29.8
29.7
-29.3
perm.DDReG
-29.9
-9.88
-29.6
-2.44
29.8
29.7
29.77
-29.3
DReG
204.7
405.3
comp.-DReG
-210.6
201.3
201.1
-21.99
-180.3
-19..6
253.5
405.3
perm.-DRteG
-210.7
204.4
201.5
215.44
180.4
194.3
-20.99
400.8
radon
DReG
-1210.5
-1210.3
1219.99
1210.3
NaN
24100.8
1624.99
1650.9
comp-DRReG
-1210.4
-1210.1
1210.2
-1210.3
N&N
-1588.0
-145--7
1650.99
perm.-DReG
-1210.4
N&N
sonar
DReG
-136.2
-126.2
-12111
-11..6
-135.4
-152.5
226.3
-29..6
comp.-DReG
-155.6
-259.6
Perm.-DReG
-136.5
-127.3
-121.3
-117.6
-115.3
-118.9
-156.4
-260.5
wells
DReG
2042.2
.2041.9
2041.88
-2041.2
-2041.9
-041.9
-041188
-2011.2
comp.-DReG
-202222
-2042.0
-2041.9
-2041.2
-2041.9
2041.99
-2041.9
-2041.2
Perm--DReG
-2012.2
-2012.0
-201.88
-20112
-20119
-20199
-2011.9
-2011.2
31


Reviews:
[{"review_id": "7dyw0HCihn", "reviewer": "Reviewer_j1dR", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n* The paper is clear and well written.\n* The contributions are novel, interesting, and useful.\n* I checked the proofs in the main text and they are correct to my knowledge.\n* The experimental evaluation is thorough and provides clear evidence to support the proposed methods.\n* Overall it is a great paper, congratulations!\n\nThe main weakness I see is a lack of discussion of 'outer averaging' vs batch averaging for variance reduction. To my understanding, the standard IWAE setup does not do any outer averaging as defined in this paper. Instead, all averaging occurs at the 'batch level', i.e. IWAE bounds for different x's are averaged, not multiple IWAE bounds for the same x. I understand that your section 2.1 is largely dedicated to arguing that practitioners should consider averaging multiple bounds for the same x, but you ignore the option of instead averaging IWAE bounds for different $x$'s. This is an important and practical baseline to consider.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper applies the theory of U-statistics to multisample bound and gradient estimators for variational inference. In many variational inference settings the user allots a budget of samples that can be drawn from the model at each gradient step. Those samples are then combined to estimate the bound and its gradient. The most popular existing method for combining the samples is, roughly speaking, to average each sample's log importance weight (IWAE). This paper suggests using the theory of U-statistics to more intelligently combine the samples and achieve a variance reduction. This can take the form of averaging multiple possibly-overlapping subsets of the weights.\n\nIn addition to these conceptual innovations, the authors propose several estimators based on U-statistics and analyze their performance theoretically and empirically, finding that they match or exceed existing estimators while incurring a small performance penalty.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "I recommend acceptance of this paper without needing any changes. However, there are a few things that would be nice to see. Mainly, it would be nice to add a discussion of batching vs averaging multiple IWAE bounds for the same x.\n\nIn addition here are a few small optional edits:\n\n* On the top of page 3 you say “we use the subscript n to denote the total number of input samples used for estimation and m for the optimized IW-ELBO objective”. Although I was able to figure it out, it's not really clear to me from this sentence what m is being used for.\n* On the bottom of page 4 you say \"can be generalized by replacing h by any other symmetric function of m variables”. While 'symmetric function of m variables' does have the technical definition 'invariant to permutations of the variables', I had to look it up just to make sure. It could be better to be more precise here, but it is probably fine as-is.\n* In Definition 3.1 you end saying “which depends only on c and not the particular s and s'” It kind of sounds like you are making a claim, not just a definition when you add the above statement. Maybe move it out of the definition?\n* It could be good to show table 2 for different values of n and m to support your claims in the text. The additional info in the appendix for the KMNIST experiments kind of gets at this this, but it could be good to support your claims more directly with experiments on the mushroom dataset.\n", "broader_impact_concerns": "No broader impact concerns.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "wWbMF2kRwR", "reviewer": "Reviewer_JAwY", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths\n\nThe idea of using U-statistics in this context is quite natural and excellent. Overall, the new objectives are practical, and often more accurate than standard IWAE-style bounds. We can hope that this work can benefit to several researchers using such bounds.\n\nThe paper is generally well-narrated, and cites most of the relevant literature.\n\nThe authors are refreshingly honest about the fact that the improvements are often slim.\n\nWeaknesses\n\nAs mentioned above, the empirical gains are often slim.\n\nAlthough the paper reads well, a few things are a bit unclear:\n\n1) The maximum function is nonsmooth, which may be problematic for gradient descent? Discussing this would be interesting.\n\n2) Some of the maths deserve a bit more details, for instance:\n    - In the beginning of Section 2, $x$ and $z$ are not properly defined, and we have no idea where they live\n    - In Def 3.1, it would be helpful to briefly say why it only depends on $c$\n    - It would be nice to add a proof of Prop 3.3 (perhaps in the appendix) for completeness. While Prop 3.2 is exactly Hoeffding's result, Prop 3.3 involves some straightforward steps that would still benefit from being spelled out\n    - It is perhaps well-known, but I was not familiar with the notation $\\binom{[n] }{m}$, maybe specifying what it means would be helful\n    - The computation needed to obtain Eqn (8) could be detailed\n    - Right after Eqn (8), you say \"for a given permutation $\\pi$, all sets in $S_\\pi$ will be independent\". I don't understand, given $\\pi$, aren't all sets actually deterministic?\n    - In Appendix B, the paragraph before the theorem is nice but explaining the rationale of the conditions a bit more would be helful. For instance, why is $\\tilde{\\mathcal{F}}$ the class of finitely supported distributions? I am not particularly familiar with Halmos's famous paper, but the link between your result and the way he phrases Theorem 5 is not that obvious.\n\n3) The use of the \"envelope\" of the curve in Figure 2 and later is quite peculiar, but not uninteresting. Is this something common? If yes, maybe you could motivate it a bit more (it is quickly explained in page 11, while Fig 2 is page 8, which is not great for the flow), and add citations of other papers that use it?\n    \nSome statements are incomplete or not fully correct:\n\na) Page 1: \"The most obvious downside is the increased computational cost (by a factor of $m$)\". Depending on the context, this can be smaller than this. For IWAEs, computing the standard bound requires a single encoding, and $m$ decodings, whereas increasing the complexity by a factor m would mean $m$ encodings, and $m$ decodings.\n\nb) Page 2: when you mention that the IWAE bound is consistant, you should mention that this is under some assumptions on the distribution of the importance weights (Burda et al. assume boundedness, Domke and Sheldon have weaker conditions)\n\nc) Page 2:  \"which simultaneously finds an approximating distribution that is close in KL divergence to p(z | x) (Domke & Sheldon, 2018)\" I believe Domke & Sheldon actually show that, when $K \\longrightarrow \\infty$, the distributions are close in $\\chi$ divergence instead of KL. I am not aware of precise results on what's happening for $K \\in \\]1,\\infty[$.\n\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The authors propose a new general recipe to derive variational bounds, in the spirit of importance weighted variational inference (IWVI), introduced by Burda et al. (2016) in the Importance-Weighted Autoencoders (IWAEs) paper. \n\nThe idea is quite simple yet novel in this context: using U-statistics instead of standard averages. Concretely, this means that the variational bound will be evaluated many times on different overlapping batches of importance weights. The novelty being the fact that the batches can be overlapping, and are much more numerous than for standard IWVI bounds. The authors use the standard theory of U-statistics to show that these new bounds have a lower variance than standard IWVI bounds. The downside of their new bounds is that computing all possible combinations of batches can be computationally intensive. The authors propose several solutions to this problem: sampling fewer batches than the complete U-statistic would, or relying on the fact that the logsumexp function can be approximated by the max function.\n\nIn several practical situations (Bayesian models and IWAEs), these new objectives are slightly more accurate than standard IWVI.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "- Clarifying/correcting the few points mentioned in the \"weaknesses\" section\n\n- There are a few points/references that would be interesting to discuss:\n    - You mention \"we are not aware of control variates specifically for IWVI\". Control variates for IWVI were developed by Mnih and Rezende (ICML 2016, https://arxiv.org/abs/1602.06725) and revisited by Liévin et al. (NeurIPS 2020, https://arxiv.org/abs/2008.01998).\n    - DReG was extended by Bauer and Mnih (ICML 2020, https://arxiv.org/abs/2101.11046)\n    - It would be nice to discuss whether or not your ideas could be applied also to the reweighted wake sleep algorithm of Bornschein and Bengio (ICLR 2015, https://arxiv.org/abs/1406.2751) and its extensions (eg Dieng and Paisley, https://arxiv.org/abs/1906.05850, Le et al., https://proceedings.mlr.press/v115/le20a.html, Kim et al., https://proceedings.mlr.press/v108/kim20b.html)\n\n\n", "broader_impact_concerns": "I have no particular concern about this work.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "2mLipLIucE", "reviewer": "Reviewer_7r3U", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "## Strengths\n- The paper is very well-written and easy to follow. Related works are well-cited and enough technical background is provided.\n- The proposed method is simple but effective. The use of U-statistics is very elegant and could motivate similar application in other related estimators for variance reduction.\n- The paper consists a nice mix of theoretical results and empirical evaluation to study the proposed estimators.\n    - Important theoretical results are established for guaranteed improvement from the proposed estimators as well as some guides of hyper-parameter selection.\n    - Empirical results are aligned with the theory and show practical improvement by using the proposed estimators.\n\n## Weaknesses\n- Only a few typos or presentation issues\n    - Missing brackets ($[$, $]$) in a few expectations ($\\mathbb{E}$) within the paragraphs between (5) and (6)\n    - The statement in Proposition 3.2 is incomplete. The first sentence should end up with something like \"... for some $\\zeta_1, \\zeta_m > 0$\" (otherwise they are not defined).\n\n## Questions\n- I wonder how the complete U-statistics estimator is implemented in practice and if any vectorization is used to accelerate it instead of a loop-based implementation. \n    - For example, if we denote the weights as $\\mathbf{e}=[e_1, e_2, e_3, e_4]$ and let $m=2$, one could write the complete U-statistics estimator as $h(\\mathbf{e}^\\top \\mathbf{M}) \\mathbf{1}  / r$ where $\\mathbf{M}$ is a $4 \\times 6$ binary matrix whose columns are indicators for each subset (e.g. the first column is $[1; 1; 0; 0]$, the second column is $[1; 0; 1; 0]$, etc.) and $h(x) = \\ln {1 \\over m} x$ is applied column-wise. This might be effectively parallelized if the vector-matrix multiplication is done on GPUs. \n    - Just throwing a random idea: Maybe this vector-matrix multiplication could be somehow approximated efficiently (say by some low-rank approximation or even random projection).", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper improves the bias-variance trade-off of IWVI by using U-statistics. The main idea is to carefully construct overlapped batches to compute individual estimate for average. Several non-trivial estimates based on this idea are introduced, including\n1. The complete U-statistics IW-ELBO estimator, which is computationally intense but achieves optimal variance reduction;\n2. The permutation-based incomplete U-statistics IW-ELBO estimator, which allows controlling the computation by a hyper-parameter $l$ and achieves variance reduction in between the optimal and an incomplete estimator based on random selection;\n3. Two approximate complete U-statistics IW-ELBO estimator, which uses the log-sum-exp approximation to allow fast approximation of the complete U-statistics.\n\nNotably, theoretical results for each of these estimators are obtained in terms of variance. Most of them have already shown the advantage (better bounds) by using the newly proposed estimators and some of them also gives practical guidance on selection of hyper-parameters such as $l$ based on Proposition 4.2.\n\nA large suite of simulation for IWVI is performed to validate the practical improvement by using the proposed estimators.\nBenchmark on the actual computation is also provided to show the disadvantage of the complete U-statistics IW-ELBO estimator, and the importance of reducing computation achieved in the other estimators.\nFinally the proposed estimators are tested in IWAEs, resulting in improved lower bound of the marginal likelihood.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "I think the paper could be accepted as it.", "broader_impact_concerns": "No concern as this work is mostly theoretical", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: All reviewers were very supportive of the paper and recommended its acceptance without revision.  I concur with this assessment and believe it to be a clearly presented paper with interesting and notable novel contributions to the literature.  
Decision: Accept (poster)