Amortized Variational Inference in Simple Hierarchical Models
Year: 2022

Abstract:
It is difficult to use subsampling with variational inference in hierarchical models since the number of local latent variables scales with the dataset. Thus, inference in hierarchical models remains a challenge at large scale. It is helpful to use a variational family with structure matching the posterior, but optimization is still slow due to the huge number of local distributions. Instead, this paper suggests an amortized approach where shared parameters simultaneously represent all local distributions. This approach is similarly accurate as using a given joint distribution (e.g., a full-rank Gaussian) but is feasible on datasets that are several orders of magnitude larger. It is also dramatically faster than using a structured variational distribution.

heading: 1 Introduction
text: Hierarchical Bayesian models are a general framework where parameters of “groups” are drawn from some shared distribution, and then observed data is drawn from a distribution specified by each group’s parameters. After data is observed, the inference problem is to infer both the parameters for each group and the shared parameters. These models have proven useful in various domains [13] including hierarchical regression amd classification [12], topic models [4, 22, 3], polling [11, 24], epidemiology [23], ecology [8], psychology [37], matrix-factorization [35], and collaborative filtering [26, 33].
A proven technique for scaling variational inference (VI) to large datasets is subsampling. The idea is that if the target model has the form p(z, y) = p(z) ∏ i p(yi|z) then an unbiased gradient can be estimated while only evaluating p(z) and p(yi|z) at a few i [29, 16, 20, 31, 30, 36, 15]. This paper addresses hierarchical models of the form p(θ, z, y) = p(θ) ∏ i p(zi, yi|θ), where only y is observed. There are two challenges. First, the number of local latent variables zi increases with the dataset, meaning the posterior distribution increases in dimensionality. Second, there is often a dependence between zi and θ which must be captured to get strong results [15, 18].
The aim of this paper is to develop a black-box variational inference scheme that can scale to large hierarchical models without losing benefits of a joint approximation. Our solution takes three steps. First, in the true posterior, the different latent variables zi are conditionally independent given θ, which suggests using a variational family of the same form. We confirm this intuition by showing that for any joint variational family q(θ, z), one can define a corresponding "branch" family q(θ) ∏ i q(zi|θ) such that inference will be equally accurate (theorem 2). We call inference using such a family the "branch" approach.
Second, we observe that if using the branch approach, the optimal local variational parameters can be computed only from θ and local data (eq. (12)). Thus, we propose to amortize the computation of the local variational parameters by learning a network to approximately solve that optimization. We
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
show that when the target distribution is symmetric over latent variables, this will be as accurate as the original joint family, assuming a sufficiently capable amortization network (claim 5).
Third, we note that in many real hierarchical models, there are many i.i.d. data generated from each local latent variable. This presents a challenge for learning an amortization network, since the full network should deal with different numbers of data points and naturally reflect the symmetry between the inputs (that is, without having to relearn the symmetry.) We propose an approach where a preliminary "feature" network processes each datum, after which they are combined with a pooling operation which forms the input for a standard network (section 6). This is closely related to the "deep sets" [39] strategy for permutation invariance.
We validate these methods on a synthetic model where exact inference is possible, and on a userpreference model for the MovieLens dataset with 162K users who make 25M ratings of different movies. At small scale (2.5K ratings), we show similar accuracy using a dense joint Gaussian, a branch distribution, or our amortized approach. At moderate scale (180K ratings), joint inference is intractable. Branch distributions gives a meaningful answer, and the amortized approach is comparable or better. At large scale (18M ratings) the amortized approach is thousands of nats better on test-likelihoods even after branch distributions were trained for almost ten times as long as the amortized approach took to converge (fig. 6).
2 Hierarchical Branched Distributions
We focus on two-level hierarchical distributions. A generic model of this type is given by
p(θ, z, y|x) = p(θ) N∏ i=1 p(zi|θ)p(yi|θ, zi, xi), (1)
where θ and z = {zi}Ni=1 are latent variables, y = {yi}Ni=1 are observations, and x = {xi}Ni=1 are covariates. As the visual representations of these models resemble branches, we refer them as hierarchical branch distributions (HBDs).
Symmetric. We call an HBD symmetric if the conditionals are symmetric, i.e., if zi = zj , xi = xj , and yi = yj , it implies that
p(zi|θ) = p(zj |θ), and p(yi|θ, zi, xi) = p(yj |θ, zj , xj). (2)
Locally i.i.d. Often local observations yi (and xi) are a collection of conditionally i.i.d observations. Then, an HBD takes the form of
p(θ, z, y|x) = p(θ) N∏ i=1 p(zi|θ) ni∏ j=1 p(yij |θ, zi, xij), (3)
where yi = {yij}nij=1 and xi = {xij} ni j=1 are collections of conditionally i.i.d observations and covariates; ni ≥ 1 is the number of observations for branch i.
No local covariates. Some applications do not involve the covariates xi. In such cases, HBDs have a simplified form of
p(θ, z, y) = p(θ) N∏ i=1 p(zi|θ)p(yi|θ, zi). (4)
In this paper, we will be using eq. (1) and eq. (3) to refer HBDs—the results extend easily to case where there are no local covariates. (For instance, in section 5, we amortize using (xi, yi) as inputs. When there are no covariates, we can amortize with just yi.)

heading: 2.1 Related Work
text: Bayesian inference in hierarchical models is an old problem. The most common solutions are Markov chain Monte Carlo (MCMC) and VI. A key advantage of VI is that gradients can sometimes be estimated using only a subsample of data [29]. Hoffman et al. [16] observe that inference in hierarchical models is still slow at large scale, since the number of parameters scales with the dataset. Instead, they assume that θ and zi from eq. (1) are in conjugate exponential families, and observe that for a mean-field variational distribution q(θ) ∏ i q(zi), the optimal q(zi) can be calculated in closed form for fixed q(θ). This is highly scalable, though it is limited to factorized approximations and requires a conditionally conjugate target model.
A structured variational approximation like q(θ) ∏
i q(zi|θ) can be used which reflects the dependence of zi on θ [15, 34, 2, 18]. However, this still has scalability problems in general since the number of parameters grows in the size of the data (section 7). To the best of our knowledge, the only approach that avoids this is the framework of structured stochastic VI [15, 18], which assumes the target is conditionally conjugate, and that for a fixed θ an optimal "local" distribution q(zi|θ) can be calculated from local data. Hoffman and Blei [15] address matrix factorization models and latent Dirichlet allocation, using Gibbs sampling to compute the local distributions. Johnson et al. [18] use amortization for conjugate models but do not consider the setting where local observations are a collection of i.i.d observations. Our approach is not strictly an instance of either of these frameworks, as we do not assume conjugacy or that amortization can exactly recover optimal local distributions [15, Eq. 7]. Still the spirit is the same, and our approach should be seen as part of this line of research.
Amortized variational approximations have been used to learn models with local variables [20, 10, 17, 5]. A particularly related instance of model learning is the Neural Statistician approach [10], where a “statistics network” learns representations of closely related datasets–the construction of this network is similar to our “pooling network” in section 6. However, in the Neural Statistician model there are no global variables θ, and it is not obvious how to generalize their approach for HBDs. In contrast, our “pooling network” results from the analysis in section 5, reducing the architectural space when θ is present while retaining accuracy. Moreover, learning a model is strictly different from our “black-box" setting where we want to approximate the posterior of a given model, and the inference only has access to log p or ∇θ,z log p (or their parts) [30, 21, 2, 1].

heading: 3 Joint approximations for HBD
text: When dealing with HBDs, a non-structured distribution does not scale. To see this, consider the naive VI objective—Evidence Lower Bound (ELBO, L). Let qJointϕ be a joint variational approximation over θ and z; we use sans-serif font for random variables. Then,
L (qJointϕ ∥p) = EqJointϕ (θ,z) [ log
p(θ, z, y|x) qJointϕ (θ, z)
] , (5)
where ϕ are variational parameters. Usually the above expectation is not tractable and one uses a Monte-Carlo estimator. A single sample estimator is given by
L̂ = log p(θ, z, y|x) qJointϕ (θ, z) , (6)
where (θ, z) ∼ qJointϕ and L̂ is unbiased. Since qJointϕ need not factorize, even a single estimate requires sampling all the latent variables at each step. This is problematic when there are a large number of local latent variables. One encounters the same scaling problem when taking the gradient of the ELBO. We can estimate the gradient using any of the several available estimators [30, 20, 31, 32]; however, none of them scale if qJointϕ does not factorize [16].

heading: 4 Branch approximations for HBD
text: It is easy to see that the posterior distribution of the HBD eq. (1) takes the form
p(θ, z|y, x) = p(θ|y, x) N∏ i=1 p(zi|θ, yi, xi).
As such, it is natural to consider variational distributions that factorize the same way (see fig. 2.) In this section, we confirm this intuition—we start with any joint variational family which can have any dependence between θ and z1, · · · , zN . Then, we define a corresponding “branch” family where z1, · · · , zN are conditionally independent given θ. We show that inference using the branch family will be at least as accurate as using the joint family. We formalize the idea of branch distribution in the next definition.
Definition 1. Let qJointϕ be any variational family with parameters ϕ. We define qBranchv,w to be a corresponding branch family if, for all ϕ, there exists (v, {wi}Ni=1), such that,
qBranchv,w (θ, z) = qv(θ) N∏ i=1 qwi(zi|θ) = qJointϕ (θ) N∏ i=1 qJointϕ (zi|θ). (7)
θ
z1 z2 zN. . .
Figure 2: qBranchv,w as in eq. (7).
Given a joint distribution, a branch distribution can always be defined by choosing w and vi as the components of ϕ that influence qJointϕ (θ) and qJointϕ (zi|θ), respectively, and choosing qw(θ) and qwi(zi|θ) correspondingly. However, the choice is not unique (for instance, the parameterization can require transformations—different transformations can create different variants.)
The idea to use qBranchv,w is natural [15, 2]. However, one might question if the branch variational family is as good as the original qJointϕ ; in theorem 2, we establish that this is indeed true.
Theorem 2. Let p be a HBD, and qJointϕ (θ, z) be a joint approximation family parameterized by ϕ. Choose a corresponding branch variational family qBranchv,w (θ, z) as in definition 1. Then,
min v,w KL (qBranchv,w ∥p) ≤ min ϕ KL (qJointϕ ∥p) .
We stress that qBranchv,w is a new variational family derived from but not identical to q Joint ϕ . Theorem 2 implies that we can optimize a branched variational family qBranchv,w without compromising the quality of approximation (see appendix D for proof.) In the following corollary, we apply theorem 2 to a joint Gaussian to show that using a branch Gaussian will be equally accurate.
Corrolary 3. Let p be a HBD, and let qJointϕ (θ, z) = N ((θ, z)|µ,Σ) be a joint Gaussian approximation (with ϕ = (µ,Σ)). Choose a variational family
qBranchv,w (θ, z) = N (θ|µ0,Σ0) N∏ i=1 N (zi|µi +Aiθ, Σi)
with v = (µ0,Σ0) and wi = (µi,Σi, Ai). Then, min v,w KL (qBranchv,w ∥p) ≤ min ϕ KL (qJointϕ ∥p) .
In the above corollary, the structured family qBranchv,w is chosen such that it can represent any branched Gaussian distribution. Notice, the mean of the conditional distribution is an affine function of θ. This affine relationship appears naturally when you factorize the joint Gaussian over (θ, z). For more details see appendix F.

heading: 4.1 Subsampling in branch distributions
text: In this section, we show that if p is an HBD and qBranchv,w is as in definition 1, we can estimate ELBO using local observations and scale better. Consider the ELBO
L (qBranchv,w ∥p) = E qBranchv,w (θ,z)
[ log
p (θ, z, y|x) qBranchv,w (θ, z)
]
= E qv(θ)
[ log p (θ)
qv (θ)
] + N∑ i=1 E qv(θ) E qwi (zi|θ) [ log p (zi, yi|θ, xi) qwi (zi|θ) ] . (8)
Without assuming special structure (e.g. conjugacy) the above expectations will not be available in closed form. To estimate the ELBO, let (θ, {zi}Ni=1) ∼ qBranchv,w . Then, an unbiased estimator is
L̂ = log p (θ) qv (θ) + N∑ i=1 [ log p (zi, yi|θ, xi) qwi (zi|θ) ] . (9)
Unlike the joint estimator of eq. (6), one can subsample the terms in eq. (9) to create a new unbiased estimator. Let B be randomly selected minibatch of indices from {1, 2, . . . , N}. Then,
L̂ = log p (θ) qv (θ) + N |B| ∑ i∈B [ log p (zi, yi|θ, xi) qwi (zi|θ) ] , (10)
is another unbiased estimator of ELBO. In figs. 3b and 3c, we present the complete pseudocodes for ELBO estimation with and without subsampling in branch distributions.
Unsurprisingly, the same summation structure appears for gradients estimators of branch ELBO, allowing for efficient gradient estimation. With subsampled evaluation and training, branch distributions are immensely computationally efficient—in our experiments, we scale to models with 103 times more latent variables by switching to branch approximations (see fig. 6).
While branch distributions are immensely more scalable than joint approximations, the number of parameters still scales as O(N). In the next section, we demonstrate that for symmetric HBDs, we can share parameters for the local conditionals (amortize) to allow further scalability.

heading: 5 Amortized branch approximations
text: In this section, we discuss how one can amortize the local conditionals of a branch approximations when the target HBD is symmetric (see eq. (2).) We first formally introduce the amortized branch distributions in the next definition and then justify the amortization for symmetric HBD.
Definition 4. Let qJointϕ be a joint approximation and let qBranchv,w be as in definition 1. Suppose netu(xi, yi) is some parameterized map (with parameters u) from local observations (xi, yi) to space of wi. Then,
qAmortv,u (θ, z) = qv(θ) N∏ i=1 qnetu(xi,yi)(zi|θ) (11)
is a corresponding amortized branch distribution.
The idea to amortize is natural once you examine the optimization for symmetric HBDs. Consider the optimization for objective in eq. (8).
max v,w L (qBranchv,w ∥p) = max v,w
[ E
qv(θ)
[ log p (θ)
qv (θ)
] + N∑ i=1 E qv(θ) E qwi (zi|θ) [ log p (zi, yi|θ, xi) qwi (zi|θ) ]]
= max v
[ E
qv(θ)
[ log p (θ)
qv (θ)
] + N∑ i=1 max wi E qv(θ) E qwi (zi|θ) [ log p (zi, yi|θ, xi) qwi (zi|θ) ]] .
The crucial observation in the above equation is that for any given v, the optimal solution of inner optimization depends only on local data points (xi, yi), i.e.,
w∗i = argmax wi E qv(θ) E qwi (zi|θ)
[ log
p (zi, yi|θ, xi) qwi (zi|θ)
] . (12)
Now, notice that if p and qBranchv,w have symmetric conditionals, then, for each i, we solve the same optimization over wi, just with different parameters yi and xi. Thus, one could replace the optimization over wi with an optimization over a parameterized function from (xi, yi) to the space of wi. Formally, when the network netu is sufficiently capable, we make the following claim.
Claim 5. Let p be a symmetric HBD and let qJointϕ be some joint approximation. Let qAmortv,u be as in definition 1. Suppose that for all v, there exists a u, such that,
netu(xi, yi) = argmax wi E qv(θ) E qwi (zi|θ)
[ log
p (zi, yi|θ, xi) qwi (zi|θ)
] . (13)
Then,
min v,u KL (qAmortv,u ∥p) ≤ min ϕ KL (qJointϕ ∥p) (14)
Note, we only amortize the conditional distribution qwi (zi|θ) and leave qv (θ) unchanged. In practice, of course, we do not have perfect amortization functions. The quality of the amortization depends on our ability to parameterize and optimize a powerful neural network. In other words, we make the following approximation
netu(xi, yi) ≈ argmax wi E qv(θ) E qwi (zi|θ)
[ log
p (zi, yi|θ, xi) qwi (zi|θ)
] . (15)
In our experiments, we found the amortized approaches work well even with moderately sized networks. Due to parameter sharing, the amortized approaches converge much faster than other alternatives (especially, true for larger models; see fig. 6 and table 2).

heading: 6 Amortized branch approximations for i.i.d. observations
text: In the previous section, we discussed how we could amortize branch approximations for symmetric HBDs. However, in some applications, the construction of amortization network netu is not as straightforward. Consider the case when we have a varying number of local i.i.d observations for each local latent variable. In this section, we highlight the problem with naive amortization for locally i.i.d HBDs, and present a simple solution to alleviate them.
Mathematically, for locally i.i.d HBDs we have yi = {yij}nij=1 and xi = {xij} ni j=1, such that the conditional over yi factorizes as
p (yi|xi, zi, θ) = ni∏ j=1 p (yij |xij , zi, θ) .
Now, if xi and yi are directly input to the amortization network netu, the input size to the network would change for different i (notice we have ni observations for ith local variable.) Another problem is that the optimal variational parameters are invariant to the order in which the i.i.d. observations are presented. For instance, consider two data points: (xi, yi) = [(xi1, yi1), . . . , (xini , yini)] and (x′i, y ′ i) = [(xini , yini), . . . , (xi1, yi1)]. A naive amortization scheme will evaluate very different conditionals for these two data points because netu(xi, yi) and netu(x′i, y ′ i) will be different.
To deal with both issues: variable length input and permutation invariance, we suggest learning a “feature network” and "pooling function" based amortization network; this is reminiscent of "deep sets" [39] albeit here intended not just to enforce permutation invariance but also to deal with inputs of different sizes. Firstly, a feature network feat_net takes each (xij , yij) pair and returns a vector of features ej . Secondly, a pooling function pool takes the collection {ej}nij=1 and achieves the two aims. First, it collapses ni feature vectors into a single fixed-sized feature e (with the same dimensions as ej). Second, pooling is invariant by construction
to the order of observations (for example, pooling function would take a dimension-wise mean or sum across j.) Finally, this pooled feature vector e is input to another network param_net that returns the final parameters wi. The pseudocode for a netu with feature networks is available in fig. 4. In table 4, in appendix, we summarize the applicability of proposed variational methods to different HBD variants.
7 Experiments
v,w , and
qAmortv,u methods.
We conduct experiments on a synthetic and a real-world problem. For each, we consider three inference methods: using a joint distribution, using a branch variational approximation, and using our amortized approach. For each method, we consider three variational approximations a completely diagonal Gaussian, a block-diagonal Gaussian (with blocks for θ and z) and a dense Gaussian; see fig. 5 for a visual. (For each choice of a joint distribution, the corresponding qBranchv,w is used for the branch variational approximation and the corresponding qAmortv,u for the amortized approach; see table 1 for details.)
We use reparameterized gradients [20] and optimize with Adam [19] (see appendix G for complete experimental details.) In appendix A, we discuss some of the observations we had during experimen-
tation involving parameterization, feat-net architecture, batch-size selection, initialization, and gradient estimators.

heading: 7.1 Synthetic problem
text: The aim of the synthetic experiment is to work with models where we have access to closed-form posterior (and the marginal likelihood.) If the variational family contains the posterior, one can expect the methods to perform close to ideal (provided we can optimize well.)
For our experiments, we use the following hierarchical regression model
p(θ, z, y|x) = N (θ|0, I) N∏ i=1 N (zi|θ, I) ni∏ j=1 N (yij |x⊤ijzi, 1), (16)
where N denotes a Gaussian distribution, and I is an identity matrix (see appendix G.2 for posterior and the marginal closed-form expressions.)
To demonstrate the performance of our methods, we experiment with three different problems scale (correspond to three different models with N = 10, 1K, and 100K.) Synthetic data is created using forward sampling for each of the scale variants independently. We avoid any test data and metrics for the synthetic problem as the log-marginal is known in closed form. The inference results are present in fig. 7 (in appendix.) In all cases, amortized distributions perform favorably when compared to branch and joint distributions.

heading: 7.2 MovieLens
text: Next, we test our method on the MovieLens25M [14], a dataset of 25 million movie ratings for over 62,000 movies, rated by 162,000 users, along with a set of features (tag relevance scores [38]) for each movie.
Table 3: Details for different scales of the MovieLens 25M problem.
Scale # of ratings # of users # ratings# of users Small 2.5K 16 156.3 Moderate 180K 1600 112.5 Large 18M 159978 112.5
Purely, to make experiments more efficient on GPU hardware, we pre-process the data to drop users with more than 1,000 ratings—leaving around 20M ratings. Also, for the sake of efficiency, we PCA the movie features to reduce their dimensionality to 10. We used a train-test split such that, for each user, one-tenth of the ratings are in the test set. This gives us ≈ 18M ratings for training (and ≈ 2M ratings for testing.) We use the hierarchical model
p(θ, z, y|x) = N (θ|0, I) N∏ i=1 N (zi|µ(θ),Σ(θ)) ni∏ j=1 B(yij |sigmoid(x⊤ijzi)), (17)
v,w , and q
Amort
v,u method.) Left to right:
small, moderate, and large scale of the MovieLens25M problem. For clarity, we plot the exponential moving average of the training ELBO trace with a smoothing value of 0.001. These traces correspond to the values reported in table 2.
where θ represents distribution over user preferences; for instance, θ might represent that users who like action films tend to also like thrillers but tend to dislike musicals; zi determine the user specific preference; xij are the features of the jth movie rated by user i; yij is the binary movie ratings; ni is the number of movies rated by user i, and B denotes a Bernoulli distribution. Here, µ and Σ are functions of θ, such that, for θ = [θµ, θΣ], we have
µ(θ) = θµ, and Σ(θ) = tril(θΣ)⊤tril(θΣ)
where tril is a function that transforms an unconstrained vector into a lower-triangular positive Cholesky factor. As movie features xij ∈ R10, we have θµ ∈ R10, θΣ ∈ R55, and zi ∈ R10. Note that as Σ depends on θ, and the likelihood is Bernoulli, the model is non-conjugate.
For inference, we use the methods as described in table 1. Note, we hold the amortization network architecture constant across the scales–the number of parameters remains fixed for qAmortv,u (for all Gaussian variants) while the number of parameters scale as O(N) for qBranchv,w (see appendix F for more details.)
In fig. 6, we plot the training time ELBO trace, and in table 2, we present the final training ELBO and test likelihood values. We approximate the test likelihood p(ytest|xtest, x, y) with Eq [p(ytest|xtest, q)], and draw a fresh batch of 10,000 samples to approximate the expectation (see appendix G for complete details.) For the smaller model, the amortized and branch approaches perform similar to the joint approach for all three variational approximations; this supports theorem 2 and claim 5. For the moderate size model, the branch and amortize approaches are very comparable to each other, while joint approaches fail to scale. For the large model (18M ratings), amortized approaches are significantly better than branch methods. We conjecture this is because parameter sharing in amortized approaches improves convergence for models where batch size is smaller compared to total iterates–true only for large model. Interestingly, the performance of amortized Dense and Block Gaussian approximation is very similar in the large and moderate setting (see qAmortv,u results in tables 2, 7 and 8.) We conjecture this is because the posterior over the global parameters is
very concentrated for this problem. As θ behaves like a single fixed value, Block Gaussian performs just as good as the Dense approach (see appendix B for more discussion.)

heading: 8 Discussion
text: In this paper, we present structured amortized variational inference scheme that can scale to large hierarchical models without losing benefits of joint approximations. Such models are ubiquitous in social sciences, ecology, epidemiology, and other fields. Our ideas can not only inspire further research in inference but also provide a formidable baseline for applications.
Acknowledgements and Disclosure of Funding
This material is based upon work supported in part by the National Science Foundation under Grant No. 1908577.

Review summaries:
This paper makes an attempt to address the problem of large-scale inference for hierarchical probabilistic graphical models that contain both global and local latent variables by proposing to use Amortized VI instead of VI in [1] for more efficient and scalable inference. The authors first define a general class of hierarchical models which they call Hierarchical Branch Distributions (HBDs). They then propose a variational family matching this class of models for performing amortized inference. The authors then perform a series of experiments on synthetic datasets and the MoveLens dataset where they evaluate different inference models in terms of the final ELBO and convergence.  

This work proposes a Gaussian naive bayes approximate posterior for a 2-level hierarchical naive bayes model. Various parameterizations for the variational posterior are explored, including a fully coupled joint posterior, a factored branch posterior, and an amortized version of the latter. Results and theory indicate that having an approximate posterior that makes the same independence assumptions as the true posterior does not hurt, and that parameter sharing improves performance at scale.

An approach vor variational inference in large-scale hierarchical models is presented that is more efficient than the base lines. This is achieved using shared parameters (amortization) and a feature pooling network. Experiments confirm the efficiency of the approach on synthetic data as well as the MovieLens dataset.

The paper deals with scalable approximate inference for hierarchical models and proposes a framework for hierarchical approximations that exhibit the same factorization as the exact posterior of the hierarchical model. This efficiently enables subsampling the data when evaluating the evidence lower bound and its gradients. The authors further propose a pooling-based strategy for amortized inference to avoid the usual linear scaling in the number of approximation parameters.  The paper is concluded with a set of numerical experiments with both synthetic and real data demonstrating the merits of the proposed framework.