Forward $\chi^2$ Divergence Based Variational Importance Sampling
Authors: Chengrui Li
Yule Wang
Weihan Li
Anqi Wu
Venue: ICLR 2024 spotlight (2024)
Keywords: Importance sampling
$\chi^2$ divergence
latent variable models
Primary area: probabilistic methods (Bayesian methods, variational inference, sampling, UQ, etc.)

Abstract: Maximizing the marginal log-likelihood is a crucial aspect of learning latent variable models, and variational inference (VI) stands as the commonly adopted method. However, VI can encounter challenges in achieving a high marginal log-likelihood when dealing with complicated posterior distributions. In response to this limitation, we introduce a novel variational importance sampling (VIS) approach that directly estimates and maximizes the marginal log-likelihood. VIS leverages the optimal proposal distribution, achieved by minimizing the forward $\chi^2$ divergence, to enhance marginal log-likelihood estimation. We apply VIS to various popular latent variable models, including mixture models, variational auto-encoders, and partially observable generalized linear models. Results demonstrate that our approach consistently outperforms state-of-the-art baselines, in terms of both log-likelihood and model parameter estimation. Code: \url{https://github.com/JerrySoybean/vis}.

Published as a conference paper at ICLR 2024
FORWARD x² DIVERGENCE BASED VARIATIONAL IM- PORTANCE SAMPLING
Chengrui Li, Yule Wang, Weihan Li & Anqi Wu School of Computational Science & Engineering Georgia Institute of Technology, Atlanta, GA 30305, USA {cnli Sonliched melichhegrri, yulewang, weihanli. i, angi www.weeegggccheccooooom edu
ABSTRACT
Maximizing the marginal log- log-likelihood is crucial aspect of learning latent vari- able models, and variational inference (VI) stands as the commonly adopted method. However, VI can encounter challenges in achieving high marginal log- likelihood when dealing with complicated posterior distributions. In response to this limitation, we introduce a novel variational importance sampling (VIS) ap- proach that directly estimates and maximizes the marginal log-likelihood. VIS leverages the optimal proposal distribution, achieved by minimizing the for- ward x2 divergence, to enhance marginal log-likelihooe estimation. We apply VIS various popular latent variable models, including mixture models, vari- ational auto-encoders, and partially observable generalized linear models. Re- sults demonstrate that our approach consistently outperforms state-of-the-art base- lines, in terms of both likeliield and model parameter estimation. Code: aattp::/gittini com wtw.uuu....o.eerryyyyyaaaaa.aaaa
INTRODUCTION
Given the latent variables and the observed variables how to find the optimal parameter set 0 that produces the maximum marginal likelihood p(æ; 0) = p(x, 2; dz is essential in a wide range of downstream applications. However, when the problem is complicated. we only know the explicit form of p(æ, z; 0) and it is intractable to compute the marginal p( æ 0) analytically. Therefore, we turn to approximation methods such as variational inference (VI) (Blei al. 2017) and importance sampling (IS) (Kloek & Van Dijk 1978) to learn the model parameter B and infer the intractable posterior p(z ().
VI uses variational distribution q( æz; () to approximate the posterior p(z| 0) with the difference as their reverse KI. divergence (T.(4(2);; Φ) [|]pp]/(] f)). where minimizing the KI. divergence equal to maximizing the evidence lower bound (LBO 0, Φ) of In p(æ; 0). However, maximizing In p(2; ()) using ELBO may not be a good choice when dealing with complex posterior distributions, such heavy-tailed or multi-modal distributions. There is chance that (L(g(=) (æ; d) 2 æx; 0)) is very small, but in fact both q(2 x Φ) and p(z|æ: x;; 0) are far from the true posterior p( 2 [z; grree, leading to higher ELBO but a lower marginal log- oigliiilloood (e.g., Section 4.1).
Although other bounds such as a divergence- lower bound (Li & Turner 2016 Hernandez- Lobato et al. 2016) and x² divergence-basee upper bound (Dieng al. 2017) can be used for better posterior approximation, a more straightforward approach is to estimate In po (z;; 0) by IS. Ideally, IS could have good estimation if choosing a proper proposal distribution q( ((2(23 ac; () and a large number of Monte Carlo samples. In practice, however, there is often a lack of clear guidance on how to choose (z [æ; σ) and limited indicators verify the quality of q(2) [2; Φ). SS & Chen (2021) showed that the variational distribution found by VI could serve as proposal distribution for IS, but is not the optimal choice (Jerfel et al. 2021; Saraswat 2014 Sason & Verdū [016 Nishiyama & Sason, 2020). Besides, Pradier et al. (2019) noticed the numerical and scalability issue in minimizing forward x divergence Finke & Thiery (2019) which should be treated rigorously.
To address these issues, we propose novel learning method named variational importance sampling (VIS). We demonstrate that an optimal proposal distribution g(> æx; σ) for can be achieved by minimizing the forward x" divergence in log space, which is numerically stable. Furthermore, with enough Monte Carlo samples, the estimated marginal likelihoodo In p(x:; 8) is an asymptotically tighter lower bound than ELBO, and hence In p(æ; 8) could be maximized more effectively. In the experiment section, we apply VIS to several models including the most general case when there is explicit decomposition pl a. 2; 0) p(x 2; 0)p(22 (2; 0), with both synthetic and real-world datasets

---

Published as a conference paper at ICLR 2024
(a)
(b)
(c)
p(z)a; 6)
lp p By (1
x (pllo))((2)
minimized by KL(glp))
minimized by xx(llg)
ELRO in6
(L(llpp)
((=(=;)))) ((=lo))))))
inp
ELBO. E. [ELBO]
family of o(z)2;))
Figure (a): The bias between the marginal oigeliilddood In pn; 8) and the expectation of its IS estimator Fy [In p(æ;; 6)]. the ELBO(#; 0, 6), and the expectation of the ELBO's estimator E, (ELLBO 0, +)]. When estimating In p(a; 0), the down-biased IS estimator
0, +) is tighter lower bound than the down-biased BLOO estimator [FLBOO (x: 0, ()). (b): Empirical visualization of the four quantities in (a) with different Monte Carlo samples K E {1,2, 3, 4, 5) Each box in (b) is based 500 repeats and the hollow circle the box is their average. An asymp- totic difference occurs when increasing K. (c): Different q x ②) are obtained by minimizing the forward x² divergence, which is optimal for doing IS v.s. by minimizing the reverse KL divergence.
to demonstrate its superiority CHIVI (Dieng et al. 2017), VBIS (Su & Chen (021), and IWAE (Burda et al. 2015). Appendix A.8 summarizes the related works and our corresponding contributions in table.
the most widely used VI and three other state-of ofettheart methods:
2 BACKGROUND OF VARIATIONAL INFERENCE
Here give brief introduction the variational inference (VI), its empirical estimator, and its bias. VI starts from the reverse KI. divergence:
(T () 0)) / ((2|00; )) In ( dz ELBO(a;; 0, + In p(2;; 0), (1) p(z
-E,, with (LBOO 0) = [n p( 0) q -) the problem maximizing In p( 0) is converted to maximizing ELBO(#; Φ). VI is often fa- vored for several reasons, such likelihood, making it numerically more stable compared to working directly with the original like- lihood; 2) when the model can be factored p(z, 2; 2) p( f)p(p)(()))) 2: 6), the ELBO can be reformulated as (LBO 0, Φ) = En p(x) z; 0)] - <L<((=/<;; Φ) [pp(;; 0)). This formulation is advantageous because the second KL term often has closed-form expression for specific choices of the prior distribution p(2; 0) and the variational distribution family {(2|ar; Φ), such as the Gaussian distribution.
az; x)). Since ELBO is lower bound of In p(ar; 0),
The ELBO is formulated in terms of expectations of log-
In practice, the target function ELBO in Eq. still requires numerical estimation, resulting in an empirical estimator
K
Inp ar, =()) In (2(+) 22; (
(2) k=1 are K Monte Carlo samples from the variational distribution q(2(a;;) Now, we convert maximizing In p(x; 0) w.r.t. maximizing FI,BO 0, w.r.t. and Φ. The score function and pathwise gradient estimator of ELBO are shown in Appendix A.I
FLBO
K
where (2(1) K
k=1
Bias of the ELBO estimator. Note that although ELBO x Φ) is an unbiased estimator of ELBO, it is strictly down-biased estimator of the marginal log likellooddod In pp; (Fig. [a)):
E, FI:BO(=;; 0) p(x; 0) FI:BO Φ) p(æ; 0) = (T((2()) ()) m. 0)).
(3)
As mentioned before, there is a chance that both <<>|<; Φ) and p(> x;; 0) are far from the true posterior p( [z; ginie), resulting in higher ELBO but lower marginal log-likelihle likelihood In p(: 0).
3 VARIATIONAL IMPORTANCE SAMPLING
To tackle this problem. we use importance sampling (IS) to estimate the marginal log likelihodd In p( 0) directly. However, the approximation quality of IS depends the choice of the pro-
2

---

Published as a conference paper at ICLR 2024
posal distribution and the number of Monte Carlo samples. We first show that using IS can get an asymptotically tighter estimator of In p(æ; 0) than (=: 0 Φ). Then, we prove that the bias and effectiveness (variance) of this estimator are both related to the forward x² divergence and the num- ber of Monte Carlo samples. This provides guidance on how to select the proposal distribution and the number of Monte Carlo samples. Finally, we derive the numerically stable gradient estimator used for obtaining the optimal proposal distribution.
Down-biased IS estimator of the marginal log-likelihood. With importance sampling (IS), the marginal is approximated via a proposal distribution g( ax; Φ). i.e., (x, (k)
p(a; 0) = pla: 2;0)
(z(k) x; Φ) = p( 0))
(4)
K
where {=(4)] J k= are K Monte Carlo samples from the proposal distribution q(= 22; 6). For numer- ical stability, we need work with it in log space,
In p(æ: 0, logsumexp In (x. (8); In (2(h) x=; +)]
In K,
(5)
where the logsumexp trick can be utilized. Appendix A.2 shows that the gradient of p(x; 0) w.r.t 0 can be estimated as
Ø Øln 0) Ø In p(a:; 0,))
(6)
00
00
Since
K
p(a,, 2.0
Eg [eppp 0,
x0 = K k=1
Σ
Eg
plxx, 2; dz = p(x; 0),
(7)
[9(2|=;))
p(x;0, 0) is an unbiased estimator of p(æ; 0). However, In(.) is a concave function, thus E, [n p(æ; 0, 6) < In [8(;;;;;; Φ)] = pn; from Jensen's inequality. This means the esti-
mator in log space In ((:;;) () is down- biased estimator of In ());
Bias of the IS estimator. Similar to ELBO x;; 0, d). the Delta method (Oehlert 1992 Struski et al. 2022).
we can derive the bias of In p( 0, o) with
p (2(1)) ax; 0)
ạ (z(A) 20; ()) )
p( O Φ) In p(:; 0) = E.g
In
K 1=
Varg g [02222 Φ)
p(z|<x;0) 0)
1
Eyg
p==|a>; 0)
R
[p(| ax; 0) g(2 Φ)
(8)
2K
2K
1
p(z 0)
1
dz 1)
2K
<(2|<;)) 6)
2K
where x" (plll) the forward x2 divergence between p and q (Fig. [)). Since Eq. 8 converges to as K + c0. In p( 0, σ) is an asymptotically tighter lower bound than ELBÖ(x; 8, Φ) (Fig. a)). Particularly when K 1, In pnx; 0, Φ) ELBO( A, Φ). To verify this relationship empirically, we repeat the estimation of in p(æ; 0, Φ) and (LBO( 6, Φ) based on K Monte Carlo samples 500 times, and plot their empirical distributions w.r.t. K in Fig. Tb). With more Monte Carlo samples K, both In and ELBÒ become stable, but the empirical expectation indicated by the hollow circle in each box implies that only pn; 0, Φ) converges to the marginal log- likelihood In p(æ; 8).
Fig. demonstrates that IS can have a much better pn; 0) estimation by setting a large K, which means using IS is more direct way to maximize In p(;; 0) than ELBO. Besides, to have a faster convergence, we also need to choose the proposal distribution q<>| Φ) that minimizes x (p(2|=x [l 2 æx; σ)) since this forward x² divergence could serve as an indicator of whether the proposal distribution is good: if the forward x² divergence is small, then the bias (the absolute value of Eq. & of the IS estimator is small.
On the other hand, we can write down the effectiveness p(:; 0, Φ). i.c.
(Fredman et al. 1998 of the estimator
1 K2
x( x: p)) 0) ((z [z|=x:00 Φ)
p(x; 0)2
Var, (1(;0.,)
K Var,
K x² (p(2) < 01] (9)
3

---

Published as a conference paper at ICLR 2024
Algorithm VIS
for I:N do
2: Sample {=(*) k=1 from g(=|æ;;) Φ).
3: Update by maximizing In (a;;; Φ) via Eq.
4: Update o by minimizing <=((= ( o)) via Eq. 12 or Eq 24
5:
end for
which is the variance of the estimator. Eq. and Eq. coincide to indicate that for a small bias of
In p(æ; 0,, () and a high effectiveness of p( 0, Φ). we want a small x (p<x=aæ;; (x2));))]]g( 0)]lg(z æ 6)) and a large K. In other words, we need as many Monte Carlo samples as possible; and the optimal choice of the proposal distribution for IS is the (=]|2x; Φ) with the minimum forward x" divergence x2 (p( 2/2;; ()((())()))))))))))))) x;; Φ)) rather than reverse KL divergence (L(((z)æ;
Φ) ()lpp((( 0)) (Fig. (()).
The algorithm of the variational importance sampling (VIS) is summarized in Alg. I We first perform IS to maximize In p(2; 0, 6) w.r.. 0, given a fixed proposal distribution g x Ø); then we fix 0 and minimize x²(((z a:; (=][= 6)) w.r.t. & to obtain a better proposal distribution for doing IS. However, minimizing x (p(22 (0p(((a;))) ar; 0) (=|=; ()) w.r.t. @ is non-trivial since we don't know p(z| z;; 0). We derive a stable gradient estimator for minimizing the forward x² divergence in the following.
Gradient estimator.
Rewrite the forward
x² divergence
1
/
p(æ, p( 2; 0)
(v(=aæ;; (1((2=
p(x;; 0)2 (2((2;;) σ) dz - 1
p(ar; 0)2 0xx;0.0) 1. (10)
So, minimizing x" (p([<x (=]/(=/2; 6)) is equivalent to minimizing V(æ; 0, Φ) = p(m.=:0)'
q(z x:인ه dz
w.r.t. Φ. It still needs to be estimated and minimized in log space for numerical stability (Pradier al. 2019 Finke & Thiery
Geffner & Domke, 2020 Yao al. 2018). In Appendix A.3
derive that in V(æ; 0, Φ) can be estimated as
In (zx; 0, o) logsumexp 2 In
(c, (k) 0) 2 Ing (x()) x;; o)
- In K = In Vo 0, 6).
(11)
The score function gradient estimator of In Vax; 0, Φ) w.r.t. @ at is
Oln (a: Φ) á In (2:0,
(12)
0o
do
When the reparameterization trick can be utilized, z| a Φ) where r(e), then have the transformation <(>/2;; Φ) dz = r(E) de (Schulman et al. 2015). Now, we can get the à In Va; 0, Φ), where we sample e ~ r(c) and use
pathwise gradient estimator 0ln (a)3..) de (k) g (e!*) Φ) in In (o; d). The derivations are shown in Appendix A.3
4 EXPERIMENTS
Baselines for comparison.
We will apply VIS on three different models and compare it with four
alternative methods:
VI: The most widely used variational inference that maximizes ELBO.
CHIVI (Dieng al. 20177: When updating o, use both an upper bound CUBO (based on forward x2 divergence) and a lower bound ELBO (based on reverse KL divergence) to squeeze the approximated posterior q( æx; Φ) the posterior p( 2/æ=;) 0).
VBIS (Su & Chen 2021 l): Use the g(2|2z Φ) learned from VI as the proposal distribution of IS. IWAE (Burda al.) (015)): The importance-weighted autoencoder. It uses IS rather than VI to learn an autoencoder. An additional competitor for the VAE model only.
Metrics. For all models and datasets, we train the model with different methods on Turin and evaluate on Frest by: marginal log-likelihooe (LL) p(2rest 0), which can be evaluated on both synthetic datasets and real-world datasets; complete log-likelihood (CLL) p( Ziest; 0), which can be only evaluated on synthetic datasets, since we have the zest when generated the data; and hidden likeliood (HLL) ( Drstti Ø), which can be only evaluated on synthetic datasets for the same reason above.
4

---

Published as a conference paper at ICLR 2024
(a) -0.6299
(b)
VI
CHIVV
V815
0.6300
6301
-30
0.65
0.63
0.55
0.65
0.0
0.60
0.60
0.80
0.53
0.50
0.50
0:50
0.50
5500
100
200
130
200
opachs
(c)
0.3
Posterior VI
Posterior CHIVI
Posterior
Posterior
00
-15
-10
pl- Digrrn)
1on 052x 0::1
pld 0:
pl0 2, ooe QU2X 1:
Figure 2: (a): LL, CLL, and HLL evaluated on the test dataset. (b): Convergence curves of the parameter set learned by different methods. The dashed curves are the true parameters used for generating the data, and the solid curves are the learned parameters. (c): The posterior distribution given z and x = 1 in different methods. The dashed curves are the true posterior p(z[2; (trt). the solid curves are the learned posterior p(z|ær; 0), and the dotted curves are the approximated posterior q== 6) learned in the variationa//proposal distribution.
4.1 A TOY MIXTURE MODEL
Model. We first use a toy mixture model to illustrate some representative behaviors of differ- ent methods. Consider the generative model p(z;0) = E T N(2; Hi, 13) with =
1- T3 ,4 J. and p(æ122; 0) Bern(r; logistic())). The parameter set is 0 { 2 the latent variable is E R, and the observed variable x E {0, 1}. Choosing the varia- tional/proposal distribution family tional/proposal parameter set is Φ = (co: C1, oo J.). The simplicity of this model enables the visualization of p(2/2; 0) for us understand the behaviors of different methods.
q(=(2;;)
N/22; Ca. Cx, o2) for
10, and the varia-
Experimental setup. Both the training set and the test set consist of 1,000 samples simulated from the p(x, 2; grrue). We Adam (Kingma & Ba 2014) the optimizer and the learning rate is set at 0.002. We run 200 epochs for each method, and in each epoch, 100 batches of size 10 are used for optimization. The number of Monte Carlo samples used for sampling the latent is K 5000. We repeat 10 times with different random seeds for each method and report the performance.
Results. Quantitatively, VIS performs consistently better than all other methods in terms of all three metrics (Fig. ((a)). In Fig. bb), we plot the convergence curves of the parameter set 0 learned by different methods. Clearly, VIS achieves a more accurate parameter estimation. This further validates that a better parameter estimation corresponds to a higher test marginal log-likel
To understand the effects of the approximated posterior g(a æ Φ) learned by different methods, we plot the true posterior p(z|z; (tTle) (dashed curves), the learned posterior p(z|2; 0) (solid curves), and the approximated posterior <=>[;; Φ) (dotted curves) conditioned on x and Jr respectively in Fig. ( First, we can tell that the true posterior p(z[a;; rz; @tree) conditioned on both x 0 and I are multi-modal shaped, with at least two distinct bumps. For example, p(=|æ 0; (iree) has one large bump centered at about
-8, one large bump centered at about 22, and one small bump centered at about z = 1 (see the purple dashed curve in Fig. c()). Then we check the learned posterior p( 2 0: 0) and the approximated posterior q(
0: o):
behavior of minimizing the reverse KL divergence in VI makes the two large bumps the left collapse into one. And the support of (=>|>: 0; Φ) only
For VI, the zero-forcee wewoo-ooring//ode-eeeeeeee
5

---

Published as a conference paper at ICLR 2024
-140
(b)
O
/
2
8
4
5
6
7
8
9
V1
0
2
8
9
5
6
7
8
9
0
1
8
9
5
6
7
8
9
160
0
1
2
4
5
6
7
8
9
-170
CHI VIS
VIS
0
1
2
4
3
7
8
9
10 epoch
15
IWAE 0
/
2
8
9
3
7
8
9
Figure 3: (a): The marginal log- likelihood on the test set w.r.t. the training epoch. (b): Examples of raw images and the reconstructed images by different methods.
covers the left large bump of p(2 x 0; (the)), which leads to p(2|z = 0; 0) have very different shape to the p(=|2: 0; (rmme.. This is the case that the reverse KI. divergence KL(g(= a:; (p( æx; 0)) is very small, but in fact both {(=| Φ) and p(z æ 8) are far from the true posterior p( zz; ftmee. leading to higher ELBO but lower marginal log-likelihood.
For VBIS, through importance sampling, the learned posterior p(=|æ - 0; 0) maintains two large bumps but the small bump centered about
is still covered by q(=|2 ())) due to the zero-foreee eroo-foring behavior of minimizing the reverse KL divergence. Besides, since the g (z; Φ) learned by minimizing the reverse KL divergence not the optimal proposal distribution for doing IS, the learned p(=|ar; 0) is not good enough to match the true pl 22; (tmee well.
For CHIVI, both the reverse KI. and the forward x2 divergences are considered, so that the support of <<(|x 0; 6) becomes much wider make sure the density under both of the two large bumps can be sampled. However, it is still not wide enough to cover the small bump centered at about 1 compared with VIS. Besides, since CHIVI optimizes ELBO rather than the marginal log- likeliidl w.r.t. 6, the learned 0 not better than VIS.
For VIS, the mass-covvrrng//men--seee makes the q<=|<r = 0; Φ) wide enough to cover both the two large bumps and the small bump centered at about
behavior of minimizing the forward x² divergence
1. Moreover, since have shown in Eq. 8 and Eq. that the q(2]a;; az; Φ) learned by minimizing the forward x2 divergence is the optimal proposal distribution for doing IS, the shape of the learned posterior p(z[æ; 0) matches the shape of the true posterior p(z (inee) the best compared with other methods.
4.2 VARIATIONAL AUTO-ENCODER
Model. The generative model of variational auto-encoder (VAE) (Kingma & Welling 2013) can be expressed as p(2; 0) = N (2; 0, I): and p(z] ;; 0) = Bern(æ; (nggisticMM.Paaar )z))). The parameter set & consists of all parameters of the MLP decoder. The variational/proposall distribution is parameterized as g( p σ) = N(=; (x),, o² (x)
where (x) and a (x) are the output of the MLP encoder given input The parameter Φ consists of all parameters of the MLP encoder.
Experimental setup. We apply the VAE model on the MMIST dataset (LeCun et al. 1998). There are 60,000 samples the training set and 10,000 samples in the test set. Each sample grayscale wrndtwritten digit, E [0. 17784, For visualization, we set z E R². Similar to KKigmaa Welling 2013), we set encoder and decoder structure as
28
MLP dee (ec(z) W Wecc22har b
hdec = tanh Wd 1 + bsec.)),
haee E 12
(ax)) = Wrhbene + by
(13)
hene = tanh (W + bene). hene E R2 128
In (z) Wohene b
We use Adam (Kingma & Ba 2014) the optimizer and the learning rate is set at 0.005. We run 20 epochs for each method. The batch size is set as 64. The number of Monte Carlo samples used for sampling the latent is K = 500. We repeat times with different random seeds for each method and report the test log-likelibhee likelihood.
Fig. 3aa plots the marginal logelikel the test set during learning. As the typical solver, VI performs roughly the same as CHIVI and VBIS, but the convergence curve of VI is a bit more stable. When comparing them with IWAE and VIS, however, IWAE is better and VIS is the best. The reconstruction images shown in Fig. 3)) also imply that VAE solved by VIS could provide good reconstructions similar the corresponding images. The learned latent manifolds by different methods are shown in Appendix A.4
Results.

---

Published as a conference paper at ICLR 2024
4.3 PARTIALLY OBSERVABLE GENERALIZED LINEAR MODELS
Model. We first present the classical generalized linear model (GLM) PPlloo et al. 2008), which studies multi-neuron interaction underlying neural spikes. We denote spike train data as Y e NTXN recorded from N neurons across T time bins, Yt,n as the number of spikes generated by the n-th neuron in the t-th time bin. When provided with Y. classic GLM predicts the firing rates ft. of the n-th neuron at the time bin as
N
1 [ Ut -t--taehh 1=1
))))) (t=ttar()))
St.n = bn
Wne-n'
with spike Vt, Poisson(Sr.n), (14)
n==1
where (-) is a non-linear function (c.g., Softplus); 6n is the background intensity (bias) of the n-th neuron whose vector form is E RN; Wne-n' is the weight of the influence from the -' neuron to the n-th neuron whose matrix form is W RNNN, ab E R4 is the pre-defined basis function summarizing history spikes from t L to 1.
The classic GLM is not a latent variable model. However, we extend GLM to a partially observable GL.M (POGLM) (Pillow & Latham] (007), which becomes a latent variable model. Specifically, POGLM studies neural interaction when the spike data is partially observable, which is often the case in neuroscience since it is usually unrealistic to collect all neurons in a target brain region. Consider group of N neurons where V of them are visible neurons and H of them are hidden neurons (with N = V + H). Given the spike train Y, we denote its left V columns as X YI.T..V E NT*V containing the visible spike train, and the right H columns as Z YIT.V++I:N NT* containing the hidden spike train. Then the firing rate
L
11 =0 bn +
Σ n'=1
Wne-nn
Σ
It In- V/
Σ n'=1+V
Wne-n'
Σ 21-1 (o-----Vi)) (15)
for both visible and hidden neurons. Since the hidden spike train is not observable, POGLM be- comes latent variable model with observed variable Fi,n and hidden variable 2t,n- The model parameter 0 is defined to be {b. W). The graphical model of POGLM is sketched in Fig. 4 a) top.
To do VI, VIS. or others on POGLM, a commonly used variationaalprooosaal distribution (Rezende & Gerstner 2014 Kajino 2021) is { 11:11,1:01) [1aa = Poisson(fe..) where ft,m is de-
fined in Eq. [5] Note that when using Eq. 15 to define the variational/proposall distribution, {b, W)
forms the variationna//roposaal parameter set Φ. The graphical model of the variational/prooosaal is
sketched in Fig. aa) bottom.
4.3.1 SYNTHETIC DATASET
Experimental setup. We randomly generate 10 different parameter sets @ of the GLM models for
data generation, corresponding to 10 trials. There are N neurons in total, where the first V =
neurons are visible and the remaining H = 2 neurons are hidden. For each trial, we simulate 40
spike trains for training and 20 spike trains for testing. The length of each spike train is T 100 time bins. The linear weights and biases of the model used for learning are all initialized as Os. We use Adam (Kingma& Ba 2014) as the optimizer and the learning rate is set at 0.01. We run 20 epochs for each method, and in each epoch, 4 batches of size 10 are used for optimization. The number of Monte Carlo samples used for sampling the hidden spikes is K = 2000. We repeat 10 times with different random seeds for each method and report the performance.
Results. From the barplot in Fig. 4b), we can see that VIS performs significantly better than the other three methods in terms of all three metrics (LL, CLL, and HLL). Similar to the toy mixture model. we can also check the parameter estimation and compare them with the true parameter set used for generating the data. The average weight and bias error are presented in the rightmost two bar plots in Fig. (b). The weight error of the VIS the smallest. For the bias error, both VBIS and VIS are the smallest and are significantly smaller than VI and CHIVI.
In Fig. 4c), we also visualize the parameter recovery results from different methods. For the bias vector, we can visually see that VI and CHIVI are worse than VBIS and VIS. For example, the bias of neuron is positive, but only VIS recovers this positive value. For the visible-to-visible weights (the lep-left block of the weight part), all four methods hidden-to-visile -iiiblle weights (the top-right block of the weight part). VI and CHIVI do not get enough gradient due to maximizing ELBO., these weights are still kept around O. For the visible-to-hidden
match the true well. For the
7

---

Published as a conference paper at ICLR 2024
(a) ptx. 2.9
(b) 23200
-00
0.215
2322
0,270
0.26
6.72
252.
-003
3,000
659
6.38
233.0
0,053
(c) linear
CHIVI
VBIS
VI5
0.2
0.0
-0.2
(d)
predictive firing rates
3.4
tem
Enue
CHVI VB5
1.0
Ttom
0.0
3.5
Figure 4: (a): Graphical model of p(X, Z; 0) and x( (b): The LL, CLL, HLL on the test set, and the average parameter error of the weights and biases the linear mapping. (c): True and estimated parameters by different methods of the first trial. For each matrix, the leftmost column is the bias b, and the remaining block the weight matrix W. The top-left block of the weight part represents visible-to-visibee the top-right block represents hidden-to-visible, the bottom-left block represents visible-to-hiddee
and the bottom-right block represents haldeen-t-hiddeen (d): Predictive firing rates on spike train from different methods. Specifically, given complete test spike train Y = X Z]. we can predict the firing rates by the complete model p(X, Z: 0) via Eq, [4 for both visible neurons (e.g., neuron 1) and hidden neurons (e.g., neuron 4). For hidden neurons (e.g., neuron 4), we can also predict the firing rates by q( x (x:)).
weights (the bottom-left block of the weight part). VI, CHIVI, and VBIS provide random-like and non-informative estimations, but VIS matches the true better. For the hidden-to-hiddee weights (the bottom-right block of the weight part), none of the four methods gives acceptable results. The worse performances on the hidden-to- to--iiibll and midden-to-hiddee blocks also reflect the limitation of the variational/proposaal distribution family.
In Fig. d), we visualize the predictive firing rates fi. learned by different methods. The top panel and the middle panel of Fig. 4d) show that the firing rates predicted by p(X, Z:0) obtained from VIS for both visible neurons and hidden neurons are the most accurate the true firing rates among all four methods. Particularly, since only VIS learns acceptable visible-to-hidder weights, the firing rates predicted by VI, CHIVI, and VBIS are significantly worse than by VIS (the middle panel of Fig. 4()). These correspond to the CLL bar plot in Fig. 4b. The bottom panel of Fig. dd) indicates that the proposal distribution of VIS can sample hidden spikes much closer to the true hidden spikes, which improves the learning effects and results in better parameter recovery. Moreover, methods except VIS in the middle panel and the bottom panel of Fig. Rdd) reveal the case that q( Z | f)) and p(Z Z x) close in terms of the reverse KL divergence, but both of them are far from the true posterior, resulting in higher ELBO but lower marginal log-likelihood than VIS.
4.3.2 RETINAL GANGLION CELL (RGC) DATASET
real neural spike train recorded from = 27 retinal ganglion neurons while mouse is performing a visual test for about 20 mins (Pillow & Scott 2012]. Neurons 1-16 are OFF cells, and neurons 17-27 are ON cells.
Dataset. We run different methods
Experimental setup. We use the first 3 segment as the training set and the remaining 3 segment the test set. The original spike train is converted to spike counts in every 50 ms time bins. For applying the stochastic gradient descent algorithm, we break the whole sequence into several pieces. The length of each piece is 100 time bins. First, we learn fully observed GLM baseline. Then,
1
8

---

Published as a conference paper at ICLR 2024
(a)
-1512
hidden
hidden
hidden
0.30
-1514
-1516
00.5
fully
0010
onserved
GLM
0.15
predicivee D
VI
CHVI
VGIS
VIS
0.25
0.2
75
1.00
75 1.00 time
100
100
Figure 5: (a): The marginal log-likelihooe on the test segment with different numbers of hidden neurons. (b): The estimated weight matrices from different methods. (c): 20 predictive firing rates generated from hidden spikes sampled from different variational/proposall distributions.
assume there are H E (1,2,3) hidden representative neurons and learn the POGLM by different methods. We use Adam (Kingma & Ba 2014) as the optimizer and the learning rate is set at 0.01. We run 10 epochs for each method. The batch size is set as 32. The number of Monte Carlo samples used for sampling the hidden are 1,000, 2,000, and 3,000 for H = 1,2,3 respectively. We repeat 10 times with different random seeds for each method and report the performance.
Results. Compared with the fully observed GLM (the dashed line in Fig. a)), adding hidden neu- rons significantly improves the capability of predicting spiking events on the test set, when learned by VBIS and VIS. This is reflected in the high test marginal log-like ledliik of VBIS or VIS shown in Fig. aJa). Particularly, VIS always obtains the highest test marginal log-likelihood compared with the three alternative methods.
We also visualize the learned weight matrix with one hidden neuron from the four methods in Fig. ().). With one hidden neuron learned by VIS, the weights from the hidden neuron nearly all OFF cells are positive, and the weights to all ON cells are negative. This implies that this hidden representative neuron behaves like an OFF cell. The signs of the weights from this hidden represen- tative neuron to the visible neurons clearly tell us the type of those visible post-synaptic neurons. All other methods do not have such a significant differentiation in the last column of the weight matrix.
Since we do not have the true hidden spike train in the real-world dataset, we sample hidden spike trains from the variationallproposaal distribution q( X Φ), and compute the corresponding firing rates that are used for sampling the hidden spike trains. In Fig. (j)). we plot 20 randomly sampled predictive firing rates of the hidden neuron in the one-hidden-neurce model. Clearly, the predic- tive firing rates generated by VIS provide wider effective support range for sampling, due to the mass-covering/maan-eeee weanngg//ean-eeeinnnggn seeking behavior of minimizing the forward x² divergence. This variability improves the effectiveness of learning In p(X;0). Compared with VIS, the variational/proposall distributions learned by VI and VBIS are very restricted and concentrative, providing less variabil- ity in sampling hidden spikes. Since CHIVI minimizes both the forward x² and the reverse KL divergences, the variability of the variational/proposall distribution is at a medium position.
5 DISCUSSION
In this paper, we introduce variational importance sampling (VIS), a novel method for efficiently learning parameters in latent variable models, based on the forward x² divergence. Unlike varia- tional inference (VI), which maximizes the evidence lower bound (ELBO), VIS directly estimates and maximizes the marginal log-likelihood to learn model parameters. Our analyses demonstrate that the quality of the estimated marginal log-liiel is assured with a large number of Monte Carlo samples and optimal proposal distribution characterized by small forward x² divergence. This highlights the statistical significance of choosing the proposal distribution. Experimental results across three different models validate VIS's ability to achieve both higher marginal log-lik likeelhooood and a better parameter estimation. This underscores VIS as a promising learning method for ad- dressing complex latent variable models. Nevertheless, it is worth noting that while this choice of the proposal distribution is statistically optimal for importance sampling, its practical significance in certain real-world applications might require further investigation and validation.

---

Published as a conference paper at ICLR 2024
REFERENCES
Ömer Deniz Akyildiz and Joaquín Míguez. Convergence rates for optimised adaptive importance samplers. Statistics and Computing. 31:1-17, 2021.
David M Blei, Alp Kucukelbir, and Jon D McAuliffe. Variational inference: A review for statisti- cians. Journal of the American statistical Association, 112(518):859-8777 877, 2017.
Yuri Burda, Roger Grosse, and Ruslan Salakhutdinov. Importance weighted autoencoders. arXiv
preprint 0w0./0090000500 000900 2015.
Adji Bousso Dieng, Dustin Tran, Rajesh Ranganath, John Paisley, and David Blei. Variational inference via X upper bound minimization. Advances Neural Information Processing Systems, 30, 2017.
Justin Domke and Daniel R Sheldon. Importance weighting and variational inference. Advances in neural information processing systems, 31, 2018.
Axel Finke and Alexandre H Thiery. On importance-weighted autoencoders. arXiv preprint vrXiv:/907.1047) 2019.
David Freedman, Robert Pisani, and Roger Purves. Statistics. w. w, 1998.
Tomas Geffner and Justin Domke. On the difficulty of unbiased alpha divergence minimization. arXiv preprint arXiv:2010.09541, 2020.
Jose Hernandez-LLobtto Yingzhen Li, Mark Rowland, Thang Bui, Daniel Hernández-Lobatte and Richard Turner. Black-box alpha divergence minimization. In International conference on ma- chine learning, pp. 151 1511-0. PMLR, 2016.
Ghassen Jerfel, Serena Wang, Clara Wong-Famnn Longg-annjiang, katherine A Heller, Yian Ma, and Michael I Jordan. Variational refinement for importance sampling using the forward kullback-leibler diver- gence. Uncertainty in Artificial Intelligence, pp. 1819-1829. PMLR, 2021.
Hiroshi Kajino. A differentiable point process with its application to spiking neural networks. In International Conference Machine Learning, pp. 5226-5235. PMLR, 2021.
Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv preprint SE1i::1412.6980 2014.
Diederik P Kingma and Max Welling. 1:Xiv:1312.6114.200 2013.
Auto-encoding variational bayes. arXiv preprint
Teun Kloek and Herman K Van Dijk. Bayesian estimates of equation system parameters: an ap- plication of integration by monte carlo. Econometrica: Journal of the Econometric Society, pp. 1-19, 1978.
Yann LeCun, Léon Bottou, Yoshua Bengio, and Patrick Haffner. Gradient-based learning applied to document recognition. Proceedings of the IEEE, 86(11):227--232, 1998.
Yingzhen Li and Richard Turner. Variational inference with rényi divergence. Statistics, 1050, 2016.
Tomohiro Nishiyama and Igal Sason. On relations between the relative entropy and x 2-divergence, generalizations and applications. Entropy, 22(5):563, 2020.
Gary W Oehlert. A note the delta method. The American Statistician, 46(1):27-29, 1992.
Jonathan Pillow and Peter Latham. Neural characterization in partially observed populations of spiking neurons. Advances Neural Information Processing Systems, 20, 2007.
Jonathan Pillow and James Scott. Fully bayesian inference for neural models with negative-binomial binomial spiking. Advances neural information processing systems, 25, 2012.
10

---

Published as a conference paper at ICLR 2024
Jonathan W Pillow, Jonathon Shlens, Liam Paninski, Alexander Sher, Alan M Litke, EJ Chichilnisky, and Eero P Simoncelli. Spatio-temporal correlations and visual signalling in complete neuronal population. Nature, 454(7207):995-999 999,99999, 2008
Melanie F Pradier, Michael C Hughes, and Finale Doshi-Velez. Challenges in computing and opti- mizing upper bounds of marginal likelihood based on chi-square divergences. In Second Sympo- sium on Advances in Approximate Bayesian Inference, 2019.
Danilo Jimenez Rezende and Wulfram Gerstner. Stochastic variational learning recurrent spiking networks. Frontiers in computational neuroscience, 8(ARTICLE):38 2014.
Ram Naresh Saraswat. Chi square divergence measure and their bounds. In 3rd International Conference on "Inovvative Approach in Applied Physical, MuthematicalSSatisiiaal" Chemical Sciences and Emerging Energy Technology for Sustainable Development, pp. 55, 2014.
Igal Sason and Sergio Verdú. f-divergence inequalities. IEEE Transactions on Information Theory, 62(11))5973-6006, 2016.
John Schulman, Nicolas Heess, Theophane Weber, and Pieter Abbeel. Gradient estimation using stochastic computation graphs. Advances in neural information processing systems, 28, 2015.
Eukasz Struski, Marcin Mazur, Pawel Batorski, Przemystaw Spurek, and Jacek Tabor. Bounding evidence and estimating log- lielikkliihodd in vae. arXiv preprint S1kiv:220060945555 2022.
Xiao Su and Yuguo Chen. Variational approximation for importance sampling. Computational Statistics,
2021.
Yuling Yao, Aki Vehtari, Daniel Simpson, and Andrew Gelman. Yes, but did it work?: Evaluating variational inference. In International Conference on Machine Learning, pp. 5581-5590, PMLR, 2018.

---

Published as a conference paper at ICLR 2024
A
APPENDIX
A.1
GRADIENT ESTIMATORS THE VARIATIONAL. INFERENCE
The derivative of (L; 0, Φ) w.r.t 0 is estimated by
(L ax; 0,
i In p( 2; 0)
[(2|2;00) o) dæ
00
00
K
Oln p (x, (R)
K
00
(16)
lpp σ, 00)
00 K
For the derivative of FL:OO(;;,0)) w.r.t Φ at Φu. the score function gradient estimator is i EEBOO(;;6) 0, σ)
In [n,p(rr 2; - [2/ 60)1 0q(=|<; Φ)
Iln (x/ Φ)
0Φ g( (o)
dz
0Φ
OΦ
/ In p(z, 2; - -n 60) (<||= 0o)
i In [=|[;; Φ)
dz
g(= () dz
dΦ
0 In (z [x;
Inp x; 20))) Inn (= x; (0))
K
Og
K
Imp 2,, =0)00 -nn (2(0) (x; 6)
86
2K
k=1
(17)
When the parameterization trick can be utilized, a;; = g(c|a;; () where e ~ r(c), then
g<||<; Φ) dz = r(E) de.
(18)
Now, get the pathwise gradient estimator,
i FLLBO 0, s)
/ ((<l=;0) [np(( 2; In (= の)] dz
OΦ
9
/ r(c) In (( Φ)) (1]0((2/1 de
0Φ
(19)
0
Σ
In x (=)))=:0): Im (g (c (z: 6)==))=0))))
[=:
dΦ
K
A.2 GRADIENT ESTIMATOR OF THE IMPORTANCE SAMPLING
The derivative of p(: 0) w.r.t. at Op is estimated by
3 In 1hpp 0)
Öp(x z: 0) dz
00
p(x:;0)
00
1 11pp(r, 0)
p(22,20 2; 00)
dz
p(=;; 0o)
00
p (n, z (k); Øln p (x, (k);
1
f( 0o) K
4 (z(4)
x ②0))
20
(20)
K
exp In
x. (+) 0) In (2(4) x:0)]
p(x; 0o) 80 K
k=1 0p(2; 8)
I 0mp((x; 8)
p(x;;0 0o)
00
20
12

---

Published as a conference paper at ICLR 2024
Due to the appearance of p(æ;; Bb) in the denominator, (l Up
is magnitude up-biased estimator
of 0 In p(=:0) However, the direction of 0 ((x00,0) is unbiased:
0e
pp(x;; 0,
1
Op (2, (b),
B
B.g
00
q
(z)
x;;
00
k=
1
Öp(æ, 2:0)
Öp(æ, 2: z)
(21)
By
dz
((z x 6)
00
00
/
0p(x; 0)
p(xx 2; dz
00
00
A.3 GRADIENT ESTIMATOR FOR UPDATING THE PROPOSAL DISTRIBUTION IN VIS
In this section, we derive the score function gradient estimator and the pathwise gradient estimator for minimizing the forward x² divergence, which equivalent minimizing in V(æ; Φ) in Eq I1 First, we show the derivation of Eq.
K
p (x, (K); 0)
1 K
In V(x; 0) () R In
<x; ()) k=
1 (z(k)
(22)
2 ln p (x, (0) 2 In (2 xx
=logsumexp
In K
=: (x; O).
The score function gradient estimator of V (x;0, 0, Φ) in Eq. is i In V(:; 0, 6)
p( z; 0)2
dz
Osb
Vax; 0, Φ0)
dΦ q( a Φ)
pla,, 2; 0)2 a
(=/ x Φ) dz
V V(x;; 0, 60)
(1 Φ0) 06
1
K
p (x, (k) 0)
Oln y (z (8)
a;; s))
V (x; 0, )0))
K
(x)
x: 0o)²
(23)
k=1
1
1 8
2np (x, (b))
In In 1=00)200)) (x:0)
exp
V(x; 0, 20)
K
2 06
k=1
0 1 In (z:
0Φ 2
When the reparameterization trick can be utilized. x <; ó g(c 0(cx::0 Φ) where e ~ r(E), then we have the transformation <=>|<;; Φ) dz r(e) de (Schulman al. 2015).
Then.
In V(x; 0,
a / ((<2|1) x; )b) p(z, 2: 8)2
((>|=; o)² dz
dΦ
V(x; O, 20) ds
V(x; 0, (0) 06 /
r(c) p(ar 2; 2)
a 0)2 de
(24)
K
exp
2 Imp
(z, (0)
2 22n 2ing(2() x; 6)]
V(x; 0, ②) 0Φ K
k=1
0
In (x:; ; 0 0)
06
13

---

Published as a conference paper at ICLR 2024
A.4 LATENT MANIFOLD OF THE MNIST DATASET
The following figures are the latent manifolds of the MNIST dataset learned by different methods.
manifold learned by VI
manifold learned by CHIVI
7
manifold learned by VBIS
manifold learned by VIS
7 7
7
manifold learned by IWAE
Figure 6: Latent manifolds of the MNIST dataset learned by different methods.
14

---

Published as a conference paper at ICLR 2024
A.5 COMPARISON OF DIFFERENT GRADIENT ESTIMATORS OF EQ.
Considering that the numerical issue minimizing the forward X² divergence is widely discovered
by lot of previous works (Pradier et al. 2019; Finke & Thiery 2019 Geffner & Domke 2020
Yao et al. 2018). we run VIS on the toy mixture model again (Sec. [4.1) using the {score function, pathwise gradient estimator in {log, original} space for minimizing the forward x² divergence. Results in Fig. / show that the score function gradient estimator is better than the pathwise gradient estimator for minimizing the forward x² divergence. Besides, it is important to estimate it in log space so that the numerical stability of the score function gradient estimator can be promised.
(a) (a)-0.50
(b)
original
log
original
1.00
1.25
-10
10-
-5
0.65
0.65
0.655
Ice 019)
ong
0.60
3.60
0.60
0.60
0.50
0,50
0.50
0.50
100
100
200
40
epochs
Posteriar
Posterior original
0.4
Posterior pethwise log
Posterior pathwise
0.4
15
giz
Figure 7: (a): LL, CLL, and HLL evaluated on the test dataset. (b): Convergence curves of the parameter set learned by different gradient estimators. The dashed curves are the true parameters used for generating the data, and the solid curves are the learned parameters. (c): The posterior distribution given x 0 and x learned by different gradient estimators. The dashed curves are the true posterior p(=|2;; (fre),, the solid curves are the learned posterior p(=|2;; 0), and the dotted curves are approximated posterior g(= [x; Φ) learned in the variational/proposall distribution.
15

---

Published as a conference paper at ICLR 2024
A.6 RUNNING TIME OF DIFFERENT METHODS
Fig. 8 shows the test LL and corresponding running time of different methods w.r.t. the number of Monte Carlo K the synthetic POGLM dataset (Sec. 4.3). In general, the running times of all methods are linear to the number Monte Carlo samples. With more Monte Carlo samples, all methods perform better, and VIS is consistently better than others especially when is large. When K is small, all methods fail because of the complex nature of the POGLM problem. This implies that for complicated graphical models and high dimensional latent space, we need enough Monte Carlo samples for all these sampling-ber based methods to become effective. Therefore, the number of Monte Carlo should be suitable to the complexity of the model/problem, rather than which method we choose.
140
method
-215
120
CHV
marginal Ing-literreeee 0.1
217
80
219
40-
221
0
250
500
750
1000
1250
1500
1750
2000
0
250
500
750
1000
1250
1500
1750
2000
number Monte Carlo K
Figure 8: Test LL (left) and corresponding running time (right) of different methods w.r.t. the number Monte Carlo K, on the synthetic POGLM dataset (Sec.
16

---

Published as a conference paper at ICLR 2024
A.7 FORWARD KL DIVERGENCE
(Jerfel et al. 2021) considers forward KL divergence as the target function for updating the proposal distribution since they noticed the drawback of the reverse KL divergence. According to (Sason &
Verdú 2016) and (Nishiyama & Sason 2020), however, KL(pllg) can be bounded by x (pll
KL(pll9) & In (plle)) < (pll))
(25)
but not vice versa. Therefore, minimizing the forward KI. divergence might not be able to get the optimal proposal distribution, which should be obtained by minimizing the forward x² divergence.
To validate this empirically, we compare minimizing the forward x² divergence (VIS) to minimizing the forward KL divergence (forward KL. on the toy mixture model again (Sec. and the results
are shown in Fig. 9
(a) -0.629
(b)
VIS
forward
SI
oowward KL
10
6300
-25
-0,631
-0.6302
0 -35
E 2.0
0.70
0.70
0.85
-0663
0.60
0.55
0.55
0.6304
4.
50
100
150
200
50
100
150
200
0.6305
epochs
0.1
6.0
-15
-10
10
-15
-10
-5
piz x D: Deghi
..... QUT. 20rpe 0,
(paj) 1:0
01/14 - 1:
plr/ 7:40)
Figure 9: (a): LL. CLL, and HLL evaluated on the test dataset. (b): Convergence curves of the parameter set 0 learned by VIS and forward KL. The dashed curves are the true parameters used for generating the data, and the solid curves are the learned parameters. (c): The posterior distribution given and x 1 learned by different gradient estimators. The dashed curves are the true posterior p(2[2: [2 gtmee), the solid curves are the learned posterior p(2|1:; 8), and the dotted curves are the approximated posterior q( [z|=;; Φ) learned in the variational/proposal distribution.
17

---

Published as a conference paper at ICLR 2024
A.8 RELATED WORKS AND CONTRIBUTIONS TABLE
Here, we aim to offer concise summary of our contributions and related works.
Table I: Contributions.
Contributions
Previous literatures
Motivate from the effectiveness of IS
[3] [5] [6] [7]
Aim learning A
[1] [3] [4] [5] [6] [7]
No restrictions the distribution families
[1] [2] [3]
Directly minimizing forward x² divergence without surrogate
[2] [3] [7]
Motivate from the bias of IS in log space
Numerically stable gradient estimator in log space
Extensive experiments cases where explicit decomposi-
tion p(æ, z; 0) = p(æ|2; p))((; 0) exists
Visualization for inferred latent and parameter 0's recovery
Motivate from the bias of IS in log space: We start by comparing the bias of the In p(x; 8, Φ) and (L:Ö 0, Φ) to analyze why doing IS and the optimal way of doing IS. And the conclusion about minimizing the forward x² divergence coincides with improving the effectiveness of the IS estimator (Fig [
Numerically stable gradient estimator in log space: Previous work already derived the gradi- ent estimator for minimizing the x² divergence in the original space but not log space. This leads to the numerical instability issue and scaling to the high dimensionality issue. We argue that it is critical to estimate its gradient in log space to obtain a numerically stable and succinct form of the gradient estimator, especially for the score function estimator (Fig. 77.
Extensive experiments on cases where no explicit decomposition p(x: z)
(()=)=;)))(;;0 0) exists: Most of the previous work only do experiments generative models with explicit decomposition p(æ, z; 0) = p(æ) z; ();(;; f), unlike the POGLM. However, when such an explicit decomposition does not exist or when the generative posterior distribution p(z æ;; 0) and the approximating posterior distribution q( ax; Φ) are not Gaussian, ELBO cannot be reformulated as (=:O 0, Φ) E, p(z zz f)] (g. =| [)p(p( f)). and hence ELBO lost lot of ad- vantages. Therefore, we do need a variety of graphical models to understand the performance of different methods.
Visualization for inferred latent and parameter 0's recovery: Although theoretical materials in this paper have shown the superiority of VIS, practical visualization of the behavior methods is still necessary for us get intuition of how and why VIS performs better than others.
of different
[1] Burda et al.
(0055
[2]
Dieng et al.
2017
[3]
Finke & T
[4]
Jerfel et al.
1)
[5]
Domke & Sheldon
2018
[6
Su & Chen (2021)
[7]
Akyildiz & Miguez
2202
18


Reviews:
[{"review_id": "g0QgmPx7E4", "reviewer": "Reviewer_84Pi", "summary": "The paper proposes a variational algorithm for learning model and latent parameters in a latent variable model which at each step first updates the variational distribution by minimizing forward chi-square divergence and then uses this distribution to estimate the log marginal likelihood using Importance Sampling. The optimization is done by gradient ascent where the gradients are estimated by MC sampling.\nThe proposed algorithm is compared against many contemporary algorithms on simulated and real world datasets, including a large scale case study on multi-neuron interaction modelled by a custom made partially observable GLM.", "questions": "1. What is the dimensionality of the POGLM model, do the authors intend to use this method as a tool for low dimensional complex posteriors because both IS and CUBO do not scale well with dimensions and even large sample size as done in this paper will not help.  \n2. Maybe include this in your conclusion section and discuss this as a limitation ? \n3. Did the authors use any other optimizers other than ADAM, did it have any effect, how did you choose the optimization algorithm hyperprameters like learning rate etc. ?\n4. Did reparameterization gradients perform better than score gradients in the case where they both were available.", "limitations": "", "rating": 5, "confidence": 3, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "1. The paper tries well to motivate itself well and the use of chi-square divergence objective as means of finding the optimal distribution for IS is theoretically sound.\n2. The paper has done experiments and analysis on may simulated and real world datasets. The experiment on GLP model for multi-neurns activation is well documented and insightful. Some of the plots look good and match the narrative. \n3. The method proposed seems sound to me and the results show that it can perform better than contemporary VI algorithms on the tasks given in the paper.", "weaknesses": "1. The paper is not polished yet, although the major parts are all there, it may require another thorough pass.\nit has too many mistakes and typos:, the notation changes from bold to normal in many places, the title has a typo: 'importane', 'log function is a convex function'. \n2. Some of the references and recent literature is missing which have looked on the quality of different divergence objectives such as CUBO and ELBO for finding the optimal sampling distribution.\n3. The theory part and the algorithm part can be emphasized more, right now it feels to compressed and dense. The figure 2 is good but it has too many colors and things to unpack, maybe use solid line for true posterior as I was thoroughly confused by the legend choices, and the use of two colors for showing modes reduced readibility for me atleast. \n4. It is the bane of chi-square divergence methods that it  does not scale well with dimensions covered in the papers here: https://arxiv.org/pdf/2010.09541.pdf and https://arxiv.org/abs/1802.02538 and it seems that this method may not scale well as it uses Chi-squared divergence minimization.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "extra_scores": {}}, {"review_id": "ztkaEzFhFe", "reviewer": "Reviewer_J8Yz", "summary": "This work proposes an adaptive importance-sampling algorithm which seeks to improve the proposal distribution by minimising the chi-square divergence from the target distribution to the proposal.", "questions": "How does the approach compare with replacing the forward $\\chi^2$- with forward KL-divergence (i.e., effectively two out of the three \"phases\" of reweighted wake--sleep)?", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "There are a large number of experiments. The proposed method seems to perform well.\n\nApproximating the gradient of the logarithm of the forward $\\chi^2$-divergence (rather than the gradient of the $\\chi^2$-divergence) seems to enhance numerical stability and seems novel. Though it should be noted that existing approaches improve numerical stability by minimising not the forward $\\chi^2$-divergence itself but rather the forward $\\chi^2$-divergence multiplied by the squared normalising constant of the target distribution: $p_\\theta(x)^2$ (see, e.g., [2]). It is not clear how this compares to the log-space approach taken here.\n\n[2] Akyildiz, Ö. D., & Míguez, J. (2021). Convergence rates for optimised adaptive importance samplers. Statistics and Computing, 31, 1-17.", "weaknesses": "Targetting the $\\chi^2$-divergence as an objective for improving the proposal distribution in adaptive importance sampling is not novel. The fact that minimising the variance of the importance weights is equivalent to minimising the chi-square divergence from the target to the proposal is well known in the importance-sampling literature and has already often been used to improve proposal distributions within adaptive importance-sampling schemes, e.g. [1, 2] and references therein.\n\nFurthermore, using such adaptive-importance-sampling approaches for variational inference is already extensively discussed in [3].: \n1. Algorithm 1 of the present paper is a special case of the generic method described in [3, Section 3] (in particular, see [3, Section 3.4]);\n2. the $\\theta$-gradient from Equation 6 of the present work is already well known (see, e.g. [3]).\n3. However, from the author's rebuttal, it is now more clear to me that their $\\phi$ gradient is slightly different than in [3] because they derive the $\\phi$-gradient in log-space.\n\n[1] Jona‐Lasinio, G., Piccioni, M., & Ramponi, A. (1999). Selection of importance weights for monte carlo estimation of normalizing constants. Communications in Statistics-Simulation and Computation, 28(2), 441-462.\n\n[2] Akyildiz, Ö. D., & Míguez, J. (2021). Convergence rates for optimised adaptive importance samplers. Statistics and Computing, 31, 1-17.\n\n[3] Finke, A., & Thiery, A. H. (2019). On importance-weighted autoencoders. arXiv preprint arXiv:1907.10477.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "extra_scores": {}}, {"review_id": "XYRRTb5XR7", "reviewer": "Reviewer_QxSX", "summary": "The paper focuses on the problem of approximate Bayesian inference in latent variable models and challenges the common wisdom of approximating the posterior with variational inference and maximizing the evidence lower bound (ELBO). The authors identify a mismatch between maximizing ELBO and maximizing loglikelihood, and posit that the latter is a better objective, as it is divorced from the quality of the approximate posterior. To tackle the maximization of loglikelihood, the paper proposes variational importance sampling (VIS), which maximizes the marginal directly via importance sampling. Upon inspecting the variance of the IS estimator, the paper also proposes minimizing the chi-squared divergence w.r.t. the approximate posterior. \n\nThe paper inspects the efficacy of this method on a series of latent variable models, and witness a consistent improvement on the inference of latent variable models.", "questions": "- The experiments presented in the paper show that VIS performs well, but it seems that a large number of Monte Carlo samples are chosen in many of the experiments. This seems that it can present a significant computational cost. Could the authors include some explanation about the effect of the number of Monte Carlo samples on the efficacy of VIS, and clarify if taking many Monte Carlo samples is practical in training latent variable models?", "limitations": "", "rating": 8, "confidence": 3, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "I think the paper warrants acceptance due to the fact that it tackles a well-motivated problem with a simple, yet empirically effective method. \n\n- Motivation: the maximization of ELBO as a proxy of likelihood maximization has long been a standard practice in latent variable, even though the mismatch can be significant. I agree with the authors that the mismatch should be inspected more carefully, and that one should separate the approximate posterior inference with the model itself. \n- Methodology: the paper proposes a simple, yet elegant solution to the question of likelihood maximization, and motivate the reasoning behind the optimization of the chi-squared divergence from the perspective of the bias and variance of the IS estimator. This approach seems novel, even though I am not very up-to-date about the research in this regard. \n- Experiments: the paper compares against alternative options in latent variable inference and showcase that the method can correctly infer the marginal likelihood, as well as outperforming competing methods in real-world datasets.", "weaknesses": "I am not entirely sure how the method proposed in this paper differs from previous work, as the practice of using Monte Carlo samples to sharpen variational bound is a topic that has been explored by previous work. I hope the authors can make some clarifications or have a small section in the manuscript that discusses the difference across different ways to combine IS with VI. \n\n- It is neat to see the same chi-squared divergence on both the bias and effectiveness of the IS estimator, but eq. 8 contains an approximation that seems to work mostly for large $K$s, and the effectiveness is evaluated at the estimator $\\hat{p}$ and not $\\log \\hat{p}$, making the connection seem a bit artificial. Is it possible to minimize a divergence in the form of, e.g., the first line in eq. 8?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "extra_scores": {}}, {"review_id": "m8TuK97Wz0", "reviewer": "Reviewer_T2mb", "summary": "This work proposes a new method of variational inference (VI) for latent variable models based on importance sampling, and shows its effectiveness with numerical experiments on toy model (synthetic), auto encoder model (real) and spike neural network (synthetic and real).\n\nThe method improves the tightness of the conventional VI by increasling the batch size of monte-carlo sampling, while giving up the exact computation/minimization of the approximation gap introduced with VI (cf. Eq. 8 and 11). This theoretical trade-off turned out to be beneficial in the experiments.", "questions": "The ground truth model of the toy experiment (Sec. 4.1) seems unidentifiable. In particular, the distribution of the hidden variable has multiple modes, which seems impossible to be identified only with the binary visible variable. This observation also explains well that despite the fact that the conventional VI substantially failed to estimate the modes, there is only a tiny difference in log-likelihood compared to other methods (order of 1e-4, see Fig. 2a). So, the question is, **why do you compare the parameter convergence, HLL and CLL**, which I think are meaningful to compare only if the model is identifiable? A similar argument also applies to the synthetic experiment of spike NN.", "limitations": "", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "The paper is original and well written. The performance gain obtained with the proposed method is significant, which is one of the main contribution of the paper.", "weaknesses": "There is no discussion on the limitation of the proposed method.\nFrom what I understood from the paper, I expect the following potential drawbacks:\n  - Increase in training time based on the additional sampling\n  - Instability in training due to the biased approximation (due to Eq. 8 and 11).\n\nDetailed discussions on this matter are highly welcomed.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (spotlight)