Year: 2022

Abstract:
We are concerned with the problem of distributional prediction with incomplete features: The goal is to estimate the distribution of target variables given feature vectors with some of the elements missing. A typical approach to this problem is to perform missing-value imputation and regression, simultaneously or sequentially, which we call the generative approach. Another approach is to perform regression after appropriately encoding missing values into the feature, which we call the discriminative approach. In comparison, the generative approach is more robust to the feature corruption while the discriminative approach is more favorable to maximize the performance of prediction. In this study, we propose a hybrid method to take the best of both worlds. Our method utilizes the black-box variational inference framework so that it can be applied to a wide variety of modern machine learning models, including the variational autoencoders. We also confirmed the effectiveness of the proposed method empirically.

heading: 1 INTRODUCTION
text: We address the problem of prediction with incomplete data, which often arises in real-life data due to the lack of data collecting resources and/or privacy concerns. For example, consider analyzing electronic health records. The input variables of the prediction are patients’ demographic characteristics and history of medical measurement data obtained with various instruments and the target variable is the survival time of patients. Thus, some of these features are not available depending on patients due to, e.g., non-standardized medical equipment, legal regulation and privacy concerns. More formally, each input is an incomplete set of features represented with a tuple u = (x̃,m) such that x̃ ∈ Rdx is the corrupted version of the complete input x ∈ Rdx and m ∈ {0, 1}dx indicates the missing entries of x̃: For all j ∈ [dx], mj = 1 if and only if the j-th feature is missing. In particular, missing entries are filled with zero and non-missing entries are identical to the one in the complete input, x̃j = (1−mj)xj . Given the incomplete features u, we want to predict the outcome y ∈ Rdy , which is modeled with a predictive distribution y ∼ p(y|u, θ). Here, the parameter θ is unknown and should be learned from data.
One straightforward way to deal with such incompletion is to incorporate generative models: The missing values in x̃ are imputed with some generative model and the target y is predicted based on the result of the imputation. Sometimes the imputation is done implicitly or simultaneously with the prediction. An advantage of this method lies in that the meaning of m is inherently incorporated. Without such context, algorithms must learn by themselves that m is actually indicating missingness (if that is relevant to the prediction task), which incurs extra generalization errors. In other words, appropriate generative models make the prediction robust to the feature incompletion.
On the other hand, the discriminative models that directly learn the conditional density pθ(y|u) sometimes exhibit superior performance (Kuhn et al., 2013) due to their alignment in terms of the objective. This is especially the case when the feature incompletion is not a dominant factor (Lasserre et al., 2006). However, it is not as straightforward as generative models to incorporate the context of missingness into discriminative models and therefore their prediction may be less robust to the feature incompletion.
To best of our knowledge, there have been no general solution to this dilemma of the generative and the discriminative approaches, i.e., generic method of discriminative inference with incomplete features for generative models, which we abbreviate as DIG. In particular, although there have been similar efforts made in the literature, the previous methods are limited in terms of the applicable models. For example, several DIG methods have been proposed for exponential families (Ghahramani & Jordan, 1994; Smola et al., 2005) and Gaussian processes (Pacheco et al., 2014), but they cannot handle modern machine-learning architectures such as deep neural networks.
This motivates us to seek for general and widely applicable frameworks that solve the dilemma. To this end, we propose a new variational approximation and a new variance reduction technique, which enable us to perform DIG within the black-box variational inference (VI) framework (Jordan et al., 1999; Zhang et al., 2018). We have also empirically confirmed the effectiveness of the proposed method through numerical experiments, employing the variational autoencoders (VAEs) (Kingma & Welling, 2013) as the base generative model.
The rest of the paper is organized as follows. In Section 2, we introduce the mathematical notation and the technical challenges in DIG with flexible models. Then, in Section 3, we present our method adopting VI in the DIG setting. In Section 4, we show the experimental results that confirm the effectiveness of the proposed method. In Section 5, we review and discuss the related work. Finally, in Section 6, we conclude the paper with some remarks and future perspectives. All the proofs for theoretical statements are given in the appendix.

heading: 2 PRELIMINARY
text: In this section, we introduce the mathematical formulation of the DIG method and then review some background of its technical challenges.

heading: 2.1 DIG FORMULATION
text: Let i ∈ [n] denote the index of instances, where n is the number of data and [n] denotes the set of integers {1, 2, . . . , n}. Let y[n] := {yi ∈ Rdy}i∈[n] and x[n] := {xi ∈ Rdx}i∈[n], be the sets of the target and the complete feature vectors, where dy and dx are the dimensionalities of vectors. We assume we can only observe the set of incomplete feature vectors u[n] := {ui := (x̃i,mi)}i∈[n] instead of x[n], where x̃i ∈ Rdx is the corrupted version of xi and mi ∈ {0, 1}dx is the indicator of corruption (‘mask vector’) for all i ∈ [n]. Each mask vector mi indicates which entry of x̃i is corrupted; We have x̃ij = (1−mij)xij for all i ∈ [n] and all j ∈ [dx].
Generative model. Suppose that the generative process of the observables (y[n],u[n]) is modeled with a parametrized family of instance-wise identical distributions. That is, the joint probability density of (y[n],u[n]) is factored as
p(y[n],u[n]|θ) = ∏ i∈[n] p(yi,ui|θ), (1)
where the function p(yi,ui|θ) on the RHS denotes the common density function across instances i ∈ [n] indexed with any parametric models θ (e.g., by neural networks).
Discriminative objective. The goal of DIG is to maximize the conditional likelihood of the corrupted data,
L[n]y|u(θ) := ln p(y [n]|u[n], θ) = ∑ i∈[n] ln p(yi|ui, θ). (2)
Now, let Liy|u(θ) := ln p(y i|ui, θ) denote the i-th summand. Following Bayes’ rule, each Liy|u(θ) is computed with a difference of log likelihoods
Liy|u(θ) = L i y,u(θ)− Liu(θ), (3) whereLiy,u(θ) := ln p(yi,ui|θ) andLiu(θ) := ln ∫ p(yi,ui|θ) dyi. Since we employ gradient-based optimization methods, it suffices to consider (approximate) evaluation of the value and the gradient
of the instance-wise objectives Liy|u(θ) separately. Thus, in the following, we limit our focus to the separate objective and omit the index i if there is no risk of ambiguity.

heading: 2.2 DIG WITH LATENT VARIABLE MODELS (DIG-LVM)
text: We give further details of DIG with a specific class of generative models, namely the latent variable models (LVMs). LVM is one of the most common and expressive generative models. In the context of DIG, the generative density function is given by
p(y,u|θ) = ∫ p(y,u, z|θ) dz, (4)
where z is the instance-wise latent variable, which is expected to capture some higher-order information of the observables. For example, the variational autoencoders (VAEs) (Kingma & Welling, 2013), the generative adversarial networks (GANs) (Goodfellow et al., 2014) and the normalizing flow models (Rezende & Mohamed, 2015) are instances of LVM.
Modeling feature-corruption process with LVM. One advantage of LVM is that complex generative processes can be modeled with relatively simple joint density p(y,u, z|θ). This is especially beneficial in modeling missing data: The full joint density p(y,u, z|θ) can be easily designed to reflect one’s belief on the corruption process such as MCAR and MNAR (Rubin, 1976) through the factorization of the density (e.g., see Collier et al. (2020)), while maintaining the expressibility of LVM. For example, the MCAR process can be modeled with
p(y,u, z|θ) = p(m) · p(z) · p(y|z, θ) · dx∏ j=1 {p(x̃j |z, θ)}1−mj , (5)
where p(m) is the true marginal density function of m, which is unknown but cancels out in the DIG objective (Equation 3), and p(z) is an arbitrary fixed prior of the latent variable (such as Gaussian). On the other hand, the MNAR process can be modeled with
p(y,u, z|θ) = p(m) · p(z) · p(y|z,m, θ) · dx∏ j=1 {p(x̃j |z,m, θ)}1−mj . (6)
Note that these are just examples and it is possible to further incorporate domain knowledge on the generative process. An example of such models is the MAR process (Rubin, 1976) and its formulation is relegated to the appendix (Section G).
Under these factorizations, the log likelihoods on the RHS of Equation 3 for LVMs are computed as Ly,u(θ) = ln ∫ p(y,u, z|θ) dz, Lu(θ) = ln ∫ p(u, z|θ) dz, (7)
where both density functions in the integrals can be analytically computed under either of Equation 5 and Equation 6.
Challenge: Inference with DIG-LVM. To maximize Equation 3 with LVM, we have to evaluate the values and gradients of the likelihoods in Equation 7. If there are only positive log integrals in the objective, this can be approximately done with the variational inference (VI) framework (Jordan et al., 1999; Zhang et al., 2018) in favor of flexibility and scalability. However, with DIG, we also have a negative log integral, i.e., −Lu(θ). The optimization of such objective is relatively difficult since there have been no equivalent of VI in terms of the flexibility and scalability (see Section 5).

heading: 3 VARIATIONAL INFERENCE FOR DIG-LVM
text: Now we present the main result, the method of VI for DIG-LVM. In Section 3.1, we derive a variational approximation to the conditional likelihood given by Equation 3. Then, in Section 3.2 and 3.3, we derive a method for optimizing the variational approximation and its variant for variance reduction, respectively. Finally, in Section 3.4, we show the prediction procedure based on the optimized parameters.

heading: 3.1 VARIATIONAL LOWER/UPPER BOUNDS
text: In this section, we discuss an approximation of Ly,u(θ) and Lu(θ) in Equation 7 individually. Note that both are log integrals of the specific form
Lv(θ) := ln ∫ p(v, z|θ) dz, (8)
where v represents either (y,u) or u. Since such integrals constitute the objectiveL(θ) with positive and negative signs, our goal is to derive both upper and lower bounds on Equation 8.
Evidence Lower Bound (ELBO). The lower bound is given in the standard way (Jordan et al., 1999; Zhang et al., 2018),
Lv(θ) ≥ Lv(θ)−DKL(z|v)(φ‖θ) = Ez∼q(z|v,φ) [ ln p(v, z|θ) q(z|v, φ) ] =:LELBO(v)(θ, φ), (9)
where φ is a variational parameter of a probability density function q(z|v, φ) and DKL(z|v)(φ‖θ) :=Ez∼q(z|v,φ)[ln q(z|v,φ)p(z|v,θ) ] is the KL divergence of the parameters φ and θ, We call LELBO(v)(θ, φ) as the evidence lower bound (ELBO). Since ELBO is an expectation of a tractable function, we approximate it with Monte-Carlo sampling.
Evidence Upper Bound (EUBO). To derive an upper bound, we start with applying the χevidence upper bound (CUBO) (Dieng et al., 2017). For any real numbers α > 1, CUBO is derived as follows:
Lv(θ) ≤ Lv(θ) + (1− α−1)Dα(z|v)(θ‖ψ)
= 1
α lnEz∼q(z|v,ψ) [ p(v, z|θ) q(z|v, ψ) ]α =:LCUBO(v)(θ, ψ), (10)
where ψ is a variational parameter of a probability density function q(z|v, ψ) and Dα(z|v)(θ‖ψ) := 1α−1 ln ∫ dz pα(z|v, θ)q1−α(z|v, ψ) denotes the α-Rényi divergence of the parameters θ and ψ. Note here, unlike ELBO, CUBO is not unbiasedly approximated because of the logarithm wrapping the expectation. To address this issue, we apply another variational approximation with a divergence function Ψα(t) :=(eαt−αt−1)/α, t ∈ R: Since Ψα(t) ≥ 0 for all t ∈ R, we have
LCUBO(v)(θ, ψ) ≤ LCUBO(v)(θ, ψ) + Ψα(LCUBO(v)(θ, ψ)− f(v; ξ))
= e−αf(v;ξ)
α Ez∼q(·|v,ψ) [ p(v, z|θ) q(z|v, ψ) ]α + f(v; ξ)− 1 α =:LEUBO(v)(θ, ψ, ξ),
(11)
where ξ is a variational parameter of a real-valued function f(v; ξ). We call the right-hand side as the evidence upper bound (EUBO). Note that the expectation in EUBO is linearized and thus can be unbiasedly approximated with Monte-Carlo estimation.
Conditional Evidence Lower Bound (CELBO). Applying ELBO on v = (y,u) and EUBO on v = u, we have a conditional evidence lower bound (CELBO),
LCELBO(y|u)(θ, φ, ψ, ξ) :=LELBO(y,u)(θ, φ)− LEUBO(u)(θ, ψ, ξ). (12)
By definition, CELBO bounds the DIG objective from below, Ly|u(θ) ≥ LCELBO(y|u)(θ, φ, ψ, ξ). The inequality is tight for any generative parameter θ owing to the tightness of ELBO and EUBO,1 i.e., for all θ, there exists a tuple of variational parameters (φ, ψ, ξ) such that the gap ∆(θ, φ, ψ, ξ) :=Ly|u(θ) − LCELBO(y|u)(θ, φ, ψ, ξ) is zero. Moreover, CELBO can be unbiasedly approximated as well as ELBO and EUBO.
1ELBO is tight if q(z|v, φ) = p(z|v, θ). EUBO is tight if q(z|v, ψ) = p(z|v, θ) and f(v; ξ) = LCUBO(v)(θ, ψ) = lnLv(θ).
Algorithm 1 Variational Inference for DIG (vDIG) Input: Data (y[n],u[n]) Output: θ, φ, ψ, ξ 1: (θ, φ, ψ, ξ)← Initialize() // Any initialization methods can be used. 2: repeat 3: Draw minibatch B ⊂ [n] 4: L ← 1|B| ∑ i∈B L̂ i CELBO(y|u)(θ, φ, ψ, ξ) // L̂CELBO(y|u) is given by Equation 13. Use Equation 15
instead for variance reduction. 5: (θ, φ, ψ, ξ)← Update((θ, φ, ψ, ξ),∇L) // Any gradient-based optimization methods can be used. 6: until converge

heading: 3.2 OPTIMIZATION ALGORITHM: VDIG
text: Since CELBO can be unbiasedly approximated, we may employ stochastic gradient-based optimization to maximize it. The full stochastic objective for CELBO is given by
L̂CELBO(y|u)(θ, φ, ψ, ξ) := ln p(y,u, zφ|θ) q(zφ|y,u, φ) − 1 α
( p(u, zψ|θ)
q(zψ|u, ψ)ef(u;ξ)
)α − f(u; ξ) + 1
α , (13)
where zφ and zψ are Monte-Carlo samples drawn from q(z|y,u, φ) and q(z|u, ψ), respectively. The gradients of L̂CELBO(y|u) is taken with any standard automatic differentiation libraries, using the reparametrization trick (Kingma & Welling, 2013) or the REINFORCE trick (Williams, 1992).
Since the actual objective function is the summation of individual losses L̂CELBO(y|u) = L̂iCELBO(y|u) over all the instances, we may draw a minibatch of instances for each iteration. We call the resulting inference algorithm as vDIG (Algorithm 1).

heading: 3.3 VARIANCE REDUCTION FOR VDIG WITH SURROGATE PARAMETRIZATION (SP)
text: The boundedness of the norm of the stochastic gradient ∇L̂CELBO(y|u)(θ, φ, ψ, ξ) is crucial for the stable and fast convergence of stochastic gradient-based algorithms like Algorithm 1. However, the stochastic CELBO contains the density ratio
wθ,ψ,ξ(u, z) := p(u, z|θ)
q(z|u, ψ)ef(u;ξ) (14)
raised to the power of α, which is problematic as the ratio wθ,ψ,ξ may have large variance and so does the gradient.
The key idea is to regularize the parameter to keep the density ratio small. Note that we have wθ,ψ,ξ(·, ·) ≡ 1 whenever the objective gap is zero. That is, constraining the parameter to satisfy supu,z wθ,ψ,ξ(u, z) ≤ 1 does not lose the model expressibility if the variational approximation is tight.
The problem is, enforcing such constraint during optimization is intractable in general. We address this issue by introducing a surrogate parametrization. Define new parameters θ′ and ξ′ formally2 by
p(y,u, z|θ′) := G(wθ,ψ,ξ(u, z)) Z(θ′) p(y,u, z|θ), f(u; ξ′) := f(u; ξ)− lnZ(θ′),
where G(w) :=(1 ∨ w)−1{1 + α ln(1 ∨ w)}1/α (w ≥ 0), a ∨ b := max{a, b}, and Z(θ′) is the normalizing constant ensuring the mass preservation of the density under θ′. We refer to the mapping TSP : (θ, φ, ψ, ξ) 7→ (θ′, φ, ψ, ξ′) as the surrogate transform, and G(w) as the gain function. The surrogate transform TSP and the gain function G(w) are designed carefully to satisfy the following two properties.
First, it preserves the effective parameters of CELBO.
2Note that these parameters are conceptual objects and there is no concrete implementation for them in the final algorithm.
Definition 1 (Effective parameters). Let Ω be a set of CELBO parameters (θ, φ, ψ, ξ). Then, we define the effective parameters of Ω by Θ0(Ω) :={(θ, φ, ψ, ξ) ∈ Ω : ∆(θ, φ, ψ, ξ) = 0}, i.e., the set of parameters inducing tight variational approximation. Proposition 2 (Effective parameter preservation). For arbitrary Ω, we have Θ0(TSP(Ω)) = Θ0(Ω).
In other words, in a sense, TSP does not alter the original model. See Section H for extended discussion.
Second, it regularizes the growth of the highly stochastic term (Equation 14). Intuitively, the effect of regularization is quantified with the gain function G(w) since it represents the ratio of stochastic term before and after the transform,
wθ′,ψ,ξ′(u, z)
wθ,ψ,ξ(u, z) = G(wθ,ψ,ξ(u, z)).
Note that G(w) is no larger than one and non-increasing for w ≥ 0, thus the stochasticity is reduced through the transform. See Figure 3 in the appendix for visualization. Consequently, TSP guarantees the boundedness of the stochastic gradient. Proposition 3 (Bounded gradient). Assume there exists K > 0 such that ‖∇ ln p(y,u, zφ|θ)‖, ‖∇ ln p(u, zψ|θ)‖, ‖∇ ln q(zφ|y,u, φ)‖, ‖∇ ln q(zψ|u, ψ)‖, ‖∇f(u; ξ)‖ ≤ K. Then, ‖∇(L̂CELBO(y|u) ◦ TSP)(θ, φ, ψ, ξ)‖ ≤ 9K.
Proposition 2 and 3 justifies optimizing the objective via the transform TSP; we get a gradientnorm bound without altering the effective parameters. The full objective function after the surrogate transform, L̂CELBO-SP(y|u) := L̂CELBO(y|u) ◦ TSP, is given by
L̂CELBO-SP(y|u)(θ, φ, ψ, ξ) = ln p(y,u, zφ|θ)G(wθ,ψ,ξ(u, zφ))
q(zφ|y,u, φ)
− 1 α {wθ,ψ,ξ(u, zψ)G(wθ,ψ,ξ(u, zψ))}α − f(u; ξ) + 1 α , (15)
which we call the stochastic CELBO-SP. The key point is that computing the stochastic CELBO-SP does not require the computation of the normalizing constant Z(θ′), which is intractable in general. This is not the case with the surrogate transformation on ELBO or EUBO alone. The CELBO-SP maximization is done by simply replacing CELBO with CELBO-SP in Algorithm 1.

heading: 3.4 PREDICTION ALGORITHM
text: Given (θ, φ, ψ, ξ) trained by Algorithm 1, we want to compute the predictive distribution on new instances given their incomplete features unew := (x̃new,mnew). Since the conditional density pθ(y|u) is intractable in general, we approximate it with the Monte-Carlo method. The approximated conditional distribution is given by
p(y|u, θ̂) := 1 kpred ∑ s∈[kpred] p(y|zsψ, θ), kpred ≥ 1, (16)
where zsψ are independently drawn from qψ(z|u), s ∈ [kpred]. This procedure is justified as follows.
Proposition 4. Let p(y|u, θ̄) :=E[p(y|u, θ̂)], where the expectation is taken with respect to the Monte-Carlo sampling. Then,
DKL(y|u)(θ‖θ̄) ≤ α
α− 1 ∆(θ, φ, ψ, ξ).
In other words, if the objective gap is small, so is the approximation error of θ̄, which is the limit of the actual predictor θ̂ with kpred →∞. In the experiment, we used kpred = 512.

heading: 4 EXPERIMENTS
text: In this section, we demonstrate the effectiveness of the proposed method, i.e., Algorithm 1, through numerical experiments. We first introduce a number of methods compared in the experiments in
Section 4.1. We then present the procedures and results of three experiments designed to show i) the effectiveness of the new variational approximation, ii) the superior performance of the proposed method, and iii) the robustness of the proposed method against feature corruption.

heading: 4.1 EXPERIMENTAL SUBJECTS
text: We employ three subject algorithms and two baselines in the experiment. The main subjects are all based on VAE (Kingma & Welling, 2013), a typical example of LVM, whose basic architectures are identical to each other: See the appendix (Section D) for more details.
Generative method: VAE, VAE*. As an example of generative approach, we employ VAE and solve missing-value imputation and regression simultaneously by regarding the target as a part of the feature, x′ ← x ⊕ y, and reconstructing x′ from its corrupted version, where the entries corresponding to y is considered as missing. In particular, we adopt the formulation of Collier et al. (2020), which comes with two variants for different missing-value processes, namely missing not at random (MNAR) and missing completely at random (MCAR). The difference between MNAR and MCAR is only in the generative model (i.e., the architecture of the decoder), represented by Equation 6 and Equation 5, respectively. We denote these variants as VAE and VAE*, respectively.
To train it to be able to reconstruct y, we double the training dataset with masking (and not masking) the target y. Note that the objective function is not exactly aligned with the prediction task because of the imputation-based formulation.
Discriminative method: CVAE. As an example of discriminative approach, we employ the conditional VAE (CVAE) proposed by Kingma et al. (2014); Sohn et al. (2015). In our setting, the conditional variables consist of the corrupted feature x̃ and the missingness indicator m. In practice, these vectors are concatenated before fed into neural networks. The generative model of CVAE in this setting is given by p(y|u, θ) := ∫ p(y, z|u, θ) dz, while the decoder represents the density in the integral and the encoder is used to approximate the posterior p(z|y,u, θ). Note that the resulting model is not informed with the context of m as VAE and VAE* are,3 but the objective is (a variational approximation of) the conditional evidence, which is aligned with the prediction task.
Proposed method: DVAE, DVAE*. We call the instantiation of vDIG on VAE as the discriminative VAE (DVAE). We implement MNAR (Equation 6) and MCAR (Equation 5) variants of DVAE as with VAE, respectively denoted by DVAE and DVAE*. The difference with VAE and its variant is that we have additional encoders corresponding to ψ and ξ, while the decoder θ and the encoder φ are common. See Section D for more details. We choose α = 2 for the parameter of the Rényi divergence.
Baseline methods: Simple, MICE. As the first baseline method, we employ a simple fullyconnected feed-forward neural network (FCFNN) with the same architecture with the encoders of the above methods, except its output is considered as a distribution of the target y instead of z. Although the output distribution is restricted to a Gaussian, it requires no variational approximation and hence the conditional likelihood ln pθ(y|u) is exactly maximized. We call this Simple. The second baseline is the method of the multiple imputation by chained equations (MICE) (Azur et al., 2011), which is a classic approach in statistics to deal with missing values. In particular, we employ the Bayesian ridge regression for the missing-value imputation task and the FCFNN of the first baseline for the subsequent regression task.

heading: 4.2 EXPERIMENTAL PROCEDURE AND RESULTS
text: The organization of the experiments is three-fold. First, in Section 4.2.1, we check the effectiveness of our variational approximation techniques, namely EUBO and the surrogate parameterization. Second, in Section 4.2.2, we compared the proposed method(s) with the existing methods in terms of the predictive performance. Finally, in Section 4.2.3, we examine the robustness of these methods against the change in the missing-value ratio. See also the appendix for the details of the experimental settings, including the information of the datasets.
3Because of this, there is no MCAR variant for CVAE.

heading: 4.2.1 EFFECT OF NEW VARIATIONAL APPROXIMATION
text: Procedure. We compare the stability of optimization of DVAE with three different variational approximation, namely, ELBO+CUBO, ELBO+EUBO and ELBO+EUBO+SP. In all cases, the Monte-Carlo approximation of ELBO is used to compute the joint evidence Ly,u(θ). The marginal evidence Lu(θ) is computed with the Monte-Carlo approximation for both CUBO (albeit biased) and EUBO. Only ELBO+EUBO+SP employs the surrogate parametrization.
Result. Figure 1 shows the training processes of DVAE with the AQ-CO dataset from UCI Machine Learning Repository (see the appendix). It is seen that ELBO+EUBO+SP is most stable in the optimization. Also note that the objective of ELBO+CUBO may be negatively biased (recall the logarithm wrapping the expectation in Equation 10), but we cannot know how much.

heading: 4.2.2 PREDICTIVE PERFORMANCE
text: Procedure. We apply seven different algorithms in Section 4.1 to six different regression tasks from UCI Machine Learning Repository, summarized in Table 3 (in the appendix). We drop the entries of feature completely at random with probability p = 0.1. In each configuration, we iterate the same procedure with five different random seed values.
Result. Table 1 shows the cross entropy, i.e., − ln pθ(y|u), averaged over the test splits. Overall, both DVAE and DVAE* almost always perform best or at least comparable to the best scores within one standard deviation, indicating the effectiveness the DIG method applied VAE. In particular, it is more clear that they outperform the others when averaged over all the datasets (‘Total’ column). Moreover, DVAE* is slightly better than DVAE as expected because the way we drop the feature entries is MCAR.

heading: 4.2.3 ROBUSTNESS AGAINST FEATURE CORRUPTION
text: Procedure. We take the result of Section 4.2.2 and examine how the missing-value ratio p affects predictive performance. For each method, we calculate the difference between the average cross entropy with p = 0.1 and p = 0.5.
Result. The difference is visualized in Figure 2. The most robust algorithms are VAE and VAE*, while DVAE* and DVAE are the second and the third most robust algorithms. This matches our expectation because both VAE and VAE* are informed with the context of m and implicitly regularized through missing-value imputation task. Similarly, DVAE and DVAE* are also informed with the meaning of m, but they are not subject to the implicit regularization. In particular, DVAE* is more robust than DVAE as the MCAR model is accurate (in this setting) and more informative than the MNAR model.

heading: 5 RELATED WORK
text: In this section, we review existing variational approximations bounding the log-integrals from above. We also discuss the relationship with the existing studies on missing data handling in the appendix.
The χ-variational upper bound (CUBO) (Dieng et al., 2017) has been proposed to estimate such upper bounds, utilizing the α-Rényi divergence (equivalently, χ-divergence). However, as Pradier et al. (2019) pointed out, the resulting estimate of CUBO is biased and in some cases numerically unstable4. This phenomenon has been also confirmed in our experiment (Section 4.2.1). Another upper bound has been proposed by (Ji & Shen, 2019) with the reverse KL divergence instead of the Rényi divergence, but it is also biased and not guaranteed to be an upper bound. We also find that the parsimonious upper bound (Mattei & Frellsen, 2018) is suitable in terms of the upper-bound guarantee, while it induces a min-max form during the upper bound minimization procedure, which implies the convergence property of the algorithm could be tricky.
The proposed upper bound, EUBO, is closely related to the one studied by Kuleshov & Ermon (2017), which is designed for estimating partition functions. Although their method also suffers from high variance, they took different approach to reduce it. In particular, they adaptively reduced the learning rate of the variational parameter ψ. However, as the authors noted, this does not solve the problem enough for scaling to large datasets.

heading: 6 CONCLUSION
text: We have proposed a novel algorithm to perform discriminative training with incomplete features for generative models, which is derived on the basis of a newly introduced variational approximation and parameter transformation. The effectiveness of the proposed method has been confirmed in terms of the stability of the new upper bound and the predictive performance and robustness of the resulting algorithm.
Possible directions of future work include the application of the surrogate transform technique to other contexts than feature incompletion.
4Another work (Lopez et al., 2020) reported that there is a case CUBO works just fine, so its seems a problem dependent phenomenon.

heading: A PROOFS
text: A.1 PROOF OF PROPOSITION 2
Proof. Θ0(TSP(Ω)) ⊃ Θ0(Ω) is trivial as ∆(θ, φ, ψ, ξ) = 0 implies wθ,ψ,ξ(·, ·) ≡ 1, which implies (θ, φ, ψ, ξ) is a fixed point of TSP, hence (θ, φ, ψ, ξ) ∈ TSP(Ω). The other direction is shown as follows. Assume ∆(θ′, φ, ψ, ξ′) = 0, where (θ′, φ, ψ, ξ′) = TSP(θ, φ, ψ, ξ) and (θ, φ, ψ, ξ) ∈ Ω. Then, we have
1 = wθ′,ψ,ξ′(u, z) = wθ,ψ,ξ(u, z)G(wθ,ψ,ξ(u, z))
for all u and z. This implies wθ,ψ,ξ(·, ·) ≡ 1 as wG(w) = 1 ⇔ w = 1. Therefore, (θ, φ, ψ, ξ) is a fixed point of TSP and (θ′, φ, ψ, ξ′) ∈ Ω.
A.2 PROOF OF PROPOSITION 3
Proof. Let g :=∇(L̂CELBO(y|u) ◦ TSP)(θ, φ, ψ, ξ) and observe
g = ∇ ln p(y,u, zφ|θ)G(w) q(zφ|y,u, φ) − ∇ α (wG(w)) α −∇f(u; ξ)
= ∇ ln p(y,u, zφ|θ) q(zφ|y,u, φ) −∇f(u; ξ) +∇H(lnw)
= ∇ ln p(y,u, zφ|θ) q(zφ|y,u, φ) −∇f(u; ξ) +H ′(lnw)∇ lnw
where w :=wθ,ψ,ξ(u, z), H(t) := lnG(et) − (etG(et)) α /α, and H ′(t) denotes the derivative of H(t). Thus, we have
‖g‖ ≤ |H ′(lnw)| ‖∇ lnw‖+ 3K ≤ 9K
as |H ′(t)| ≤ 2 and ‖∇ lnw‖ ≤ 3K.
A.3 PROOF OF PROPOSITION 4
Proof. According to the information processing inequality, we have DKL(y|u)(θ‖θ̄) ≤ DKL(z|u)(θ‖ψ). Moreover, by the construction of LMEUBO, we have (1 − α−1)Dα(z|u)(θ‖ψ) ≤ ∆(θ, φ, ψ, ξ). The desired result is seen by combining these two inequalities with the fact that the α-Rényi divergence dominates the KL divergence for all α > 1.

heading: B EXPERIMENTAL CONFIGURATIONS
text: Computing infrastructure. The computing infrastructure used in the experiments is summarized in Table 2.
Preprocessing. All datasets used in the experiment are standardized before fed into algorithms so that the variables of the feature and the target have zero mean and unit empirical variance on the training split. More precisely, we apply the following transformation
x̃i ← (1−mi) (x̃i − µ̂x) σ̂x, yi ← (yi − µ̂y) σ̂y, i ∈ [n],
where µ̂x,j := ∑n i=1(1−m i j)x i j∑n
i=1(1−mij) and σ̂ 2x,j :=
∑n i=1(1−m i j)(x i j−µ̂x,j)
2∑n i=1(1−mij)
for j ∈ [dx], and µ̂y := 1 n ∑n i=1 y i and σ̂ 2y := 1 n ∑n i=1(y
i − µ̂y) 2. Here, and denotes the elementwise multiplication and division of vectors, respectively. Note that we refill the missing entries of x̃i with zero after the affine transformation, which is corresponding to filling them with the means in the original dataset. In the prediction step, the input feature x̃ undergoes the same transformation, x̃i ← (1−m) (x̃− µ̂x) σ̂y , and the predictive distribution p(y|u, θ) for the standardized target is made based on it. Then, the distribution is transformed inversely to predict the non-standardized target, p(y|u, θ)← ( ∏dy j=1 σ̂ −1 j,y ) · p((y − µ̂y) σ̂y|u, θ).
Validation split. We randomly hold 20% of the training split out for validation. Algorithms receive as the input only the remaining 80%.
Initialization. All the parameters involved in the experiments are contained in the Linear module of Pytorch and initialized with its default initialization method, i.e., independently subject to the uniform distribution on the interval [− 1√
k , 1√ k ], where k is the number of input variables for the
respective modules.
First-order optimizer. In all the experiments, we employ the AdaBelief optimizer (Zhuang et al., 2020) with its default setting,5 which is recently proposed as an improved alternative for the Adam optimizer (Kingma & Ba, 2014). The size of minibatch is always taken to be 512. In particular, for sampling-based methods (i.e., VAE, VAE*, CVAE, DVAE, DVAE*), the stochastic gradient is computed with the reparametrization trick.
Number of Iterations. We iterate the loop 2000 times for each configuration, i.e., 2000 × 512 samples are seen in one run irrespective to the data size.
Evaluation. In the above iteration, we save the models at 20 predefined points t ∈ {0, 1, 2, 3, 5, 8, 13, 19, 29, 43, 64, 94, 138, 203, 298, 436, 638, 934, 1366, 1998}, which are (approximately) equally spaced in the log scale. Finally, the best model among these 20 checkpoints is chosen with respect to the cross entropy on the validation split and evaluated on the test split. If the subject algorithm is based on LVM, then we sample kpred = 512 instances from the posterior distribution q(z|u, ψ) and make prediction with Equation 16.

heading: C DATASETS
text: All the datasets used in the experiments are taken from UCI Machine Learning Repository. AQ-CO, AQ-NMHC and AQ-NOx are taken from the AirQuality dataset (De Vito et al., 2008) and corresponding to different targets denoted by their suffixes. YearPred is a part of the YearPredicitonMSD dataset (Bertin-Mahieux, 2011), where, to accommodates the fast iteration of the experiments, the number of records is restricted to 10,000 by separate random sampling for training (8000 records) and test (2000 records) splits. Boston and Diabetes are respectively taken from the dataset of the same names. The basic statistics are summarized in Table 3.
5See https://github.com/juntang-zhuang/Adabelief-Optimizer/tree/update 0.2.0

heading: D METHODS COMPARED
text: In this section, we show the detailed architecture of the methods compared in the experiments.
Both the encoder and the decoder are fully-connected feed-forward neural networks (FCFNN) with single hidden layer of width 256 and ReLU activation that take a vector input t ∈ Rk, k ≥ 1, and emit a diagonal Gaussian distribution depending on t. To facilitate the specification of the architectures, we first define some basic building blocks.
Let every distinct occurrence of Lineard(·) denote a distinct linear layer with the output dimension d ≥ 1 with its own parameters (W,b) such that Lineard(t) = W t+b with a weight matrix W ∈ Rd×∗ and an offset vector b ∈ Rd for t ∈ R∗. The Gaussian layer is defined by the composition of the Gaussian density function and the linear layers,
Gaussd(t|v) := Nd[t |Lineard(v), + ln(1 + Lineard(v))]
where := 10−8 is a small constant and Nd[t|v, s] denotes the d-product of the probability density functions of the Gaussian distributions with mean vj and standard deviation sj evaluated at tj , j ∈ [d], for t,v ∈ Rd and s ∈ Rd>0. We also denote the ReLU activation function by ReLU(t) := max {0, t}. Finally, we define the composition of a Gaussian layer and a Linear layer with k hidden neurons and ReLU activation function by
Gauss-NNd,k(t|v) := Gaussd(t|(ReLU ◦ Lineark)(v)).
Now, we are ready to describe the architecture specifications.
Simple. The objective of the Simple baseline is given by
− ln p(y|u, θ) := − lnGauss-NNdy,256(y|x̃⊕m),
where θ denotes all the parameter implicitly involved in the RHS.
MICE. MICE is implemented with IterativeImputer from scikit-learn (ver. 0.24.2) with its default argument except with the burn-in period changed from 10 to 40. After the imputation, only the final imputation result is passed to the Simple baseline to reduce the computation time.
VAE, VAE*. The objective of VAE is given slightly modifying Collier et al. (2020). The difference is that the mask vector m is modeled as a conditional variable rather than a random variable. The reason of this modification is because we are not interested in the likelihood of m and conditioning on such variables simplifies the resulting model. Namely,
L̂VAE(θ, φ) := − ln p(y,u|z, θ)− lnN10[z; 0, 1] + ln q(z|y,u, φ),
where
p(y,u|z, θ) := p(m) · Gaussdy (y|v) · ∏
j∈[dx]:mj=0
Gauss1(x̃j |v),
q(z|y,u, φ) := Gauss-NN10,256(z|y ⊕ x̃⊕m),
v := (ReLU ◦ Linear256)(z ⊕m) and z ∼ q(z|y, x̃,m, φ). For VAE*, replace v with v∗ := (ReLU ◦ Linear256)(z). Here, p(m) is an arbitrary density function of m, which is constant with respect to the parameters (θ, φ) and ignored in the gradient computation.
CVAE. The objective of CVAE is given according to Sohn et al. (2015),
L̂CVAE(θ, φ) := − ln p(y|z,u, θ)− lnN10[z; 0, 1] + ln q(z|y,u, φ),
where
p(y|z,u, θ) := Gauss-NNdy,256(y|z⊕ x̃⊕m), q(z|y,u, φ) := Gauss-NN10,256(z|y ⊕ x̃⊕m),
and z ∼ q(z|y,u, φ).
DVAE, DVAE*. The objective of DVAE is given by Equation 13, where p(y,u|z, θ) := p(m) · Gaussdy (y|v) · ∏
j∈[dx]:mj=0
Gauss1(x̃j |v),
q(z|y,u, φ) := Gauss10(z|g256(y,0dy , x̃,m;ω)), q(z|u, ψ) := Gauss10(z|g256(0dy ,1dy , x̃,m;ω)), f(u; ξ) := ln p(m) + Linear1(g256(0dy ,1dy , x̃,m;ω)),
gk(t (1), t(2), t(3), t(4);ω) := (ReLU ◦Lineark)(t(1)⊕ t(2)⊕ t(3)⊕ t(4)) for t(1), t(2) ∈ Rdy and t(3), t(4) ∈ Rdx , ω denotes the shared parameter of φ, ψ, ξ, and v := (ReLU◦Linear256)(z⊕m). For DVAE*, replace v with v∗ := (ReLU ◦ Linear256)(z). Here, p(m) is an arbitrary density function of m, which cancels out in the final objective of CELBO or CELBO-SP.

heading: E ADDITIONAL RESULTS
text: We show in Table 4 the results of the same experiment as in Section 4.2.2 except with the increased missing-value ratio p = 0.5.

heading: F ADDITIONAL DISCUSSION ON RELATED WORK: MISSING-DATA HANDLING
text: Arguably the most classic approach to the missing value problem is the two-step approach, also known as the imputation method (Enders, 2010). This category includes traditional listwise or pairwise deletion methods, single imputation methods and multiple imputation methods. The key feature of this approach is that one processes the incomplete features to get estimates of complete ones x̂ ≈ x in the first step and then performs prediction based on x̂. The score-based methods (Jaakkola & Haussler, 1999; Holub et al., 2005; Perina et al., 2009; Śmieja et al., 2018) can be considered as a generalization of the two-step approach developed in the machine learning literature. With a score-based method, one first learns the generative distribution pθ(x), which is used to extract some information called score, s. The score is then fed to predictive models pθ(y|s). Several drawbacks stem from the two-step nature of these methods. In case of the imputation method, the predictor loses the information whether each element of x̂ is original or imputed. This makes it difficult to estimate the uncertainty resulted from the feature corruption. Moreover, even though the uncertainty
problem can be addressed with the score-based method, learning good score representations requires solving optimization problems possibly irrelevant to the original prediction problem and thus it may compromise the performance.
Another category of missing feature handling fully utilizes generative models. One of such method is referred to as the full information maximum likelihood methods (see also Secion 4, Enders (2010)), where the joint complete-data distribution pθ(y,x) is explicitly modeled and marginalized over the corrupted elements of features to obtain the objective ln pθ(y,u). This approach may suffer from unnecessarily performance degradation as in the score-based approach, since the objective contains the generative term of u, ln pθ(y,u) = ln pθ(y|u) + ln pθ(u). The third approach is the discriminative approach. Specifically, tree-based models such as gradient boosting trees are able to naturally handle incomplete features (Twala et al., 2008). Moreover, it is also recommended in (Kuhn et al., 2013) to encode missingness as a distinct feature, i.e., treat the concatenation of the corrupted feature and the mask vectors, x̃ ⊕m, as the input to the predictive models. The advantage of this approach is that the resulting objective function is coherent with the goal of predictive risk minimization, i.e., there is no generative term unlike in the two-step and generative approach. However, there is no trivial way to inform the model that m actually indicates missing features. Therefore it may take extra samples to learn the meaning of m by itself.
Finally, the DIG approach can be thought of as a hybrid of the generative and discriminative approaches. With DIG, the data is modeled with joint distribution pθ(y,u), but the learning objective is the conditional evidence ln pθ(y|u), which is computed from Bayes’ rule. Therefore, it naturally incorporates the information of missingness and is directly trained to maximize the predictive performance. The first application of the DIG strategy in the context of incomplete feature is Ghahramani & Jordan (1994), which was followed by Smola et al. (2005) with a kernel-based generalization. It is crucial in their results that the complete-data model pθ(y,x) is a exponential family so that the objective function is optimized with the EM algorithm (Dempster et al., 1977). A tractable approximation algorithm for the case of Gaussian processes is derived by Pacheco et al. (2014) with the combination of the variational inference (Jordan et al., 1999) and the expectation propagation (Minka, 2001) framework. As opposed to their method, our focus is on black-box variational inference algorithm applicable to a variety of models.
See Table 5 for the summary of the comparison.

heading: G MODELING THE MAR PROCESS
text: In this section, we give a method of modeling the third type of the feature-corruption processes, the MAR process. Let O ⊂ [dx] be an index set on which the feature is always observed (i.e., mj = 0 for all j ∈ O) and define x̃O := {x̃j : j ∈ O}. Then, the MAR process is modeled with
p(y,u, z|θ) = p(x̃O,m) · p(z) · p(y|z, x̃O, θ) · ∏
j∈[dx]\O
{p(x̃j |z, x̃O, θ)}1−mj , (17)
where p(x̃O,m) will cancel out as well as p(m) in the MNAR and MCAR processes.

heading: H ADDITIONAL JUSTIFICATION AND LIMITATION OF THE SURROGATE TRANSFORM
text: We present another justification of the CELBO-SP maximization (i.e., optimization of Equation 15) as an alternative to the CELBO maximization (i.e., optimization of Equation 13). Recall that
CELBO-SP is derived by the surrogate transform TSP and Proposition 2 shows a desirable property of TSP viewing it as an operator acting on the parameters. In this section, we show another desirable property of TSP viewing it as an operator acting on the objective function. More specifically, we view TSP as a mapping from L̂CELBO(y|u) to L̂CELBO-SP(y|u) = L̂CELBO(y|u) ◦ TSP and discuss the relationship of these objective functions.
Let ζ := (φ, ψ, ξ) denote the tuple of the variational parameters for brevity and Ω be the set of the parameters (θ, ζ) on which we perform the optimization. Let θ∗ denote the true generative parameter of (y,u), which is not necessarily contained in Ω. Finally, define the discrepancy function of (θ, ζ) ∈ Ω with respect to θ∗ by
δ(θ∗‖θ, ζ) := DKL(y|u)(θ∗‖θ) + ∆(θ, ζ),
where DKL(y|u)(θ∗‖θ) := E(y,u)∼θ∗ [ln p(y|u,θ ∗)
p(y|u,θ) ] is the conditional Kullback–Leibler divergence of θ∗ and θ and ∆(θ, ζ) = ∆(θ, φ, ψ, ξ) is the gap function defined just after Equation 12. Note that it is nonnegative and takes zero if and only if both θ and ζ are correct with respect to θ∗, i.e., δ(θ∗‖θ, ζ) = 0 if and only if
p(y|u, θ) = p(y|u, θ∗), q(z|y,u, φ) = p(z|y,u, θ), q(z|u, ψ) = p(z|u, θ), f(u; ξ) = ln p(u|θ),
for all y, u and z. Note that the correctness of (θ, ζ) is measured with respect to the predictive form of the true model p(y|u, θ∗), not the generative form p(y,u|θ∗), which is an essence of DIG. In other words, δ(θ∗‖θ, ζ) measures the predictive discrepancy of (θ, ζ) from θ∗.
Justification of CELBO maximization. Observe that
E(y,u)∼θ∗ [L̂CELBO(y|u)(θ, ζ)] = −h(y|u)− δ(θ∗‖θ, ζ),
where h(y|u) := Eθ∗ [− ln p(y|u, θ∗)] is the differential entropy of y given u. Since h(y|u) is independent of (θ, ζ), the CELBO maximization is justified as the discrepancy-function minimization,
maximize (θ,ζ)∈Ω E(y,u)∼θ∗ [L̂CELBO(y|u)(θ, ζ)] ⇔ minimize (θ,ζ)∈Ω
δ(θ∗‖θ, ζ).
Justification of CELBO-SP maximization as surrogate. Similarly, we have
E(y,u)∼θ∗ [L̂CELBO-SP(y|u)(θ, ζ)] = −h(y|u)− δSP(θ∗‖θ, ζ),
where δSP(θ∗‖θ, ζ) := δ(θ∗‖TSP(θ, ζ)) is referred to as the transformed discrepancy function. Then, the CELBO-SP maximization is seen as the transformed-discrepancy-function minimization,
maximize (θ,ζ)∈Ω E(y,u)∼θ∗ [L̂CELBO-SP(y|u)(θ, ζ)] ⇔ minimize (θ,ζ)∈Ω δSP(θ ∗‖θ, ζ).
Moreover, the transformed discrepancy function is consistent with the original discrepancy function in the following sense. Proposition 5. For all θ∗, θ and ζ,
δ(θ∗‖θ, ζ) = 0 ⇔ δSP(θ∗‖θ, ζ) = 0. (18)
Proof. It is shown as a corollary of Proposition 2.
In other words, δSP(θ∗‖θ, ζ) also measures the predictive discrepancy of (θ, ζ) from θ∗ in a different way. This justifies the CELBO-SP maximization as a surrogate of the CELBO maximization.
Limitation of CELBO-SP maximization as surrogate. A property stronger than the consistency (Equation 18) is the domination of the discrepancy function. Here, we say δSP dominates δ if there exists C <∞ such that
δ(θ∗‖θ, ζ) ≤ CδSP(θ∗‖θ, ζ). (19)
for all θ∗ and (θ, ζ) ∈ Ω. If the domination holds, then the CELBO-SP maximization implies (in expectation) the CELBO maximization up to a multiplicative constant.
Unfortunately, however, this is not the case in general. It is partly by design since the transformed discrepancy function δSP(θ∗‖θ, ζ) is derived as a result of suppressing the divergence of the density ratio wθ,ψ,ξ(u, z), which also causes the divergence of the original discrepancy function δ(θ∗‖θ, ζ). The following proposition shows that there exists no such constant C <∞ satisfying Equation 19. Proposition 6. There exists a triple (θ∗, θ, ζ) such that δ(θ∗‖θ, ζ) =∞ and δSP(θ∗‖θ, ζ) <∞.
Proof. It suffices to take (θ∗, θ, ζ) such thatwαθ,ψ,ξ(u, z) is not integrable with the density q(z|u, ψ), but {wθ,ψ,ξ(u, z)G(wθ,ψ,ξ(u, z))}α is integrable with the same density. For example, assume z takes a value in a Euclid space and take p(z|u, θ) ∝ exp(−‖z‖) and q(z|u, ψ) ∝ exp(−‖z‖2).

Review summaries:
The paper addresses the problem of predictive modelling with missing input features. The authors formulate the problem as a latent variable model, and in addition to the standard variational lower bound (ELBO) propose to use a variational upper bound based on CUBO (Dieng at al 2017) modified by an exponential divergence to solve the MC estimation of CUBO. They further propose surrogate parametrization to reduce the variance in the gradients. The experimental evaluation over standard regression UCI datasets with randomly dropped features shows marginal improvements over existing baselines.

The paper considers the problem of prediction with missing (incomplete) features. The authors propose a class of generative models that includes missingness of features, and develop a discriminative learning algorithm that maximises the conditional (posterior) log-likelihood of the training data approximately. Experiments show that the method is competitive compared with existing approaches, and in particular with approaches based on VAEs. 

The authors propose a new method, DIG, for discriminative tasks with missing input features. It uses latent variable models to marginalize out the label and the missing part of the features given the latent variable, in order to compute the objective, the conditional log likelihood of the label given the corrupted features. As this objective is intractable, the paper builds a conditional evidence lower bound (CELBO) that can be unbiasedly approximated using Monte Carlo samples. CELBO consists of the regular ELBO as the lower bound for the log joint probability of the label and the observed features, and an evidence upper bound (EUBO) that bounds the log marginal probability of the observed features. The derivation of EUBO involves the alpha-renyi divergence and the exponential divergence. The stochastic CELBO contains a density ratio that can lead to large variance in the stochastic gradients during optimization, so the authors propose a surrogate parameterization to bound the gradient norm. Experiments on real datasets justify the effectiveness of the variational approximations to stabilize the optimization. When compared with VAE, CVAE, and MICE, the DIG algorithm shows better or comparable predictive performance and robustness against feature corruption.

This paper proposes a new method for learning with missing data. Compared with previous approaches, the authors choose to perform discriminative learning with generative modeling so as to borrow the benefits from these two types of methods. To optimize with the underlying intractable loss function, the authors start from the traditional variational lower bound ELBO and one upper bound CUBO from a previous work (the $\chi$-divergence lower bound [1]) and derives a lower bound for the original loss function. To solve the issue with the estimation bias as well as the potential huge variance, the authors change the divergence function in CUBO as well as add the surrogate parameterization so that the Monte Carlo estimation of the loss can be unbiased and (potentially) with smaller variances. Experiment results show the proposed methods run stably and perform comparably or better compared to baseline methods.  References:  [1] Adji Bousso Dieng, Dustin Tran, Rajesh Ranganath, John Paisley, and David Blei. Variational inference via χ upper bound minimization. In Advances in Neural Information Processing Systems, pp. 2732–2741. 2017.