Variational Inference for Discriminative Learning with Generative Modeling of Feature Incompletion
Authors: Kohei Miyaguchi
Takayuki Katsuki
Akira Koseki
Toshiya Iwamori
Venue: ICLR 2022 Oral (2022)
Keywords: Black-box variational inference
missing values
evidence upper bound

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.


Published as a conference paper at ICLR 2022
VARIATIONAL INFERENCE FOR DISCRIMINATIVE LEARNING WITH GENERATIVE MODELING OF FEA- TURE INCOMPLETION
Kohei Miyaguchi, Takayuki Katsuki, Akira Koseki & Toshiya Iwamori IBM Research- - Tokyo m miyaguchi gibm. com, {kats, akoseki, iwamori iip. ibm. com
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-valus imputation and regression, simultancously or sequen- tially, which 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 favor- able 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 varia- tional inference framework so that it can be applied wide variety of modern machine learning models, including the variational autoencoders. We also con- firmed the effectiveness of the proposed method empirically.
INTRODUCTION
We address the problem of prediction with incomplete data, which often arises in real-life data due 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 charac- teristies 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, c.g., non-standardized medical equipment, legal regulation and privacy concerns. More formally, each input is an incomplete set of features represented with tuple that & E dd is the corrupted version of the complete input x E Rdd and m E f0. 1 1]dd indicates the missing entries of x: For all [dz], my 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;
(X, m) such
predict the outcome y Rdd, which is modeled with a predictive distribution ~ p(y u, 0). Here, the parameter 0 is unknown and should be learned from data.
Xxy) Given the incomplete features
One straightforward way to deal with such incompletion is to incorporate generative models: The missing values 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 that the meaning of is inherently incorporated. Without such context, algorithms must learn by themselves that 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 the feature incompletion.
On the other hand, the discriminative models that directly learn the conditional density pe (y|J)) sometimes exhibit superior performance (Kuhn et al. 2013) due to their alignment in terms of the ob- jective. This is especially the case when the feature incompletion is not a dominant factor (Lasserre et al. 2006). However, it is not straightforward generative models to incorporate the context of missingness into discriminative models and therefore their prediction may be less robust to the feature incompletion.

---

Published as a conference paper at ICLR 2022
To best of our knowledge, there have been general solution to this dilemma of the generative and the discriminative approaches, Le., generic method of discrininative inference with incomplete features for generative models, which abbreviate 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 (Ghahra- mani & Jordan [1994] Smola et al. 2005) and Gaussian processes (Pacheco et al. 2014). but they cannot handle modern machine-learning architectures such 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 -ox 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 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 we review and discuss the related work. Finally, in Section [ conclude the paper with some remarks and future perspectives. All the proofs for theoretical statements are given in the appendix.
PRELIMINARY
In this section, we introduce the mathematical formulation of the DIG method and then review some background of its technical challenges.
2.1 DIG FORMULATION
Let E [n]] denote the index instances, where n is the number of data and [n]] denotes the set of
integers (1, n.).
Let yll == {y² R& heja and xn == (x' E Rd Jeju] be the sets of the target and the complete feature vectors, where dy and d, are the dimensionalities of vectors. We assume we can only observe the set of incomplete feature vectors u!d) {u' (x² m'))) instead of [() where x' E Rd is the corrupted version of x' and m E {0, 1} de is the indicator of corruption ('mask vector') ) for all i E [n]. Each mask vector m² indicates which entry of x' is corrupted; We have xj. = (1 m )x] for all i E [n] and all j [dz].
Generative model. Suppose that the generative process of the observables (yl"), 1l") modeled with a parametrized family of instance- wisee identical distributions. That is, the joint probability density of (yfml, u!"l) is factored as
p(ybn ml u 00 = II ply u'10),
(I)
sE(n]
where the function p(y', u' u'l0) on the RHS denotes the common density function across instances i E [n]] indexed with any parametric models 0 (c.g., by neural networks).
Discriminative objective. rupted data,
The goal of DIG to maximize the conditional likelihood of the cor-
(n) ( (0) In iny(y! In .0
In Inp(y"lu',on
(2)
SE|n
Now, let (0) == Inp(y"lu' 8) denote the i-th summand. Following Bayes' rule, each 2 (1) is computed with a difference of log likelihoods
f (0) = (.) ())
where C2 (0): := In ply', u (8) and C (() := In / p(y", u'|0) dy Since we employ gradient- based optimization methods, it suffices consider (approximate) evaluation of the value and the gradient
2

---

Published as a conference paper at ICLR 2022
of the instance- wise objectives (0) separately. Thus, in the following, limit our focus to the separate objective and omit the index if there is risk of ambiguity.
2.2 DIG WITH LATENT VARIABLE MODELS (DIG-LVM)
We give further details of DIG with specific class of generative models, namely the latent variable models (LVMs). LVM one of the most common and expressive generative models. In the context of DIG, the generative density function is given by
ply, u / ply, u, z00) dz,
where z is the instance-wise latent variable, which is expected to capture some higher-ordes order infor- mation of the observables. For example, the variational autoencoders (VAFs) (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 gen- erative processes can be modeled with relatively simple joint density p(y, u, z z()). This is especially beneficial in modeling missing data: The full joint density ply. u. 2|0) can be easily designed to reflect one's 's belief the corruption process such 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
ply. u, 2 p(m) p(z) p(y|2,0) 0) II (z(x)12,0))1
(5)
where p(m) is the true marginal density function of m, which is unknown but cancels out in the DIG objective (Equation[3 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
de
ply, u, 0 = p(m) p(z) p(y|z, m. 0) II {p(x, (0,zzz.m.00 0)) j=1
1)
Note that these are just examples and 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 relegated to the appendix (Section G
Under these factorizations, the log likelihoods on the RHS of Equation 3 for LVMs are computed as
(0) =In / p(y, u, 14,200 dz,
(.) In p(u, z00) dz,
where both density functions in the integrals can be analytically computed under either of Equation 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 ct al. 1999 Zhang et al., 2018) in favor of flexibility and scalability. However, with DIG, also have a negative log integral, i.e., -Lu (f). The optimization of such objective is relatively difficult since there have been equivalent of VI in terms of the flexibility and scalability (see Section 5)
VARIATIONAL. INFERENCE FOR DIG-LVM
Now we present the main result, the method of VI for DIG-LVM. In Section 3.1 tional approximation to the conditional likelihood given by Equation 3 Then, in Section 3.2 and 3.3 derive method for optimizing the variational approximation and its variant for variance reduc- tion, respectively. Finally, in Section 3.4 we show the prediction procedure based on the optimized parameters.
we derive varia-
3

---

Published as a conference paper at ICLR 2022
3.1 VARIATIONAL LOWER/UPPER BOUNDS
In this section, we discuss an approximation of L0, and Ln(0) in Equation Zindividually. Note that both log integrals of the specific form
L (0) p(v, 2(V,
(8)
where v represents either (y, u) or u. Since such integrals constitute the objective (()) with positive and negative signs, our goal is to derive both upper and lower bounds Equation
Evidence Lower Bound (ELBO). The lower bound is given in the standard way (Jodann et al. 1999
Zhang et al. 2018).
(v) 2 (0) (0110)
In p(v,z00)
= g(
(ELBO
(9)
((llv,
where ④ is variational parameter of probability density function q(zlv, Φ) and g(xlv.e) (e) [In p(a is the KL divergence of the parameters Φ and 8, We call (xl(( (01(0)== (oll() LELBO(v) (0, Φ) as the evidence lower bound (ELBO). Since ELBO is expectation of tractable function, we approximate it with Monte-Carlo sampling.
Evidence Upper Bound (EUBO). To derive upper bound, start with applying the x- evidence upper bound (CUBO) (Dieng et al. 2017). For any real numbers a > 1, CUBO is derived follows:
(0) < (0) + (1 (a
1
plv. ((2/v, v))
In Eeng(el)))))
LcuBo(v) (+)
(10)
where 10 is variational parameter of a probability density function q(zlv, v)) and D Da(aiv) -- In pz (ziv, denotes the a-Rényi divergence of the pa- rameters and wh Note here, unlike ELBO, CUBO is not unbiasedly approximated because of the logarithm wrapping the expectation. To address this issue, we apply another variational approxi- mation with a divergence function (t) we have
:- at - 1)/a,t E R: Since a 2 0 for all f E I,
(cUBO(v) (0, < fcUBO(v) (0, Va (a(Ccuuo(v) - f(v:E))
e f(v:(
p(v,2 g(z(v, v))
Ba (-vaa (.
+ =: LEUROI (0,4.
(11)
where E is 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 = (y, u) and EUBO on V = u, we have conditional evidence lower bound (CELBO),
(CELBO(y) (0, o, v. C0 LEUBO( E)
(12)
By definition, CELBO bounds the DIG objective from below, Lyyu(() 2 LCELBOCY (0, di th, E The inequality is tight for any generative parameter 0 owing to the tightness of ELBO and EUBO i.c., for all 6, there exists a tuple of variational parameters (ó, v. such that the gap A o, v, E) Cy
LCELBOCY u) (0, á, 1h, E) zero. Moreover, CELBO can be unbiasedly approximated as well as ELBO and EUBO.
'ELBO is tight if q(llv, +)) 0). EUBO is tight if g(zlv, p(z|v,0) 0) and f(v:E) CeuBorv) v)) In L.
4

---

Published as a conference paper at ICLR 2022
Algorithm 1 Variational Inference for DIG (vDIG)
Input: Data (ylmm un))
Output: 0, 3
I: (0, Φ, () Initialize( // Any initialization methods can be used.
repeat
3: Draw minibatch C [n]
<- TT' EiEn Limaooty (0, E) (сι is given by Equation 13 Use Equation
IS
instead for variance reduction.
(0, Φ E) +- e (10, o, v5 E): VL) // Any gradient- optimization methods can be used.
6: until converge
3.2 OPTIMIZATION ALGORITHMM: VDIG
Since CELBO be unbiasedly approximated, we may employ stochastic gradient-based optimiza- tion to maximize it. The full stochastic objective for CELBO is given by
(CELBO (0, E) =n ply, 2 (0)
p(u, z210) g(zz) u, ()
a
((: + + (13)
g(zøly, u1 ()
where 2 and Z are Monte-Carlo samples drawn from gizly, u, Φ) and g(zlu, 16), respectively.
The gradients of LCELbory 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 (CELBOGy)J)) LCELBOGyj over all the instances, we may draw minibatch of instances for each iteration. We call the resulting inference algorithm vDIG (Algorithm) D
3.3 VARIANCE REDUCTION FOR VDIG WITH SURROGATE PARAMETRIZATION (SP)
The boundedness of the norm the stochastic gradient V (0, Φ, +) is crucial for the stable and fast convergence of stochastic gradient-based algorithms like Algorithm However, the stochastic CELBO contains the density ratio
wo.v.c (u, = p(u, z
(14)
q(z) u, (u:E)
raised to the power of which is problematic the ratio wa 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 have wo,v.e(;:)) = whenever the objective gap is zero. That is, constraining the parameter to satisfy
sup,,x W8. (u, < 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 surrogate parametrization. Define new parameters 6' and g formall by
ply. =,, G(wo v (0.,e((,))) (),2)) ply, u1 20)),
f(u:g') f(a:) è - In 2(0').
Z(0")
where G(w)= :1 v w) + a In(1 V w)) (w 2 0), a avvbb= == max(a, [a,b], b}. and Z(6') is the normalizing constant ensuring the mass preservation of the density under 0' We refer to the mapping TsP (0, o. v, {) + (0', Φ, w. E') as the surrogate transform, and G(w) as the gain function. The surrogate transform Tsp and the gain function G(w) are designed carefully satisfy the following two properties.
First, it preserves the effective parameters of CELBO.
Note that these parameters are conceptual objects and there is concrete implementation for them in the final algorithm.
5

---

Published as a conference paper at ICLR 2022
Definition 1 (Effective parameters), Let Ω be of CELBO parameters (9, E).. Then, define the effective parameters of S2 by Po(2))=((0 o, v, E Q ( v, v, = 0), i.e. the set of parameters inducing tight variational approximation.
Proposition 2 (Effective parameter preservation). For arbitrary Q we have
Oo(
In other words, in sense, Tsp does alter the original model. See Section for extended discussion.
Second, regularizes the growth of the highly stochastic term (Equation 14p. 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,
Wo0 & (u, 2) G(ue, wi 2)).
W0,0ig (u,z)
Note that G(w) is no larger than one and non -oniincreasing for w 2 O, thus the stochasticity is reduced through the transform. See Fn the appendix for visualization. Consequently, TsP guarantees the boundedness of the stochastic gradient.
Proposition (Bounded gradient). Assume there exists | In p an 10)1, || g(zoly d)1l, g[ v)ll, [V; S [V(CCLLBo(y uj (O, o, vh E) s 9K.
such that [| Inply. u, 2.z00))
K. Then,
Proposition 2 and 3 justifies optimizing the objective via the transform Tsp; we get a gradient- norm bound without altering the effective parameters. The full objective function after the surrogate transform, (CELBO-SP(Y()) = LCELBO(y| TsP, given by
ply, 2 0 ( 26))
m () In
q y, u)
( (u, Zu G( (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 2(0'), which is intractable in general. This is not the case with the surrogate transformation ELBO or EUBO alone. The CELBO-SP maximization is done by simply replacing CELBO with CELBO-SP Algorithm
3.4 PREDICTION ALGORITHM
Given (8, ó. sb, {) trained by Algorithm I we want to compute the predictive distribution on new instances given their incomplete features new = (x"ew mow). Since the conditional density po(y/u) u) is intractable in general, we approximate with the Monte-Carlo method. The approxi- mated conditional distribution is given by
p(ylu, 0)
Σ sE(Kymel
ply|z 0),
Xipred 2 1,
(16)
Repred
where zi. independently drawn from (z)u), & [kiree).. This procedure is justified follows.
Proposition 4. Let ply 0) Elp(y/u, 0)}, where the expectation is taken with respect to the Monte-Carlo sampling. Then,
DKLLyju < a
In other words, if the objective gap is small, so is the approximation error of 0, which is the limit of the actual predictor 0 with kireed + oc. In the experiment, we used kpred 512.
4 EXPERIMENTS
In this section, we demonstrate the effectiveness of the proposed method, i.c., Algorithm I through numerical experiments. We first introduce number of methods compared in the experiments in

---

Published as a conference paper at ICLR 2022
Section 4.1 We then present the procedures and results of three experiments designed to show i) the effectiveness of the new variational upproximation, ii) the superior performance of the proposed method, and iii) the robustness of the proposed method against feature corruption.
4.1 EXPERIMENTAL SUBJECTS
We employ three subject algorithms and two baselines in the experiment. The main subjects are all based on VAE (Kingma & Welling, 2013). 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 a part of the feature, x' < x ① y, and reconstructing x' from its corrupted version, where the entries cor- responding to y is considered (2020), which comes with two variants for different missing-value processes, namely missing not random (MNAR) and missing completely random (MCAR). The difference between MNAR and MCAR is only in the generative model (i.c., the architecture of the decoder), represented by Equation and Equation 5
missing. In particular, we adopt the formulation of Collier et al.
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-basee formulation.
Discriminative method: CVAE. As an example of discriminative approach, we employ the con- ditional VAE (CVAE) proposed by Kingma et al. (2014); [ohn al. (2015). In setting, the conditional variables consist of the corrupted feature x and the missingness indicator m. In prac- tice, these vectors are concatenated before fed into neural networks. The generative model of CVAE in this setting is given by p [,, 0) == f p(y; z z/u, 0) dz, while the decoder represents the density in the integral and the encoder used approximate the posterior p(z y, 0). Note that the sulting model is not informed with the context of m as VAE and VAE* are 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 discrimina- tive VAE (DVAE). We implement MNAR (Equation () and MCAR (Equation (5) variants of DVAE 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 0 and the encoder common. See Section D for more details. We choose 2 for the parameter of the Rényi divergence.
Baseline methods: Simple, MICE. As the first baseline method, we employ simple fully- connected feed-forward neural network (FCFNN) with the same architecture with the encoders of the above methods, except its output considered as distribution of the target y instead of z. Although the output distribution is restricted to Gaussian, it requires no variational approximation and hence the conditional likelihood In p 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 statisties to deal with missing values. In particular, we em- play the Bayesian ridge regression for the missing-valle value imputation task and the FCFNN of the first baseline for the subsequent regression task.
4.2 EXPERIMENTAL PROCEDURE AND RESULTS
The organization of the experiments is three-fold. First, Section 4.2.1 we check the effective- ness of our variational approximation techniques, namely EUBO and the surrogate parameteriza- tion. 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.
Because of this, there is MCAR variant for CVAE.

---

Published as a conference paper at ICLR 2022
ELBO=CUBO CUBO
ELBO-EUBO
ELBO+EUB0++P
Taking L2St
4-
10
20
10
20
10
20
Time Eserl
Time (sec]
Time [sec]
Figure 1: Each of three panels corresponds different variational approximation. Lines in each panel denote the training objectives during 5 independent training runs.
AQ-CO
AQ-NMHC
AQ-NOx
Boston
Diabetes
YearPred
Total
CVAE
0.41 (0.02)
(0.14)
5.33 (0.02)
3.01 (0.18)
5.46 (0.03)
3.49 (0.06)
3.82 (0.04)
DVAE
0.45 (0.03)
(0.10)
5.26 (0.02)
2.91 (0.23)
5.44 (0.08)
3.36 (0.04)
3.76 (0.05)
DVAE*
0.44 (0.03)
(0.11)
5.27 (0.03)
2.92 (0.21)
5,42 (0.07)
3.35 (0.02)
3.76 (0.04)
Simple
0.49 (0.03)
(0.11)
5.40
3.13 (0.18)
5.47 (0.04)
3.62 (0.04)
3,93 (0.04)
MICE
0.46 (0.06)
5.34 (0.23)
5.38 (0.02)
2.94 (0.10)
5,46 (0.03)
3.61 (0.08)
3,87 (0.05)
VAE
0.47 (0.01)
5.21 (0.08)
5.49 (0.06)
2.95 (0.23)
5,48 (0.10)
3.59 (0.02)
3.87 (0.04)
VAE*
0.46 (0.02)
5.20 (0.09)
5.47 (0.05)
2.81 (0.20)
5,46 (0.05)
3.58 (0.03)
3.83 (0.04)
Test per-sample cross entropies of different algorithms for different datasets. The num- bers represent the mean (and standard deviation in parenthesis) of five independent runs for each configuration. For each column, the best score is indicated with bold face.
Table 1:
4.2.1 EFFECT OF NEW VARIATIONAL APPROXIMATION
Procedure. We compare the stability optimization of DVAE with three different variational approximation, namely, ELBO+CUBO, ELBO+EUBO and ELBO+EUBO+SP. In all cases, the Monte-Carte Carlo approximation of ELBO is used to compute the joint evidence (0). The marginal evidence Lu (0) is computed with the Monte-Carlo approximation for both CUBO (albeit biased) and EUBO. Only ELBO+EUBO+SP employs the surrogate parametrization.
Result. Figure shows the training processes of DVAE with the AQ-CO dataset from UCI Ma- chine Learning Repository (see the appendix). It is 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 Equation 10. but cannot know how much.
4.2.2 PREDICTIVE PERFORMANCE
Procedure. We apply seven different algorithms in Section 4.1 to six different regression tasks from UCI Machine Learning Repository, summarized in Table (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 I shows the entropy, i.e., - In po(yjuu u), averaged over the test splits. Overall, both DVAE and DVAE* almost always perform best or at least comparable to the best scores within standard deviation, indicating the effectiveness the DIG method applied VAE. In particular, is more clear that they outperform the others when averaged over all the datasets ('Total' column). Moreover, DVAE * is slightly better than DVAE expected because the way we drop the feature entries is MCAR.
4.2.3 ROBUSTNESS AGAINST FEATURE CORRUPTION
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.
8

---

Published as a conference paper at ICLR 2022
0.400
Diference Cror Entepy
0.350
0.325
0.275
0:200
DVAE
Simple
MICE
Figure 2: The average performance loss when the missing-value ratio is changed from = to p .5. Lower values imply more robustness against feature incompletion. The error bars represent the estimated standard deviation of the expectations.
Result. The difference is visualized in Figure 2] The most robust algorithms 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* informed with the context of and implicitly regu- larized 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.
5 RELATED WORK
In this section, review existing variational approximations bounding the fog-integrals from above. We also discuss the relationship with the existing studies on missing data handling in the appendix.
The x-variational upper bound (CUBO) (Dieng al. 2017) has been proposed to estimate such upper bounds, utilizing the a- a-Rényi divergence (equivalently, x-divergence). However, as Pradierr cct al. (2019) pointed out, the resulting estimate of CUBO is biased and in some cases numerically unstable This phenomenon has been also confirmed in experiment (Section 4.2.1p. Another upper bound has been proposed by (Ji Shen, 2019) with the reverse KL divergence instead of the Rényi divergence, but it is also biased and not guaranteed to be upper bound. We also find that the parsimonious upper bound (Mattei & Frellsen 2018) is suitable in terms of the upper-bound guarantee, while induces 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 (017)). which is designed for estimating partition functions. Although their method also suffers from high variance, they took different approach reduce it. In partieular, they adaptively reduced the learning rate of the variational parameter v). However, as the authors noted, this does not solve the problem enough for scaling large datasets.
CONCLUSION
We have proposed novel algorithm to perform discriminative training with incomplete features for generative models, which is derived on the basis of newly introduced variational approximation and parameter transformation. The effectiveness 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.
Another work (Lopez al. 2020 reported that there is case CUBO works just fine, its seems
problem dependent phenomenon.

---

Published as a conference paper at ICLR 2022
REFERENCES
Melissa J Azur, Elizabeth A Stuart, Constantine Frangakis, and Philip J Leaf. Multiple imputation by chained equations: what it and how does it work? International journal of methods in psychiatric research, 20(1):40-49, 2011.
T.
YearPrediicionMISD. UCI Machine Learning Repository, 2011.
Mark Collier, Alfredo Nazabal, and Christopher KI Williams. Vaes in the presence of missing data. arXiv preprint (rXXi:2006605300 2020.
De Vito, E. Massera, M. Piga, L.. Martinotto, and G. Di Francia. field calibration of an elec- tronic nose for benzene estimation in an urban pollution monitoring scenario. Sensors and Actua- tors B: Chemical, 129())750-757, 2008. ISSN 0925-4005. doi: https///bbi.rrg//000006...
2007.09.060 URL https://wwws www.::/wcenee.iicc rect com http:////eeeeec/eatt.cllle/////e
2092540050700765
Arthur P Dempster, Nan M Laird, and Donald Rubin. Maximum likelihood from incomplete data via the algorithm. Journal of the Royal Statistical Society: Series B (Methodological), 39(1): 1-22, 1977.
Adji Bousso Dieng, Dustin Tran, Rajesh Ranganath, John Paisley, and David Blei. Variational in- ference via X upper bound minimization. In Advances in Neural Information Processing Systems, pp. 2732-2741. 2017.
Craig K Enders. Applied missing data analysis. 2010.
Zoubin Ghahramani and Michael I Jordan. Supervised learning from incomplete data via an em approach. In Advances Neural Information Processing Systems, pp. 120-127, 1994.
lan Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio. Generative adversarial nets. Advances in neural information processing systems, 27, 2014.
Alex D Holub, Max Welling, and Pictro Perona. Combining generative models and fisher kernels for object recognition. In Tenth IEEE International Conference Computer Vision (ICCV'05) Volume 1, volume 1, pp. 136-143. IEEE, 2005.
Tommi Jaakkola and David Haussler. Exploiting generative models in discriminative classifiers. In Advances in neural information processing systems, pp. 487-493, 1999.
Chunlin Ji and Haige Shen. Stochastic variational inference via upper bound. In NeurIPS workshop on Bayesian Deep Learning, 2019.
Michael Jordan, Zoubin Ghahramani, Tommi S Jaakkola, and Lawrence K Saul. An introduction to variational methods for graphical models. Machine learning, 37(2):183-233, 1999.
Diederik Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980 2014.
Diederik P Kingma and Max Welling. Auto-encoding variational bayes. In International Conference on Learning Representations, 2013.
Durk P Kingma, Shakir Mohamed, Danilo Jimenez Rezende, and Max Welling. Semi-supervised learning with deep generative models. Advances in neural information processing systems, 27: 3581-3589, 2014.
Max Kuhn, Kjell Johnson, et al. Applied predictive modeling, volume 26. Springer, 2013.
Volodymyr Kuleshov and Stefano Ermon. Neural variational inference and learning undirected graphical models. In Advances in Neural Information Processing Systems, pp. 6734-6743, 2017.
Julia A Lasserre, Christopher M Bishop, and Thomas P Minka. Principled hybrids of generative and discriminative models. In 2006 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR106), volume 1, pp. 87-94. IEFE, 2006.
10

---

Published as a conference paper at ICLR 2022
Romain Lopez, Pierre Boyeau, Nir Yosef, Michael Jordan, and Jeffrey Regier. Decision-maiee makinn with auto-encoding variational bayes. arXiv preprint arXiv:2002.07221, 2020.
Mattei and Jes Frellsen. Leveraging the exact likelihood of deep latent variable models. In Advances in Neural Information Processing Systems, pp. 3855-3866, 2018.
Thomas Minka. Expectation propagation for approximate bayesian inference. In Proceedings of the Seventeenth conference on Uncertainty in artificial intelligence, pp- 362-369, 2001.
Ricardo Andrade Pacheco, James Hensman, Max ZwieBele, and Neil D Lawrence. Hybrid discriminative-gennenee approach with gaussian processes. In Artificial Intelligence and Statis- fics, pp. 47-56, 2014.
Alessandro Perina, Marco Cristani, Umberto Castellani, Vittorio Murino, and Nebojsa Jojic. A hybrid generative/discriiinatve classification framework based on free-energy terms. In 2009 IEEE 12th International Conference on Computer Vision. pp. 2058-2065. IEEE, 2009.
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 Symposium Advances in Approximate Bayesian Inference, 2019.
Danilo Jimenez Rezende and Shakir Mohamed. Variational inference with normalizing flows. arXiv preprint arXiv: arXiv:1505.00777 2015.
Donald B Rubin. Inference and missing data. Biometrika, 63(3):581-592, 1976.
Marek Smieja, Lukasz Struski, Jacek Tabor, Bartosz Zieliñski, and Przemystaw Spurek. Processing of missing data by neural networks. Advances Neural Information Processing Systems, pp- 2719-27201 2018.
Alexander J Smola, SVN Vishwanathan, and Thomas Hofmann. Kernel methods for missing vari- ables. AISTATS, 2005.
Kihyuk Sohn, Honglak Lee, and Xinchen Yan. Learning structured output representation using deep conditional generative models. Advances in neural information processing systems, 28: 3483-3491, 2015.
BETH Twala, MC Jones, and David Hand. Good methods for coping with missing data in decision trees. Pattern Recognition Letters, 29(7)9550-956, 2008
Ronald J Williams. Simple statistical gradient-following algorithms for connectionist reinforcement learning. Machine learning, 8(3-4):229-25) 256. 1992.
Cheng Zhang, Judith Bütepage, Hedvig Kjellström, and Stephan Mandt. Advances in variational inference. IEEE transactions on pattern analysis and machine intelligence, 21(8:22088. 2018.
Juntang Zhuang, Tommy Tang, Yifan Ding, Sekhar Tatikonda, Nicha Dvornek, Xenophon Pa- pademetris, and James S Duncan. Adabelief optimizer: Adapting stepsizes by the belief in ob- served gradients. arXiv preprint arXiv:2010.07466 2020.

---

Published as a conference paper at ICLR 2022
Giwl
Figure 3: The gain function (ww) with a = 2 (the orange line). The green line shows the size of the stochastic term wa, (u, z) after the transformation, whereas the blue dashed line shows the original size.
A PROOFS
A.1 PROOF OF PROPOSITION
Proof. (o(Tsp((2) S Oo(2) is trivial as (0, 6, b, {) = 0 implies wa.w.e = 1, which implies (0, o, th, is fixed point of Tsp, hence (6, o, v. {) E Tsr. The other direction is shown as follows. Assume A Φ, E.) = 0, where (6', Ø, t, E') = (sp o, (. and (0, Φ, t/, E Ω. Then, have
= wo w.e (u, 2) wo.se (u, 2) G ((
for all u and z. This implies Wa (t. = as wG(w) = 1 + w = 1. Therefore, (0, 0 th {) is fixed point of TsP and (0', 27 v, E Ω.
A.2 PROOF PROPOSITION 3
Proof. Let g= V(CCEBBO u) Tsp)((, o, )) and observe
g Vin ply, 11 2. (w)) V (Z y u. ()
(wG(w)) V)(
Vin In ply, u, (0
+) (H(nw)
g(z y, u, 0)
Vin ply, 00)
vf(u + H'In
((z y. u, 6)
where w:= wa (u, 2). H(t) := In (e'G(el))) /a, and (t(t) denotes the derivative of (t).. Thus, we have
|sl|| < |/ (I) [Vv +3K < 9K
as (7' S and | w|l SK.
A.3 PROOF OF PROPOSITION 4
Proof. According the information processing inequality, DKL(z/u) (0]l)). Moreover, by the construction of LMEUBO. we have (1 ) (0 (0l))) A0 é, v. E). The desired result is seen by combining these two inequalities with the fact that the a-Rényi divergence dominates the KL divergence for all a > 1.
have DKL(y|u) (01
B EXPERIMENTAL CONFIGURATIONS
Computing infrastructure. The computing infrastructure used the experiments is summarized
in Table
12

---

Published as a conference paper at ICLR 2022
CPU
RAM
GPU
PyTorch
Intel(R) (0-6670) CPU @ 4.00GHz
64GB
NVIDIA TITAN X
1.9.0
Table 2: Summary of computing infrastructure.
Preprocessing. All datasets used in the experiment are standardized before fed into algorithms so that the variables of the feature and the target have mean and unit empirical variance on the training split. More precisely, we apply the following transformation
8² + (1 - m²) 3 (x' (a)) 2 āz:
y' - (y' - (1y 2 6g
E [m]],
-1= ال (1-mm=) and 002 2n (1 m M(x -pe)) 002 for j E [dz]. and By ==
where Paj
z.j
Eis (I m))
1 EE=1/ 1² and 00² = 2 (v/ - Ay) Here, O and @ denotes the elementwise multi- plication and division of vectors, respectively. Note that we refill the missing entries of x 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 < (1 - m) @ (x - (Az) @ oy. and the predictive distribution p(y i,, 0) for the standardized target made based on it. Then, the distribution is transformed inversely predict the -onnstandddd
target, p(ylu u, 0) (II) ) p((y - (by) @ 8g u, 0).
Validation split. We randomly hold 20% of the training split out for validation. Algorithms re- ceive as the input only the remaining 80%.
Initialization. All the parameters involved in the experiments are contained in the Linear mod- of Pytorch and initialized with default initialization method, i.e., independently subject to the uniform distribution on the interval [- Vie JEl. where k is the number of input variables for the respective modules.
First-order optimizer. In all the experiments, employ the AdaBelief optimizer (Zhuang al, 2020) with its default setting 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... 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 (0, 2, 5, 8, 13, 19, 29, 43, 64, 94, 138, 203, 298, 436, 638, 934, 1366. 1998), which are (approx- imately) equally spaced in the log scale. Finally, the best model among these 20 checkpoints 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 knned 512 instances from the posterior distribution g( gz//u, v)) and make prediction with Equation 16
C DATASETS
All the datasets used in the experiments 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 corre- sponding to different targets denoted by their suffixes. YearPred is a part the YearPrediciooMMSD 20IT). where, to accommodates the fast iteration of the experiments, the number of records is restricted 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 Table 3
dataset
See https://///// / thub. wwt.ub.comjjunett..e. www.t//ttnnooooooooooo/eeooooaee whwpp/nnttnnnnnngggngaea zhuang/ Adabel fef- ii zer wwtpi//e//pdaaaaaaa 0.2.0
13

---

Published as a conference paper at ICLR 2022
AQ-CO AQ-NMHC
AQ-NOx
Boston
Diabetes
YearPred"
72
6139
731
6174
404
353
10000
10
10
10
13
10
90
Table 3: Summary statistics of datasets. 'x' indicates subsampled dataset.
D METHODS COMPARED
In this section, we show the detailed architecture of the methods compared the experiments.
Both the encoder and the decoder cally-connected feed-forward neural networks (FCFNN) with single hidden layer of width 256 and ReLU activation that take a vector input E IR., 2 1, and emit diagonal Gaussian distribution depending architectures, we first define some basic building blocks.
To facilitate the specification of the
Let every distinct occurrence of Li (a))) denote distinet linear layer with the output dimension > 1 with its own parameters (W, b) such that Li dineara(t = Wt + with a weight matrix W E Rdx* and offset vector b IRd for R The Gaussian layer is defined by the composition of the Gaussian density function and the linear layers,
Gaussa(tlv) Nait Lineara(v), In(1 + Lineara(v)))
where e = 1018101800 is small constant and Na[tlv. denotes the d-product of the probability density functions of the Gaussian distributions with mean V, and standard deviation s, evaluated ty. j E [d], for t, v E Rd and E 1R We also denote the ReLU activation function by ReLU(t) == max {0, t}. Finally, we define the composition of Gaussian layer and a Linear layer with & hidden neurons and ReLU activation function by
Gauss Causs-NNa, (t = Gauss, o Lineary((v))
Now, we are ready to describe the architecture specifications.
Simple. The objective of the Simple baseline is given by
Inp(y/u, 0) In Mauss-NNa,,red (y/ + m).
where & denotes all the parameter implicitly involved in the RHS.
MICE. MICE is implemented with Iterative Imputer from scikit- (ver. 0.24.2) with its default argument except with the burn-in period changed from 10 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. (020). The difference that the mask vector m is modeled as conditional variable rather than a random variable. The reason of this modification is because on such variables simplifies the resulting model. Namely,
not interested in the likelihood of and conditioning
LVAE σ) - Inply, ulz, 0) In 0, In u.
where
ply, u/z, 0) = p(m) Gaussa, (vlv).
II
(
jE,d,]m,,
q(zly, u, Φ) = Gauss e m).
(ReLJ Lineargss m) and q(zyy, o). For VAE*, replace with v = (ReLU o Linear Here, p(m) is an arbitrary density function of m, which is constant with respect the parameters (0, o) and ignored in the gradient computation.
14

---

Published as a conference paper at ICLR 2022
AQ-CO
AQ-NMHC
AQ-NOx
Boston
Diabetes
YearPred
Total
CVAE
0.76 (0.04)
5.62 (0.13)
5.82 (0.05)
3.32 (0.21)
5.56 (0.05)
3.57 (0.10)
4.11 (0.05)
DVAE
0.76 (0.03)
5,53 (0.08)
5.74()
3.31 (0.24)
5.53 (0.02)
3,44 (0.02)
4.05 (0.04)
DVAE*
0.76 (0.04)
(0.08)
5.75 (0.03)
3.23 (0.15)
5.50 (0.02)
3,43 (0.02)
4.03 (0.03)
Simple
0.87 (0.04)
5,87 (0.08)
5.96 (0.07)
3,43 (0.21)
5.57 (0.02)
3.71 (0.03)
4.23 (0.04)
MICE
0.95 (0.06)
6.07 (0.70)
6.01 (0.07)
3.80 (0.59)
5.61 (0.06)
3.68 (0.02)
4.35 (0.15)
VAE
0.80 (0.03)
5.53 (0.07)
5.88 (0.03)
3.37 (0.37)
5.57 (0.03)
3.55 (0.03)
4.12 (0.06)
VAE*
0.81 (0.04)
5.49 (0.10)
5.87 (0.04)
3.28 (0.40)
5.50 (0.04)
3.53 (0.04)
4.08 (0.07)
Table 4: The result of the same experiment of Section 4.2.2 except the missing vassueg-value ratio is = 0.5
instead of p = 0.1.
CVAE.
The objective of CVAE is given according to Sohn et al. (015).
LoVAE(0, In plyl2, 0) In Iz; 0, 1] In (z/(z u, o),
where
p 12, u, 0) == Gauss 250 (y S O m), g u, o) = Gauss @ x m).
and g(zly, u, Φ).
DVAE, DVAE*.
The objective of DVAE is given by Equation
13 where
ply, u|z, 0) == p(m) Gaussa, (ylv).
II
Gaussi (x,lv),
je sml
q u, Φ) == Gauss10(2|9256 (y, O2 x, m: w)),
(2/uu, v/= = (2) dy 1d, m;w)),
f(u: {) = In p(m) + Linear O0,y 1d, X. m))),
(t()(1) (2), (3) ww) = (
@ (2) @ (3) a (1) for (()))(()))))))))))))))))))))))
((2) E Rd, and
((() (4)) E Rdd, 3 denotes the shared parameter of d th, E. and v = (ReLUOLi nearsse) (z m)). For DVAE*, replace v with v (ReLU o Linear (z). Here, p(m) is an arbitrary density function of m, which cancels out in the final objective of CELBO or CELBO-SP.
E ADDITIONAL RESULTS
We show in Table 4 the results of the same experiment as in Section 4.2.2 except with the increased missing-var value ratio p 0.5.
F ADDITIONAL. DISCUSSION ON RELATED WORK: MISSING-DATA HANDLING
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 pair- wise deletion methods, single imputation methods and multiple imputation methods. The key fea- ture of this approach is that one processes the incomplete features to get estimates of complete ones & R x in the first step and then performs prediction based on & The score-based methods (Jaakkola)
& Haussler 1999: Holub et al. 2005 Perina et al. 2009 Smieja 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 pa (x), which is used to extract some information called score, s. The score is then fed to predictive models (yy (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 & is original or imputed. This makes it difficult to estimate the uncertainty resulted from the feature corruption. Moreover, even though the uncertainty
15

---

Published as a conference paper at ICLR 2022
Approach
Modeling
Objective
Missingness-Aware
Objective Alignment
Two-step
s(u), po(y]s) po u
p (s())
+
Generative
In p (y, u) In py
Discriminative
po(y/u) Po u)
DIG
(ypuly
Table 5: Summarized comparrson of related work.
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 such method is referred to as the full information maximum likelihood methods (see also Secion 4, Enders (2010), where the joint complete-dee corrupted elements of features obtain the objective In po (y,y, 1y, u). This approach may suffer from unnecessarily performance degradation as in the score-based approach, since the objective contains In (g)).
distribution po x) is explicitly modeled and marginalized over the
the generative term of u, In pe (y, u) In pe ( ())
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, 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, & O m, as the input 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 no generative term unlike in the two-step and generative approach. However, there is no trivial way to inform the model that 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 hybrid of the generative and discriminative ap- proaches. With DIG, the data is modeled with joint distribution pa u), but the learning objective is the conditional evidence In po(y/u), which is computed from Bayes' rule. Therefore, it naturally incorporates the information of missingness and is directly trained to maximize the predictive perfor- mance. The first application of the DIG strategy in the context of incomplete feature is Ghahramani & Jordan| 11944) which was followed by Smola et al. (2005) with a kernel-based generalization. It is crucial in their results that the complete-data model po xo.y.) is 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 derived by Pacheco et al. (2014) with the combination of the variational inference (Jordan et al.] 1999) and the expectation propaga- tion (Minka (001) framework. As opposed to their method, our focus on black-box variational inference algorithm applicable to variety of models.
See Table for the summary of the comparison.
G MODELING THE MAR PROCESS
In this section, we give method of modeling the third type of the (eature-corruption processes, the MAR process. Let O C [dz] be an index set on which the feature is always observed (i.e., my 0
for all j O) and define X = {x, O). Then, the MAR process is modeled with
ply, u, z ()) p(xo, m) p(z) p(y|z, xo, 0)-
II
{p(x, z, z,xo, 8)) 1mm
(17)
jE[d,IVO
where pxxo, m) will cancel out as well as pmm) in the MNAR and MCAR processes.
H ADDITIONAL JUSTIFICATION AND LIMITATION OF THE SURROGATE TRANSFORM
We present another justification of the CELBO-SP maximization (i.c., optimization of Equation [5] an alternative to the CELBO maximization (i.e., optimization of Equation 13). Recall that
16

---

Published as a conference paper at ICLR 2022
CELBO-SP is derived by the surrogate transform Tsp and Proposition 2 shows desirable prop- enty of Tsp viewing it as an operator acting on the parameters. In this section. we show another desirable property of Tsp viewing it an operator acting the objective function. More specif- ically, view TSP mapping from (CELBO(y)u) LCELBO-SP(yJu) LCELBO(ylu) Tsp and discuss the relationship of these objective functions.
Let <= v5, <) denote the tuple of the variational parameters for brevity and be the set of the parameters () on which we perform the optimization. Let 0* denote the true generative parameter of (y, u), which is not necessarily contained in f2. Finally, define the discrepancy function of (0, () E 52 with respect to 0 by
8(0" (), DKL(yin) (0" (10) ((,,C)
p(ylu.0) is the conditional Kullback-Leblerr divergence of 0* and and ((,, () A(0, Φ, 1h, 5) is the gap function defined just after Equation 12 Note that it is nonnegative and takes zero if and only if both 0 and are correct with respect to 0", i.e., 8(0 |10, = 0 if and only
where DKL(yluj (0 = -) (In p(ylu,")]
p(ylu 0) ply u. 0*),
aizly, Φ) - p(zly, 0) f(u: = In p(u
q u, vb) = p(z|u, 0),
for all y, and z. Note that the correctness of (0, is measured with respect to the predictive form of the true model p(y u, 8"), not the generative form p(y, u u which is an essence of DIG. In other words, 6(0* (|e, measures the predictive discrepancy of () from l".
Justification of CELBO maximization. Observe that
uyy [ACELBOCy]u) (0, --y u) 800"110,
where (y) = Eg-1 |- lly/ 0*)] is the differential entropy of given Since h(y) u) is inde- pendent of (0, <), the CELBO maximization is justified as the discrepancy-fumat funncion minimization,
maximize E (yyya) (ACELBOTy] ( (0, O) (0,C)ER
minimize 8(0" (10, C). (0,C)E
Justification of CELBO-SP maximization surrogate. Similarly, we have
By (P (0, {)] u-(// dsp (0"110
the transformed discrepancy function. Then, the CELBO-SP maximization is seen as the cuunforredddiisceppaee--tteeeeed
where dgp (8 (10 8(0" TTsr (0. ()) is referred
minimization,
maximize Ey.,u)= (cCELBO (0, () (9,))
minimize dsp (6" [|A, C)- (9,C)E
Moreover, the transformed discrepancy function is consistent with the original discrepancy function in the following sense.
Proposition 5. For all 0", 9 and Ç.
8(0" () =
8sp (00 0.
(18)
Proof. It is shown as a corollary of Proposition
In other words, (6" [l+CC)
(00 () also measures the predictive discrepancy of (0, () from 6* in a different way. This justifies the CELBO-SP maximization surrogate of the CELBO maximization.
Limitation of CELBO-SP maximization as surrogate. A property stronger than the consis- tency (Equation 18) is the domination of the discrepancy function. Here, we say dgr dominates 5 if there exists oc such that
5(0" (10,c) () < Cōsr 10
(19)
for all 8* and (0, <) E Q. If the domination holds, then the CELBO-SP maximization implies (in expectation) the CELBO maximization up to multiplicative constant.
17

---

Published as a conference paper at ICLR 2022
Unfortunately, however, this is not the case in general. It is partly by design since the transformed
discrepancy function (f [[,,5 <) is derived as a result suppressing the divergence of the density ratio wa (u, z), which also causes the divergence of the original discrepancy function 8(8" ||0, () The following proposition shows that there exists no such constant C' < oo satisfying Equation 19
Proposition There exists a triple (0", 0, () such that 8(0" 10, <) = oo and dgp (0" |10, <) oo.
Proof. It suffices take (0",0, () such that whe (u, z) is not integrable with the density g(z)u, u, v²),
but (we.vve (u, () w. (u, } ) is integrable with the same density. For example, assume
takes value in Euclid space and take p(z 0) x exp(- ||z]]) and q(z u, 1h) x exp(- ||z]])). O
18


Reviews:
[{"review_id": "wMV5F4agxNA", "reviewer": "Reviewer_WZsz", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "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": "", "main_review": "Strengths:\n\n1. The idea of learning missing data using discriminative learning together with generative modeling is interesting. As mentioned in the paper, performing such kind of learning will resulting a loss function as a subtraction on two integrals with respect to the latent variables, which makes it harder to derive a lower bound compared to the traditional variational inference cases. To solve the term being subtracted, the authors found a upper bound that could be estimated in an unbiased way with Monte Carlo methods.\n\n2. The exponential function in the first version of CELBO will potentially has a bigger variance when estimated using Monte Carlo. To solve this issue, the authors propose adding a regularization to the loss function while remain the optimal solution unchanged under the zero gap case. This also helps a lot in making training more stable, as shown in the experiments.\n\nWeakness:\n\n1. The topic of this paper is to focus on missing data. However, this paper does not put enough efforts on learning various cases of data missing patterns. As in the experiments, the authors only test the case where the data are missing completely at random (MCAR), which may not be the most common case in reality. MNAR case might be a more interesting situation to study with. The proposed method is mostly focusing on solving the variational upper bound, while overlooks modeling the missing patterns. Suggest the authors could add some experiments with MNAR data. Also it would be better if the authors could add some modeling part on the missing pattern into the loss. For example, let the mask $m$ to depends on $(x, y, z)$. This will make the proposed model more useful in practice.\n\n2. From Proposition 2, we know that the surrogate parameterization can make the optimal solution of CELBO remains under the zero gap case. However, it is nearly impossible to reach the zero gap case in reality since it is unlikely to select variational distributions (i.e. $q(\\cdot)$'s) to perfectly estimate the model posterior distributions. What about the \"sub-optimal\" cases? Is the CELBO-SP optimal solution close to the CELBO optimal solution in a small gap (but not zero gap) case? I understand this might not be easy, but it will be better if the authors could add some theoretical analysis on this.\n\n3. The performance metrics in Table 1 is not showing the proposed methods can outperform the baselines with big gaps, meaning that the proposed methods is not much better compared to previous approaches empirically. However, I think there are many ways that the authors could try to improve the performances. For example, the authors could try a different divergence function, a better way to add the regularization, etc.", "summary_of_the_paper": "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.\n\nReferences:\n\n[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.", "summary_of_the_review": "The author propose an interesting discriminative learning approach with generative modeling to solve the missing data modeling problem, by extending the traditional variational lower bound (ELBO), with a novel and stable upper bound that can be estimated without bias with Monte Carlo estimation. It is better if the author can study more on the missing data pattern (both empirically and theoretically) and the optimal solution preservation (under sub-optimal cases). Also, the empirical performance still has some space to improve.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "lwnRX4lZh3", "reviewer": "Reviewer_unNG", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "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": "", "main_review": "(+) pros / (-) cons\n-------------------\n\n(+) Predictive modelling over inputs with missing features is an important problem arising often in many application domains. This paper contributes this somewhat underexplored field.\n\n(-) The rather low documented performance benefits over simpler baselines do not justify the use of the complex model (combining 5? networks) proposed here as opposed to simpler VAE or CVAE variants.\n\n(+) The method and the various bounds introduced are mathematically intriguing, well motivated and potentially useful in follow-up research, however, ... \n(-) the paper is difficult to  follow and at places the reader is left guessing what the authors meant. This should be improved. Concretely\n1. last para of section 2 - the optimization of negative L \"is relatively difficult\". Why? What makes it difficult?\n2. last para of section 2 - \"... have been no equivalent of VI ... \" What about CUBO and its variants pick up from in your work? These have some specific flaws for which they do not qualify here?\n3. before equation 3 - \"exponential divergence\". You mean the Bregman exponential divergence? A citation to help the reader?\n4. $p(u, z | \\theta)$ in equations (11) and (12) seem to use the same parameters $\\theta$ though for (11) $u = (y, \\tilde{x})$ and for (12) it is $u = \\tilde{x}$. Is this in practice the same network with two outputs?\n5. But then in equation (15) these use different $z_{\\theta}$ and $z_{\\psi}$ samples. How is this designed and trained in practice?\n6. page 5 - clarify notation for and explain the gain function; what is the intuition / purpose for it?\n7. Def 1 - effective parameters are those with gradient zero. \".. i.e. the set of parameters inducing tight variational approximation.\" How zero gradient achieves this in a complex non-convex problem, i.e. can't this be a local non-tight extremum?\n8. page 7 DVAE/ DVAE* - you say these are MNAR and MCAR model variants as in Collier at al. 2020. Can you clarify how these translate into your rather more complex model formulation and what specifically changes in the loss (especially the EUBO part)?\n\nFurther questions for clarification/discussion\n-----------------------------------------------\n\n1) Why do you condition $y$ and $x$ on $m$ in equation (1). These are the complete $x$ data so they should not depend on the masking so that $p(y, x | m) = p(y, x)$. Or is this not true? Or is it the $y$ that depends on $m$? Or is it rather the $m$ which depends on $x$? (As in some values being more likely to be masked?)\n\n2) You introduce two latent variable models in equations (7) and (8). My understanding is that the latent $z$ is shared as (8) is just a marginalization of (7) over $y$. You then formulate to approximate posteriors $q(z| y, \\tilde{x})$ and $q(z| \\tilde{x})$ the first learned through ELBO maximization, the 2nd through EUBO minimization. There is currently no link between the two (approximate) posteriors. Would it make sense to somehow link them? (Sorry, I don't know how and may not be possible, or not easily.)\n\n3) You used the Bayes rule to decompose the predictive conditional log likelihood into two terms in equation (4) of which one you are bounding from the bottom (ELBO) and the oher from top (EUBO). What is the effect on the predictive conditional $p(y | x)$? Is it somehow sandwiched or not really due to splitting and modelling the two non-conditional log likelihoods separately? \n\n\nMinor text problems / typos\n---------------------------\n\n1. This first proposition in in page 6 is numbered 2 (not 1) - confusing", "summary_of_the_paper": "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.", "summary_of_the_review": "The paper contributes to practically very important yet relatively little explored area of research - that of predictive modelling with missing data. The proposed method is rather complex, composed of multiple steps adding onto each other to solve a problem arising in the previous steps. These are all well motivated, however, overall the current presentation of the method is difficult to follow and should be improved to help the reader (see main review). Moreover, the documented performance benefits seem to be rather little to justify the use of such a complex method over simpler baseline. These two (lack of clarity, low performance given the complexity of method) are for me the reasons not to consider the paper for this conference. ", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "1: The contributions are neither significant nor novel."}, "extra_text": {}}, {"review_id": "hq2bSm0sYK2", "reviewer": "Reviewer_vurp", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "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": "", "main_review": "**Strengths**\nThe paper is overall clearly written. The issue it tries to solve, discriminative tasks with missing input features, has great impact for a wide range of practical machine learning problems in real life. Technically, the paper has quite some novelty including the creation of a rigorous lower bound to the true objective using recent advances in the variational inference area, and designed an effective surrogate parameterization to stabilize the optimization.\n\n**Weaknesses and Questions**\n1. Sec. 3.1: More detailed explanation of the exponential divergence would be beneficial. Is there a reference for it? What role does $f(\\boldsymbol{u}; \\xi)$ play? If it can be any real-valued function, why was it chosen to be a Gaussian pdf, as shown in the appendix?\n2. Sec. 3.2: Which standard automatic differentiation library was used? The submission mentioned both the reparameterization trick and the REINFORCE trick - which one was actually used in the experiments?\n3. Sec. 3.3: I don't quite understand how this part works. All I can see that in Eq(17) the problematic ratio term is multiplied by $G$, which is always smaller than 1 and non-increasing based on Figure 3. But why $G$ was defined in that math format? What does $\\vee$ mean? Why does $G$ represent the ratio before and after the transform (equation between proposition 2 and 3)?\n4. Sec. 4.1: Missing value processes (MCAR and MNAR). What do they mean? Are they different ways to decide what values are missing in the feature, and thus leading to different versions of a dataset? If so, shouldn't we also add CVAE*, Simple* and MICE*? Could you give more explanation for the last sentence of Section 4 (saying DVAE* is more robust than DVAE)?\n5. Size of the test datasets. Based on Table 3, the datasets are all quite small, ranging from 353 to 10k data points. And we further split these points into training and test, which makes the training sets even smaller. In the appendix it's said the minibatch size is 521 -- what if the entire training set is smaller than 512? How long did the algorithm take to run on YearPred? Is the algorithm able to be easily extended to larger datasets?\n6. How does DIG works compared with other more recent imputation baselines such as MIWAE (Mattei & Frellsen, 2019) and GAIN (Yoon et al., 2018)?", "summary_of_the_paper": "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.", "summary_of_the_review": "Given the strengths of the paper listed above, I would recommend acceptance for this paper, if the authors can figure out a clear feedback for the questions I summarized when reading the paper.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "Beg2TKjfty0", "reviewer": "Reviewer_7Sj3", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "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": "", "main_review": "The ab-initio generative model class proposed by the authors for handling predictions with missing features is convincing. It has the advantage that the involved distributions are simple (factorising), however at the price of introducing latent variables. To learn it discriminatively requires to maximise a difference of concave functions. The first term is lower bounded by ELBO as in VAEs. The second term requires a tractable upper bound. The authors develop a  novel upper bound (starting from alpha-Rényi divergence) that admits a stochastic gradient estimator.  They further introduce a data dependent surrogate reparametrisation in order to achieve an estimator with low variance. The technical part of the paper is concisely written and correct. \n\nThe authors prove that the transformation used for the reparametrisation preserves the effective parameter subset, i.e. the subset of parameter combinations for which the overall bound is tight. This is indeed a desirable property, but is in my view not sufficient. The reason is, that this effective subset can be very small and cover a subset of simple models only. Moreover, there is no guarantee that the respective gap will become small during learning.\n\nThe experimental section first analyses the learning properties of the method in an ablation study. The authors then show competitiveness of their method by comparing it with existing approaches on a subset of tasks taken from the  UCI Machine Learning Repository. The description of the experiments is clear and reproducible. The experiments are however not fully convincing w.r.t. the scalability of the approach. All networks used for the model and bound construction have only one fully connected hidden layer. This seems to be sufficient for the considered tasks from the  UCI repository. However, this would be not sufficient e.g. for image classification tasks where the involved networks are usually deep CNNs.\n\nFurther comments:\n- You mention earlier works  (Ghahramani & Jordan, 1994; Smola et al., 2005), noting that their applicability is restricted to exponential families. Please explain whether these approaches are / are not applicable for the model class analysed in your work. As I understand it, the models p(y,x,z|m) considered by you are exponential families, but of course after marginalising over z, the resulting mixture model p(y,x|m) is not any more. It remains however unclear to me, whether a DCA (difference of convex functions algorithm) for learning p(y|x,z,m) can be somehow generalised for learning p(y|x,m).\n\n- I would suggest to drop the data instance superscript earlier in the text, e.g. starting from subsection 2.2. the latest. This would in my view improve readability and reduce clutter.", "summary_of_the_paper": "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. ", "summary_of_the_review": "The conceptual part, i.e. the model and the proposed learning approach are in my view concise and sufficiently novel. This outweighs the missing scalability analysis in the experimental part. I would however expect the authors to clearly address the raised conceptual questions.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}]

Author rebuttal:

Meta-review: While generative model can be used to input data, this work propose to a novel discriminative learning approach to optimize this data imputation phase by deriving a discriminative version of the traditional variational lower bound (ELBO). The resulting bound can be estimated without bias with Monte Carlo estimation leads to a practical approach, leading to encouraging experimental performances.

The reviewers recognised the novelty and suggest that this approach, given its novelty and wide applicability, could be considered for an oral presentation.
Decision: Accept (Oral)