# VARIATIONAL INFERENCE FOR DISCRIMINATIVE LEARNING WITH GENERATIVE MODELING OF FEA## TURE INCOMPLETION

**Kohei Miyaguchi, Takayuki Katsuki, Akira Koseki & Toshiya Iwamori**
IBM Research – Tokyo
miyaguchi@ibm.com, {kats,akoseki,iwamori}@jp.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-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.

1 INTRODUCTION

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 ∈ R[d][x] is the corrupted version of the complete input x ∈ R[d][x] and m ∈{0, 1}[d][x] 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, ˜xj = (1 _−_ **mj)xj. Given the incomplete features u, we want to predict the outcome y ∈** R[d][y],
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.

2 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 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][]] := **y[i]** R[d][y] _i_ [n] and x[[][n][]] := **x[i]** R[d][x] _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][]] := {u[i] := (˜x[i], m[i])}i∈[n]
instead of x[[][n][]], where ˜x[i] _∈_ R[d][x] is the corrupted version of x[i] and m[i] _∈{0, 1}[d][x]_ is the indicator
of corruption (‘mask vector’) for all i ∈ [n]. Each mask vector m[i] indicates which entry of ˜x[i] is
corrupted; We have ˜x[i]j [= (1][ −] **[m]j[i]** [)][x][i]j [for all][ i][ ∈] [[][n][]][ and all][ j][ ∈] [[][d][x][]][.]

**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][]]|θ) =_


_p(y[i], u[i]|θ),_ (1)
_iY∈[n]_


where the function p(y[i], u[i]|θ) 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,


_Ly[[][n]|[]]u[(][θ][) := ln][ p][(][y][[][n][]][|][u][[][n][]][, θ][) =]_


ln p(y[i]|u[i], θ). (2)
_iX∈[n]_


Now, let Ly[i] _|u[(][θ][) := ln][ p][(][y][i][|][u][i][, θ][)][ denote the][ i][-th summand. Following Bayes’ rule, each][ L][i]y|u[(][θ][)]_
is computed with a difference of log likelihoods

**y** **u[(][θ][) =][ L]y[i]** _,u[(][θ][)][ −L][i]u[(][θ][)][,]_ (3)
_L[i]_ _|_

where Ly[i] _,u[(][θ][) := ln][ p][(][y][i][,][ u][i][|][θ][)][ and][ L][i]u[(][θ][) := ln]_ _p(y[i], u[i]|θ) dy[i]. Since we employ gradient-based_
optimization methods, it suffices to consider (approximate) evaluation of the value and the gradient
R


-----

of the instance-wise objectives Ly[i] _|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.

2.2 DIG WITH LATENT VARIABLE MODELS (DIG-LVM)

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)
Z

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


_dx_

_p(˜xj_ **z, θ)** _,_ (5)
_{_ _|_ _}[1][−][m][j]_
_j=1_

Y


_p(y, u, z|θ) = p(m) · p(z) · p(y|z, θ) ·_


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


_dx_

_p(˜xj_ **z, m, θ)** _._ (6)
_{_ _|_ _}[1][−][m][j]_
_j=1_

Y


_p(y, u, z|θ) = p(m) · p(z) · p(y|z, m, θ) ·_


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

**y,u(θ) = ln** _p(y, u, z_ _θ) dz,_ **u(θ) = ln** _p(u, z_ _θ) dz,_ (7)
_L_ _|_ _L_ _|_
Z Z

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., **u(θ). The optimization of such objective is relatively difficult**
_−L_
since there have been no equivalent of VI in terms of the flexibility and scalability (see Section 5).

3 VARIATIONAL INFERENCE FOR DIG-LVM

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.


-----

3.1 VARIATIONAL LOWER/UPPER BOUNDS

In this section, we discuss an approximation of **y,u(θ) and** **u(θ) in Equation 7 individually. Note**
_L_ _L_
that both are log integrals of the specific form

**v(θ) := ln** _p(v, z_ _θ) dz,_ (8)
_L_ _|_
Z

where v represents either (y, u) or u. Since such integrals constitute the objective L(θ) 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][|][θ][)]_ =: ELBO(v)(θ, φ), (9)
_∼_ _|_ _q(z_ **v, φ)** _L_

 _|_ 

where _φ_ is a variational parameter of a probability density function _q(z|v, φ)_ and
_DKL(z_ **v)(φ** _θ) := Ez_ _q(z_ **v,φ)[ln** _[q]p[(]([z]z[|][v]v[,φ],θ)[)]_ []][ is the KL divergence of the parameters][ φ][ and][ θ][, We call]
_|_ _∥_ _∼_ _|_ _|_

ELBO(v)(θ, φ) as the evidence lower bound (ELBO). Since ELBO is an expectation of a tractable
_L_
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)(θ∥αψ)_

_p(v, z_ _θ)_

= [1] _|_ =: CUBO(v)(θ, ψ), (10)

_α_ [ln][ E][z][∼][q][(][z][|][v][,ψ][)] _q(z_ **v, ψ)** _L_

 _|_ 

where _ψ_ is a variational parameter of a probability density function _q(z_ **v, ψ)** and
1 _|_
_Dα(z|v)(θ∥ψ) :=_ _α−1_ [ln] dz p[α](z|v, θ)q[1][−][α](z|v, ψ) denotes the α-R´enyi divergence of the pa
rameters θ and ψ. Note here, unlike ELBO, CUBO is not unbiasedly approximated because of
R
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

CUBO(v)(θ, ψ) CUBO(v)(θ, ψ) + Ψα( CUBO(v)(θ, ψ) _f_ (v; ξ))
_L_ _≤L_ _L_ _−_

_α_

_p(v, z_ _θ)_

= _[e][−][αf]_ [(][v][;][ξ][)] Ez _q(_ **v,ψ)** _|_ + f (v; ξ)

_α_ _∼_ _·|_ _q(z_ **v, ψ)** _−_ _α[1]_ [=:][ L][EUBO(][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)(θ, ψ) = ln Lv(θ)._


-----

**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:4: _LDraw minibatch ←_ _|B1_ _|_ _i∈B BL[ˆ] ⊂CELBO[i]_ [n](y|u)[(][θ, φ, ψ, ξ][)][ //][ ˆ]LCELBO(y|u) is given by Equation 13. Use Equation 15

instead for variance reduction.

P

5: (θ, φ, ψ, ξ) ← Update((θ, φ, ψ, ξ), ∇L) // Any gradient-based optimization methods can be used.

6: until converge


3.2 OPTIMIZATION ALGORITHM: VDIG

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


_α_
_f_ (u; ξ) + [1] (13)
_−_ _α_ _[,]_



_LˆCELBO(y|u)(θ, φ, ψ, ξ) := ln_ _q[p]([(]z[y]φ[,][ u]y[,],[ z] u[φ], φ[|][θ][)])_ _α_

_|_ _[−]_ [1]


_p(u, zψ_ _θ)_
_|_

_q(zψ_ **u, ψ)e[f]** [(][u][;][ξ][)]
_|_


where zφ and zψ are Monte-Carlo samples drawn from q(z|y, u, φ) and q(z|u, ψ), respectively.
The gradients of [ˆ]CELBO(y **u) is taken with any standard automatic differentiation libraries, using**
_L_ _|_
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[ˆ]CELBO[i]_ (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).

3.3 VARIANCE REDUCTION FOR VDIG WITH SURROGATE PARAMETRIZATION (SP)

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

_p(u, z_ _θ)_
_wθ,ψ,ξ(u, z) :=_ _|_ (14)

_q(z|u, ψ)e[f]_ [(][u][;][ξ][)]

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 ξ[′] formally[2] by

_p(y, u, z_ _θ[′]) :=_ _[G][(][w][θ,ψ,ξ][(][u][,][ z][))]_ _p(y, u, z_ _θ),_ _f_ (u; ξ[′]) := f (u; ξ) ln Z(θ[′]),
_|_ _Z(θ[′])_ _|_ _−_

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, [ˆ]CELBO-SP(y **u) := [ˆ]CELBO(y** **u)** _TSP, is given by_
_L_ _|_ _L_ _|_ _◦_

_LˆCELBO-SP(y|u)(θ, φ, ψ, ξ) = ln_ _[p][(][y][,][ u][,][ z]q[φ]([|]z[θ][)]φ[G]y[(],[w] u[θ,ψ,ξ], φ)_ [(][u][,][ z][φ][))]

_|_

(15)

_−_ _α[1]_ _α_ _[,]_

_[{][w][θ,ψ,ξ][(][u][,][ z][ψ][)][G][(][w][θ,ψ,ξ][(][u][,][ z][ψ][))][}][α][ −]_ _[f]_ [(][u][;][ ξ][) + 1]

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.

3.4 PREDICTION ALGORITHM

Given (θ, φ, ψ, ξ) trained by Algorithm 1, we want to compute the predictive distribution on new
instances given their incomplete features u[new] := (˜x[new], m[new]). Since the conditional density
_pθ(y_ **u) is intractable in general, we approximate it with the Monte-Carlo method. The approxi-**
_|_
mated conditional distribution is given by


_p(y|u,_ _θ[ˆ]) :=_


_p(y|z[s]ψ[, θ][)][,]_ _kpred ≥_ 1, (16)
_s∈X[kpred]_


_kpred_


where z[s]ψ [are independently drawn from][ q][ψ][(][z][|][u][)][,][ s][ ∈] [[][k][pred][]][. 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._

4 EXPERIMENTS

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.

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), 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 posteriorsulting model is not informed with the context ofR **m as VAE and VAE* are, p(z|y, u[3]** _, θbut the objective is (a). Note that the re-_
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 R´enyi_
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.

4.2 EXPERIMENTAL PROCEDURE AND RESULTS

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.


-----

|ELBO+CUBO|Col2|Col3|Col4|
|---|---|---|---|
|||||
|||||
|||||
|||||
|||||

|Col1|ELBO+EUBO+SP|Col3|Col4|
|---|---|---|---|
|||||
|||||
|||||
|||||
|||||


10 20

Time [sec]


10 20

Time [sec]

|Col1|ELB|O+|Col4|Col5|
|---|---|---|---|---|
||||||
||||||
||||||
||||||
||||||
||10 Tim|e|||


Figure 1: Each of three panels corresponds to a different variational approximation. Lines in each
panel denote the training objectives during 5 independent training runs.

|Col1|AQ-CO AQ-NMHC AQ-NOx Boston Diabetes YearPred|Total|
|---|---|---|


|CVAE DVAE DVAE* Simple MICE VAE VAE*|0.41 (0.02) 5.19 (0.14) 5.33 (0.02) 3.01 (0.18) 5.46 (0.03) 3.49 (0.06) 0.45 (0.03) 5.16 (0.10) 5.26 (0.02) 2.91 (0.23) 5.44 (0.08) 3.36 (0.04) 0.44 (0.03) 5.16 (0.11) 5.27 (0.03) 2.92 (0.21) 5.42 (0.07) 3.35 (0.02) 0.49 (0.03) 5.47 (0.11) 5.40 (0.02) 3.13 (0.18) 5.47 (0.04) 3.62 (0.04) 0.46 (0.06) 5.34 (0.23) 5.38 (0.02) 2.94 (0.10) 5.46 (0.03) 3.61 (0.08) 0.47 (0.01) 5.21 (0.08) 5.49 (0.06) 2.95 (0.23) 5.48 (0.10) 3.59 (0.02) 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.82 (0.04) 3.76 (0.05) 3.76 (0.04) 3.93 (0.04) 3.87 (0.05) 3.87 (0.04) 3.83 (0.04)|
|---|---|---|


Table 1: Test per-sample cross entropies of different algorithms for different datasets. The numbers 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.


4.2.1 EFFECT OF NEW VARIATIONAL APPROXIMATION

**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 **y,u(θ). The marginal**
_L_
evidence **u(θ) is computed with the Monte-Carlo approximation for both CUBO (albeit biased)**
_L_
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.


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 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.


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.


-----

0.400

0.375

0.350

0.325

0.300

0.275

0.250

0.225

0.200


VAE* VAE DVAE* DVAE CVAE Simple MICE

Figure 2: The average performance loss when the missing-value ratio is changed from p = .1 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 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.

5 RELATED WORK


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 α-R´enyi divergence (equivalently, χ-divergence). However, as Pradier
et al. (2019) pointed out, the resulting estimate of CUBO is biased and in some cases numerically
unstable[4]. 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 R´enyi 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.

6 CONCLUSION


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.


-----

REFERENCES

Melissa J Azur, Elizabeth A Stuart, Constantine Frangakis, and Philip J Leaf. Multiple imputation
by chained equations: what is it and how does it work? _International journal of methods in_
_psychiatric research, 20(1):40–49, 2011._

T. Bertin-Mahieux. YearPredictionMSD. UCI Machine Learning Repository, 2011.

Mark Collier, Alfredo Nazabal, and Christopher KI Williams. Vaes in the presence of missing data.
_arXiv preprint arXiv:2006.05301, 2020._

S. De Vito, E. Massera, M. Piga, L. Martinotto, and G. Di Francia. On field calibration of an electronic nose for benzene estimation in an urban pollution monitoring scenario. Sensors and Actua_tors B: Chemical, 129(2):750–757, 2008. ISSN 0925-4005. doi: https://doi.org/10.1016/j.snb._
2007.09.060. [URL https://www.sciencedirect.com/science/article/pii/](https://www.sciencedirect.com/science/article/pii/S0925400507007691)
[S0925400507007691.](https://www.sciencedirect.com/science/article/pii/S0925400507007691)

Arthur P Dempster, Nan M Laird, and Donald B Rubin. Maximum likelihood from incomplete data
via the em 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 inference via χ 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 in Neural Information Processing Systems, pp. 120–127, 1994.

Ian 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 Pietro Perona. Combining generative models and fisher kernels
for object recognition. In Tenth IEEE International Conference on 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 I 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 P 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 in 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 (CVPR’06), volume 1, pp. 87–94. IEEE, 2006._


-----

Romain Lopez, Pierre Boyeau, Nir Yosef, Michael I Jordan, and Jeffrey Regier. Decision-making
with auto-encoding variational bayes. arXiv preprint arXiv:2002.07217, 2020.

Pierre-Alexandre 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 P 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 Zwießele, and Neil D Lawrence. Hybrid
discriminative-generative approach with gaussian processes. In Artificial Intelligence and Statis_tics, pp. 47–56, 2014._

Alessandro Perina, Marco Cristani, Umberto Castellani, Vittorio Murino, and Nebojsa Jojic. A
hybrid generative/discriminative 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 optimizing upper bounds of marginal likelihood based on chi-square divergences. In Symposium on
_Advances in Approximate Bayesian Inference, 2019._

Danilo Jimenez Rezende and Shakir Mohamed. Variational inference with normalizing flows. arXiv
_preprint arXiv:1505.05770, 2015._

Donald B Rubin. Inference and missing data. Biometrika, 63(3):581–592, 1976.

Marek Smieja, Łukasz Struski, Jacek Tabor, Bartosz Zieli´[´] nski, and Przemysław Spurek. Processing
of missing data by neural networks. In Advances in Neural Information Processing Systems, pp.
2719–2729, 2018.

Alexander J Smola, SVN Vishwanathan, and Thomas Hofmann. Kernel methods for missing variables. In 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 J Hand. Good methods for coping with missing data in decision
trees. Pattern Recognition Letters, 29(7):950–956, 2008.

Ronald J Williams. Simple statistical gradient-following algorithms for connectionist reinforcement
learning. Machine learning, 8(3-4):229–256, 1992.

Cheng Zhang, Judith B¨utepage, Hedvig Kjellstr¨om, and Stephan Mandt. Advances in variational
inference. IEEE transactions on pattern analysis and machine intelligence, 41(8):2008–2026,
2018.

Juntang Zhuang, Tommy Tang, Yifan Ding, Sekhar Tatikonda, Nicha Dvornek, Xenophon Papademetris, and James S Duncan. Adabelief optimizer: Adapting stepsizes by the belief in observed gradients. arXiv preprint arXiv:2010.07468, 2020.


-----

|w G(w|)|Col3|Col4|Col5|Col6|Col7|
|---|---|---|---|---|---|---|
|wG(|w)||||||
||||||||
||||||||
||||||||
||||||||
||||||||


_w_

_G(w)_

_wG(w)_


Figure 3: The gain function G(w) with α = 2 (the orange line). The green line shows the size
of the stochastic term wθ,ψ,ξ(u, z) after the transformation, whereas the blue dashed line shows the
original size.

A PROOFS

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][(ln][ w][)]

= ln _[p][(][y][,][ u][,][ z][φ][|][θ][)]_
_∇_ _q(zφ_ **y, u, φ)**

_|_ _[−∇][f]_ [(][u][;][ ξ][) +][ H] _[′][(ln][ w][)][∇]_ [ln][ w]

where w := wθ,ψ,ξ(u, z), H(t) := ln G(e[t]) (e[t]G(e[t]))[α] _/α, and H_ _[′](t) denotes the derivative of_
_−_
_H(t). Thus, we have_

_∥g∥≤|H_ _[′](ln w)| ∥∇_ ln w∥ + 3K ≤ 9K

as |H _[′](t)| ≤_ 2 and ∥∇ ln w∥≤ 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
_α-R´enyi divergence dominates the KL divergence for all α > 1._

B EXPERIMENTAL CONFIGURATIONS

**Computing infrastructure.** The computing infrastructure used in the experiments is summarized
in Table 2.


-----

|CPU|RAM GPU PyTorch|
|---|---|


|Intel(R) Core(TM) i7-6700K 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 zero mean and unit empirical variance on the
training split. More precisely, we apply the following transformation

**x˜[i]** (1 **m[i])** (˜x[i] _µˆx)_ _σˆx,_ **y[i]** (y[i] _µˆy)_ _σˆy,_ _i_ [n],
_←_ _−_ _⊙_ _−_ _⊘_ _←_ _−_ _⊘_ _∈_

_ni=1[(1][−][m]j[i]_ [)][x][i]j _ni=1[(1][−][m]j[i]_ [)(][x][i]j _[−]µ[ˆ]x,j_ )[⊙][2]

where ˆµx,j := P _ni=1[(1][−][m][i]j_ [)] and ˆσx,j[⊙][2] [:=] P _ni=1[(1][−][m][i]j_ [)] for j ∈ [dx], and ˆµy :=

_n1_ _ni=1_ _[y][i][ and][ ˆ]σy[⊙]P[2]_ := _n1_ _ni=1[(][y][i][ −]_ _µ[ˆ]y)[⊙][2]. Here,P ⊙_ and ⊘ denotes the elementwise multiplication and division of vectors, respectively. Note that we refill the missing entries of ˜x[i] with
P P
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, θ) ← ([Q][d]j=1[y] _σ[ˆ]j,y[−][1][)][ ·][ 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_ _[,]_ _√1k_ []][, 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 (approx-_
imately) 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.

C DATASETS

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


-----

AQ-CO AQ-NMHC AQ-NOx Boston Diabetes YearPred[∗]

_n_ 6139 731 6174 404 353 10000
_d_ 10 10 10 13 10 90

Table 3: Summary statistics of datasets. ‘∗’ indicates subsampled dataset.

D METHODS COMPARED

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 ∈ R[k], 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**
_≥_ _∈_
R[d][×∗] and an offset vector b ∈ R[d] for t ∈ R[∗]. The Gaussian layer is defined by the composition of
the Gaussian density function and the linear layers,

Gaussd(t **v) :=** _d[t_ Lineard(v), ϵ + ln(1 + Lineard(v))]
_|_ _N_ _|_

where ϵ := 10[−][8] is a small constant and _d[t_ **v, s] denotes the d-product of the probability density**
_N_ _|_
functions of the Gaussian distributions with mean vj and standard deviation sj evaluated at tj,
_j ∈_ [d], for t, v ∈ R[d] and s ∈ R[d]>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, θ) :=** ln Gauss-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,

ˆVAE(θ, φ) := ln p(y, u **z, θ)** ln 10[z; 0, 1] + ln q(z **y, u, φ),**
_L_ _−_ _|_ _−_ _N_ _|_

where

_p(y, u_ **z, θ) := p(m)** Gaussdy (y **v)** Gauss1(˜xj **v),**
_|_ _·_ _|_ _·_ _|_

_j∈[dxY]:mj_ =0

_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.


-----

|Col1|AQ-CO AQ-NMHC AQ-NOx Boston Diabetes YearPred|Total|
|---|---|---|


|CVAE DVAE DVAE* Simple MICE VAE VAE*|0.76 (0.04) 5.62 (0.13) 5.82 (0.05) 3.32 (0.21) 5.56 (0.05) 3.57 (0.10) 0.76 (0.03) 5.53 (0.08) 5.74 (0.03) 3.31 (0.24) 5.53 (0.02) 3.44 (0.02) 0.76 (0.04) 5.54 (0.08) 5.75 (0.03) 3.23 (0.15) 5.50 (0.02) 3.43 (0.02) 0.87 (0.04) 5.87 (0.08) 5.96 (0.07) 3.43 (0.21) 5.57 (0.02) 3.71 (0.03) 0.95 (0.06) 6.07 (0.70) 6.01 (0.07) 3.80 (0.59) 5.61 (0.06) 3.68 (0.02) 0.80 (0.03) 5.53 (0.07) 5.88 (0.03) 3.37 (0.37) 5.57 (0.03) 3.55 (0.03) 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.11 (0.05) 4.05 (0.04) 4.03 (0.03) 4.23 (0.04) 4.35 (0.15) 4.12 (0.06) 4.08 (0.07)|
|---|---|---|


Table 4: The result of the same experiment of Section 4.2.2 except the missing-value ratio is p = 0.5
instead of p = 0.1.

**CVAE.** The objective of CVAE is given according to Sohn et al. (2015),


ˆCVAE(θ, φ) := ln p(y **z, u, θ)** ln 10[z; 0, 1] + ln q(z **y, u, φ),**
_L_ _−_ _|_ _−_ _N_ _|_


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)** Gauss1(˜xj **v),**
_|_ _·_ _|_ _·_ _|_

_j∈[dxY]:mj_ =0

_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)]** _∈_ R[d][y] and
**t[(3)], t[(4)]** _∈_ R[d][x], ω 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.

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-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 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; 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 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


-----

|Approach|Modeling Objective Missingness-Aware Objective Alignment|
|---|---|


|Two-step Generative Discriminative DIG|s(u), p θ(y|s) ln p θ(y|s(u)) -/+ - p (y, u) ln p (y, u) + - θ θ p θ(y|u) ln p θ(y|u) - + p θ(y, u) ln p θ(y|u) + +|
|---|---|


Table 5: Summarized comparison 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 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.

G MODELING THE MAR PROCESS

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 ) and define ˜x := **x˜j : j** . Then, the MAR process is modeled with
_∈O_ _O_ _{_ _∈O}_

_p(y, u, z|θ) = p(˜xO, m) · p(z) · p(y|z, ˜xO, θ) ·_ _{p(˜xj|z, ˜xO, θ)}[1][−][m][j]_ _,_ (17)

_j∈[Ydx]\O_

where p(˜x _, m) will cancel out as well as p(m) in the MNAR and MCAR processes._
_O_

H ADDITIONAL JUSTIFICATION AND LIMITATION OF THE SURROGATE
TRANSFORM

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 [ˆ]CELBO(y **u) to** [ˆ]CELBO-SP(y **u) = [ˆ]CELBO(y** **u)** _TSP and_
discuss the relationship of these objective functions.L _|_ _L_ _|_ _L_ _|_ _◦_

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]p[(]([y]y[|]|[u]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 inde-**
_|_ _−_ _|_ _|_
pendent 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 that wθ,ψ,ξ[α]_ [(][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]).


-----

