# GRADIENT IMPORTANCE LEARNING FOR INCOMPLETE OBSERVATIONS

**Qitong Gao[∗]** **Dong Wang[∗]** **Joshua D. Amason[∗]** **Siyang Yuan[∗]** **Chenyang Tao[∗]**

**Ricardo Henao[∗]** **Majda Hadziahmetovic[∗]** **Lawrence Carin[∗][,][†]** **Miroslav Pajic[∗]**

ABSTRACT

Though recent works have developed methods that can generate estimates (or imputations) of the missing entries in a dataset to facilitate downstream analysis, most
depend on assumptions that may not align with real-world applications and could
suffer from poor performance in subsequent tasks such as classification. This is particularly true if the data have large missingness rates or a small sample size. More
importantly, the imputation error could be propagated into the prediction step that
follows, which may constrain the capabilities of the prediction model. In this work,
we introduce the gradient importance learning (GIL) method to train multilayer
perceptrons (MLPs) and long short-term memories (LSTMs) to directly perform
inference from inputs containing missing values without imputation. Specifically,
we employ reinforcement learning (RL) to adjust the gradients used to train these
models via back-propagation. This allows the model to exploit the underlying information behind missingness patterns. We test the approach on real-world time-series
(i.e., MIMIC-III), tabular data obtained from an eye clinic, and a standard dataset
(i.e., MNIST), where our imputation-free predictions outperform the traditional
_two-step imputation-based predictions using state-of-the-art imputation methods._

1 INTRODUCTION

We consider learning from incomplete datasets. Components of the data could be missing for various
reasons, for instance, missing responses from participants, data loss, and restricted access issues. This
phenomenon is prevalent in the healthcare domain, mostly in the context of electronic health records
(EHRs), which are structured as patient-specific irregular timelines with attributes, e.g., diagnosis,
laboratory tests, vitals, thus resulting in high missingness across patients for any arbitrary time
point. Such missingness introduces difficulties when developing models and performing inference in
real-world applications such as Kam & Kim (2017); Scherpf et al. (2019). Existing works tackle
this problem by either proposing imputation algorithms (IAs) to explicitly produce estimates of the
missing data, or by imposing imputation objectives during inference, e.g., withholding observed
values as being the ground-truth and learning to impute them. However, some of these either require
additional modeling assumptions for the underlying distributions (Bashir & Wei, 2018; Fortuin et al.,
2020), or formatting of the data (Yu et al., 2016; Schnabel et al., 2016). Other applications depend on
domain knowledge for pre-processing and modeling such as Yang et al. (2018); Kam & Kim (2017),
or introduce additional information-based losses (Cao et al., 2018; Lipton et al., 2016), which are
usually intractable with real-world data. Moreover, generative methods (Mattei & Frellsen, 2019;
Yoon et al., 2018) could result in imputations with high variation (e.g., low confidence of the output
distribution), when data have high missingness rates or small sample sizes. Figure 1 illustrates
imputations generated by two state-of-the-art deep generative models on the MNIST dataset with
50%, 70% and 90% of the pixels set as missing (i.e., masked out). It is observed that both approaches
suffer from inaccurate reconstruction of the original digits as the missingness rate increases, which
is manifested as some of the imputed images not being recognizable as digits. Importantly, the
error introduced by the imputation process can be further propagated into downstream inference and
prediction stages.

Imputing the missing data (either explicitly or implicitly) is, in most cases, not necessary, more so,
considering that sometimes where and when the data are missing can be intrinsically informative.

_∗Duke University, USA. Contact: {qitong.gao, miroslav.pajic}@duke.edu._
_†King Abdullah University of Science and Technology, Saudi Arabia._
[Code available at https://github.com/gaoqitong/gradient-importance-learning.](https://github.com/gaoqitong/gradient-importance-learning)


-----

Consider a scenario in which two patients, A and B, admitted to the intensive care unit (ICU) suffer
from bacterial and viral infections, respectively, and assume that the healthcare provider monitors the
status of the patients by ordering two (slightly) different blood tests periodically, namely, a culture
test/panel specific to bacterial infections and a RT-PCR test/panel specific to viral infections. Hence,
patient A is likely to have much fewer orders (and results) for viral tests. In both cases, the missingness
_patterns are indicative of the underlying condition of both patients. Moreover, such patterns are_
more commonly found in incomplete data caused by missing at random (MAR) or missing not at
random (MNAR) mechanisms, which usually introduce additional difficulties for inference (Little &
Rubin, 2019). Inspired by this, we propose the gradient importance learning (GIL) method, which
facilitates an imputation-free learning framework for incomplete data by simultaneously leveraging
the observed data and the information underlying missingness patterns.

Specifically, we propose to use an importance
matrix to weight the gradients that are used
to update the model parameters in the backpropagation process, after the computational
graph of the model is settled, which allows
the model to exploit the information underlying missingness without imputation. However,
these gradients cannot be tracked by automated
differentiation tools such as Tensorflow (Abadi
et al., 2016). So motivated, we propose to gener- Figure 1: MNIST digits imputed by state-of-the-art
ate the importance matrix using a reinforcement imputation methods MIWAE (Mattei & Frellsen,
learning (RL) policy on-the-fly. This is done by 2019), GAIN (Yoon et al., 2018).
conditioning on the status of training procedure
characterized by the model parameters, inputs and model performance at the current training step.
Concurrently, RL algorithms are used to update the policy by first modeling the back-propagation
process as an RL environment, or a Markov decision process (MDP). Then, by interacting with the
environment, the RL algorithm can learn the optimal policy so the importance matrix can aid the
training of the prediction models to attain desirable performance. Moreover, we also show that our
framework can be augmented with feature learning techniques, e.g., contrastive learning as in Chen
et al. (2020), which can further improve the performance of the proposed imputation-free models.

The technical contributions of this work can be summarized as follows: (i) to the best of our
knowledge, this is the first work that trains deep learning (DL) models to perform accurate imputation_free predictions with missing inputs. This allows models to effectively handle high missing rates, while_
significantly reducing the prediction error compared to existing imputation-prediction frameworks.
(ii) Unlike existing approaches that require additional modeling assumptions or that rely on methods
(or pre-existing knowledge) that intrinsically have advantages over specific types of data, our method
_does not require any assumptions or domain expertise. Consequently, it can consistently achieve top_
performance over a variety of data and under different conditions. (iii) The proposed framework
also facilitates feature learning from incomplete data, as the importance matrix can guide the hidden
layers of the model to capture information underlying missingness patterns during training, which
results in more expressive features, as illustrated in Figure 3.

2 GRADIENT IMPORTANCE LEARNING (GIL)


In this section, we introduce GIL; a method that
facilitates training of imputation-free prediction
models with incomplete data. In Section 2.2, we
show how the gradient descent directions can
be adjusted using the importance matrix A. In
Section 2.3, we introduce an RL framework to
generate the elements of A which can leverage
the information underlying missingness during
training, as illustrated in Figure 2.


Figure 2: Overview of the GIL framework.


2.1 PROBLEM FORMULATION
We consider a dataset containing N observations X = (X1, X2, . . ., XN ) where each Xj can be
represented by a vector xj ∈ R[d] (for tabular data) or a matrix Xj ∈ R[T][j] _[×][d]_ (for time-series) with


-----

_Tj denoting the time horizon of the sequence Xj. Since the focus of this work is principally on_
time-series data, recurrent neural networks will receive significant attention. We also define the set
of missing indicatorsdepending on the dimension of M = (M X1, Mj, and 0’s (or 1’s) correspond to the entries in2, . . ., MN ), where Mj = mj ∈{0, 1}[d] or X Mj that are missingj ∈{0, 1}[T][j] _[×][d]_
(or not missing), respectively. Finally, we assume that eachthus = (y1, y2, . . ., yN ) denotes the set of labels. We define the problem as learning a model Xj is associated with some label yj ∈Y;
_Y_
parameterized by W to directly predict yj by generating ˆyj that maximizes the log-likelihood
_Nj=1_ [log][ p][(]y[ˆ]j|Xj, Mj, W), without imputing the missing values in Xj.
P2.2 GRADIENT IMPORTANCE

Missingness is a common issue found in tabular data and time series, where multi-layered perceptron
(MLP) (Bishop, 2006) and long short-term memory (LSTM) (Hochreiter & Schmidhuber, 1997)
models, respectively, are often considered for predictions. Below we illustrate how the importance
_matrix can be applied to gradients, which are produced by taking regular stochastic gradient descent_
(SGD) steps, toward training MLP models with incomplete tabular inputs. Note that this idea can
be easily extended to sequential inputs with LSTM models and the corresponding details can be
found in Appendix A. First, the definitions of the dataset and missing indicators can be reduced
torespectively. Then we define an MLP with Xtab = (x1, x2, . . ., xN ) and Mtab = ( k hidden layers, form1, m2, . . ., m kN _≥), with2, as follows xj ∈_ R[d] and mj ∈ R[d],

**yˆ = φout(Woutφk(Wk . . . φ2(W2φ1(W1x)))),** (1)

where x ∈ R[d] is the input, Wi and φi are the weight matrix and activation functions for the i-th
hidden layer, respectively, and we have omitted the bias terms for notational convenience. The
first layer can be interpreted as the encoding layer Wenc = W1 that maps inputs to features
_ζ = fenc(x|Wenc), while the rest are the inference layers Winf = {W2, . . ., Wk, Wout} that_
map the features to prediction ˆy = finf (ζ **Winf** ). In the following proposition, we show that the
_|_
gradients of some commonly used loss functions, E(yˆ, y), such as cross entropy or mean squared
errors, w.r.t. W1, can be formulated in the form of an outer product. The proof and details can be
found in Appendix E.

**Proposition 1 Given a MLP and a smooth loss function E(yˆ, y), the gradients of E w.r.t. the encod-**
_ing layer can be formulated as ∂E/∂W1 = ∆1x[⊤], where ∆1 contains the gradients propagated_
_from all the inference layers, x ∈Xtab and y ∈Y._

From this Proposition it can be seen that the gradients that are used to train the encoding layers,
following regular SGD solvers, can be formulated as the outer product between the gradients ∆
propagated from the inference layers and the input x as

(∂E/∂Wenc)SGD = ∆ _· x[⊤],_ (2)

wherethe input. To simplify notation, note that henceforth we use Wenc ∈ R[e][×][d], ∆ _∈_ R[e], x ∈ R[d], e is the dimension of the features, and x to refer to an individual tabular data d is the dimension of
**xj** _tab, regardless of its index._
_∈X_

The corresponding SGD updates for training Wenc, using learning rate α, are Wenc **Wenc**
_α∆_ **x[⊤]. As a result, it is observed from (2) that the j-th (j** [1, d]) column of (∂E/∂ ←Wenc)SGD −

_·_ _∈_
is weighted by the j-th element in x given ∆. However, given that some entries in x could be missing
(according to the missing indicator m) and their values are usually replaced by a placeholder ξ ∈ R,
it may not be ideal do directly back-propagate the gradients in the form of (2), as the features ζ
captured by the encoding layers may not be expressive enough for the inference layers to make
accurate predictions. Instead, we consider introducing the importance matrix A ∈ [0, 1][e][×][d] to adjust
the gradient descent direction as

**Wenc** **Wenc** _α(∂E/∂Wenc)SGD_ **A,** (3)
_←_ _−_ _⊙_

which not only trains the model to capture information from the observed inputs that are most useful
for making accurate predictions, but also to leverage the information underlying the missingness in
the data. The elements of A can be generated using the RL framework introduced in the next section.

Note that above we have omitted all the bias terms to simplify our presentation and note that:
_i) gradients w.r.t. to biases do not conform to the outer product format, and ii) more importantly,_
these gradients do not depend on the inputs – thus in practice, the importance matrix A is only applied


-----

to the gradients of the weights. Moreover, we do not consider convolutional neural network (CNN)
models (LeCun et al., 1999) in this work because of i). Though the proposed framework could still
be applied to CNNs, it may not be as efficient as for MLPs or LSTMs, where the search space for
**A is significantly reduced by taking advantages of the outer product format, as shown in (4) below.**
Importantly, our focus is to address the missingness in tabular data and time-series, in which case
MLPs and LSTMs are appropriate, while the missingness in images is usually related to applications
in other domains such as compressed sensing (Wu et al., 2019) or image super-resolution (Dong et al.,
2015), which we plan to explore in the future.

2.3 RL TO GENERATE IMPORTANCE MATRIX
Now we show how to use RL to generate the importance matrix A. In general, RL aims to learn an
optimal policy that maximizes expected rewards, by interacting with an unknown environment (Sutton
& Barto, 2018). Formally, the RL environment is characterized by a Markov decision process (MDP),
with details introduced below. Each time after the agent chooses an action a = π(s) at state s
following some policy π, the environment responds with the next state s[′], along with an immediate
reward r = R(s, a) given a reward function R(·, ·). The goal is to learn an optimal policy that
maximizes the expected total reward defined as J(π) = Eρπ [[P][T]i=0 _[γ][i][r][(][s][i][, a][i][)]][, where][ T][ is the]_
horizon of the MDP, γ ∈ [0, 1) is the discounting factor, and ρπ = {(s0, a0), (s1, a1), . . . |ai =
_π(si)_ is the sequence of states and actions drawn from the trajectory distribution determined by π.
_}_

In our case, the elements of A are determined on-the-fly following an RL policy πθ, parameterized by
_θ, conditioned on some states that characterize the status of the back-propagation at the current time_
step. The policy πθ is updated, concurrently with the weights Wenc, by an RL agent that interacts
with the back-propagation process modeled as the MDP defined as M = (S, A, P, R, γ), with each
of its elements introduced below.

**State Space S.** The state space characterizes the agent’s knowledge of the environment at each
time step. To constitute the states we consider the 4-tuple s = (x, m, ζ, ˆy) including current training
input x ∈ R[d], the missing indicator m ∈ R[d], the feature ζ, i.e., the outputs of the encoding layer of
MLPs or the hidden states ht of LSTMs, and the predictions ˆy produced by the inference layers.

**Action Space A.** We consider combining the importance matrix A ∈ [0, 1][e][×][d] with the parameter
gradients (∂E/∂Wenc)SGD when updating Wenc following

**Wenc[′]** (4)

_[←]_ **[W][enc]** _[−]_ _[α][(][∂E/∂][W][enc][)][SGD]_ _[⊙]_ **[A][ =][ W][enc]** _[−]_ _[α][∆]_ _[·][ (][x][⊤]_ _[⊙]_ **[a][⊤][);]**
where here the equality follows from (2) and Proposition 1, such that all rows of A can be set to be
equal to the importance a[⊤] [0, 1][d], which is obtained from the policy following a = πθ(s).
_∈_

**Transitions P : S × A →S.** The transition dynamics determines how to transit from a current
state s to the next state s[′] given an action a in the environment. In our case, the pair (x, m) is sampled
from the training dataset at each step, so the transitions x → **x[′]** and m → **m[′]** are determined by
how training samples are selected from the training batch during back-propagation. The update
oftransitions Wenc → ζ =W fenc[′] _enc([is conditioned on the importance]x|Wenc) →_ _ζ_ _[′]_ = fenc(x[′]|Wenc[′] **[ a][)][. The update of][ as captured in (4), which results in the][ W][inf]** _[→]_ **[W]inf[′]** [follows]
from the regular SGD updates. Then, the transition ˆy = f (x|W) → **yˆ[′]** = f (x[′]|W[′]) follows
from the other transitions defined above. Finally, we can define the transition between states as
**s = (x, m, ζ, ˆy) →** **s[′]** = (x[′], m[′], ζ _[′], ˆy[′])._

**Reward Function R.** After the RL agent takes an action a at the state s, an immediate reward
_r = R(s, a) is returned by the environment which provides essential information to update πθ (Silver_
et al., 2014; Lillicrap et al., 2015). We define the reward function as R(s, a) = −E(yˆ, y).

We utilize actor-critic RL (Silver et al., 2014; Lillicrap et al., 2015) to update πθ, which outputs the
importance a that is used to concurrently update Wenc. Specifically, we train the target policy (or
actor) πθ along with the critic Qν : S × A → R, parameterized by ν, by maximizing Jβ(πθ) =
Es _ρβ_ [Qν(s, πθ(s))] which gives an approximation to the expected total reward J(π). Specifically,
_∼_
the trajectories ρβ = {(s0, a0), (s1, a1), . . . |ai = β(si)} collected under the behavioral policy
_β : S →A are used to update θ and ν jointly following_
_ν[′]_ _ν + ανδ_ _νQν(s, a),_ _θ[′]_ _θ + αθ_ _θπθ(s)_ _aQν(s, a)_ **a=πθ(s);** (5)
_←_ _∇_ _←_ _∇_ _∇_ _|_
where β is usually obtained by adding noise to the output of πθ to ensure sufficient exploration of the
state and action space, δ = r + γQν(s[′], πθ(s[′])) _Qν(s, a) is the temporal difference error (residual)_
_−_
in RL, and αθ, αν are the learning rates for θ, ν, respectively.


-----

**Algorithm 1 Gradient Importance Learning (GIL).**
**Input: X**, Y, M, Wenc, Winf, πθ, Qν, αθ, αν, α, E
**Begin:**

1: Initialize Wenc and Winf, actor πθ and critic Qν
2: Sample x from X and obtain the corresponding label y from Y
3: Obtain the feature ζ _fenc(x_ **Wenc) and prediction ˆy = finf** (ζ **Winf** ) from the encoding
_←_ _|_ _|_
and inference layers, respectively

4: s ← (x, m, ζ, ˆy)

5: for iter in 1 : max_iter do
6: Obtain importance from a behavioral policy a = β(s _πθ)_
_|_

7: Train the encoding layer following Wenc[′]

8: Train the inference layers following regular gradient descent,[←] **[W][enc]** _[−]_ _[α][∆]_ _[·][ (] i.e.[x][⊤],[⊙]_ **[a][⊤][)][ as in (4)]**
**Winf[′]**

9: Obtain the prediction following the updated weights[←] **[W][inf][ −]** _[α][(][∂E/∂][W][inf]_ [)][SGD] **ˆy ←** _f_ (x|Wenc[′] _[,][ W]inf[′]_ [)]

10: Obtain the reward r ← _R(s, a)_

11: Get a new sample x[′] from X and obtain the corresponding label y[′] from Y

12: Obtain the feature ζ _[′]_ _←_ _fenc(x[′]|Wenc[′]_ [)][ and prediction][ ˆ]y[′] = finf (ζ _[′]|Winf[′]_ [)][ from the encod-]
ing and inference layers, respectively

13: **s[′]** _←_ (x[′], m[′], ζ _[′], ˆy[′])_

14: Update the actor πθ and critic Qν using (s, a, r, s[′]) following (5)

15: **s** **s[′], Wenc** **Wenc[′]** _[,][ W][inf]_ _inf_
_←_ _←_ _[←]_ **[W][′]**

16: end for

We summarize the GIL approach in Algorithm 1 and the detailed descriptions can be found in
Appendix B. Note that the missing indicator m ∈ R[d] would be a good heuristic to replace the
importance a as it could prevent the gradients produced by the missing dimensions from propagation.
However, it does not train the model to capture the hidden information underlying the missing data,
which results in its performance being dominated by GIL as demonstrated in Section 4.

2.4 EXTENSIONS OF THE FRAMEWORK
The proposed GIL framework uses RL agents to guide the model toward minimizing its objective
_E(ˆy, y), which is usually captured by general prediction losses such as cross entropy or mean_
squared error. Below we use an example to show how our method can be extended to adapt ideas
from related domains, which can augment training of the imputation-free prediction models by
altering the reward function. Recent works in contrastive learning, e.g., Chen et al. (2020), can train
DL models to generate highly expressive features by penalizing (or promoting) the distributional
difference D(ζ [+], ζ _[−]) among the features associated with inputs that are similar to (or different from)_
each other, where ζ [+] denotes the set of features corresponding to one set of homogeneous inputs
_con, and ζ_ _[−]_ are the features generated by data that are significantly different from the ones in _con._
_X_ _X_
However, such methods may not be directly applied to the missing data setting considered in this
work, as their loss D is usually designed toward unsupervised training and might need to be carefully
re-worked in our case. This idea can be adapted to our framework by defining ζ [+] as the features
mapped from the inputs (by the encoding layers Wenc) associated with the same label y and ζ _[−]_
_∈Y_
the features corresponding to some other label y[′] _∈Y. Then we can define the new reward function_
_R(s, a) = −E(ˆy, y) + c · D(ζ_ [+], ζ _[−]), c > 0, which does not require D to be carefully crafted,_
provided that Wenc is not trained by directly propagating gradients from it. In Section 4 it will be
shown in case studies that by simply defining D as the mean squared error between ζ [+] and ζ _[−]_ can
improve the prediction performance.

3 RELATED WORK

**Missing Data Imputation.** Traditional mean/median-filling and carrying-forward imputation methods are still used widely, as they are straightforward to implement and interpret (Honaker & King,
2010). Recently, there have been state-of-the-art imputation algorithms (IAs) proposed to produce
smooth imputations with interpretable uncertainty estimates. Specifically, some adopt Bayesian methods where the observed data are fit to data-generating models, including Gaussian processes (Wilson
et al., 2016; Fortuin & Rätsch, 2019; Fortuin et al., 2018), multivariate imputation by chained equations (MICE) (Azur et al., 2011), random forests (Stekhoven & Bühlmann, 2012), etc., or statistical


-----

optimization methods such as expectation-maximization (EM) (García-Laencina et al., 2010; Bashir
& Wei, 2018). However, they often suffer from limited scalability, require assumptions over the
underlying distribution, or fail to generalize well when used with datasets containing mixed types of
variables, i.e., when discrete and continuous values exist simultaneously (Yoon et al., 2018; Fortuin
et al., 2020). There also exist matrix-completion methods that usually assume the data are static,
_i.e., information does not change over time, and often rely on low-rank assumptions (Wang et al.,_
2006; Yu et al., 2016; Mazumder et al., 2010; Schnabel et al., 2016). More recently, several DL-based
IAs have been proposed following advancements in deep generative models such as deep latent
variable models (DLVMs) and generative adversarial networks (GANs) (Mattei & Frellsen, 2019;
Yoon et al., 2018; Fortuin et al., 2020; Li et al., 2019; Ma et al., 2018). DLVMs are usually trained to
capture the underlying distribution of the data, by maximizing an evidence lower bound (ELBO),
which could introduce high variations to the output distributions (Kingma & Welling, 2013) and
cause poor performance in practice, if the data have high missingness (as for the example shown in
Figure 1). There also exist end-to-end methods that withhold observed values from inputs and impose
reconstruction losses during training of prediction models (Cao et al., 2018; Ipsen et al., 2020). In
addition, data-driven IAs are developed specifically toward medical time series by incorporating prior
domain knowledge (Yang et al., 2018; Scherpf et al., 2019; Calvert et al., 2016; Gao et al., 2021).
On the other hand, a few recent works attempt to address missing inputs through an imputation-free
manner (Morvan et al., 2020; Sportisse et al., 2020), however, they are limited to linear regression.

**Attention.** The importance matrix used in the GIL is somewhat reminiscent of visual attention
mechanisms in the context of image captioning or multi-object recognition (Mnih et al., 2014;
Ba et al., 2014; Xu et al., 2015) as they both require training of the prediction models with RL.
We briefly discuss the connections and distinctions between them, with full details provided in
Appendix C. Visual attentions are commonly used to train CNNs to focus on specific portions of
the inputs that are most helpful for making predictions. They are determined by maximizing an
evidence lower bound (ELBO), which is later proved to be equivalent to the REINFORCE algorithm
in the RL literature (Sutton & Barto, 2018; Mnih et al., 2014). However, these methods cannot be
applied directly to our problem, as they require features to be exclusively associated spatially with
specific parts of the inputs, which is attainable by using convolutional encoders with image inputs
but intractable with the type of data considered in our work. Instead, our approach overcomes this
issue by directly applying importance weights, generated by RL, into the gradient space during
back-propagation. Such issues motivate use of the term importance instead of attention. Lastly, our
method does not require one to formulate the learning objective as an ELBO, and as a result GIL
can adopt any state-of-the-art RL algorithm, without being limited to REINFORCE as in Mnih et al.
(2014); Ba et al. (2014); Xu et al. (2015). Besides, other types of attention mechanisms are also
proposed toward applications in sequence modeling and natural language processing (NLP) such
as Cheng et al. (2016); Vaswani et al. (2017).

4 EXPERIMENTS

We evaluate the performance of the imputation-free prediction models trained by GIL against the
existing imputation-prediction paradigm on both benchmark and real-world datasets, where the
imputation stage employs both state-of-the-art and classic imputation algorithms (IAs), with the
details introduced in Section 4.1. We also consider variations of the GIL method proposed in Section 2
to illustrate its robustness and flexibility. The datasets we use include i) MIMIC-III (Johnson et al.,
2016) that consists of real-world EHRs obtained from intensive care units (ICUs), ii) a de-identified
ophthalmic patient dataset obtained from an eye center in North America, and iii) hand-written digits
MNIST (LeCun & Cortes). We also tested on a smaller scaled ICU time-series from 2012 Physionet
challenge (Silva et al., 2012) and these results can be found in Appendix D.4. Some of the data are
augmented with additional missingness to ensure sufficient missing rates and the datasets we use
cover all types of missing data, i.e., missing complete at random (MCAR), MAR and MNAR (Little
& Rubin, 2019). We found that the proposed method not only outperforms the baselines on all three
datasets under various experimental settings, but also offers better feature representations as shown in
Figure 3. We start with an overview of the baseline methods, in the following sub-section, and then
proceed to present our experimental findings.

4.1 VARIANTS OF GIL AND BASELINES
Our method (GIL) trains MLPs or LSTMs (depending on the type of data) to directly output the
predictions given incomplete inputs without imputation. Following from Section 2.4, we also test on


-----

Table 1: Accuracy and AUC obtained from the MIMIC-III dataset.

|GIL -D -H|GAIN MIWAE GP-VAE BRITS MICE Mean CF kNN MF EM|
|---|---|
|Acc. 93.32 93.09 89.17 Var-l. AUC 96.10 96.79 92.96|90.32 88.71 - - 92.17 88.02 87.32 84.79 75.81 68.20 95.57 94.28 - - 95.97 92.56 91.78 91.86 81.73 75.23|
|Acc. 91.47 91.01 88.25 88.48 86.18 76.50 80.24 90.09 86.41 86.87 85.48 78.11 70.51 Fix-l. AUC 95.29 95.57 92.99 91.94 93.10 81.47 92.13 94.02 91.69 91.98 92.38 84.54 79.97||



binary classification tasks using GIL-D which includes the distributional difference term D(ζ [+], ζ _[−]),_
captured by mean squared errors, into the reward function. For baselines, we start with another
variant of GIL that uses a simple heuristic – the missing indicator m – to replace the importance a
in GIL; we denote it as GIL-H. This helps analyze the advantages of training the models using the
importance obtained from the RL policy learned by GIL, versus a heuristic that simply discards the
gradients produced by the subset of input entries that are missing.

For other baselines, following the imputation-prediction framework, the imputation stage employs
state-of-the-art IAs including MIWAE (Mattei & Frellsen, 2019) which imputes the missing data
by training variational auto-encoders (VAEs) to maximize an ELBO, GP-VAE (Fortuin et al., 2020)
that couples VAEs with Gaussian processes (GPs) to consider the temporal correlation in time-series,
and the generative adversarial network-based method GAIN (Yoon et al., 2018). Further, some
classical imputation methods are also considered, including MICE, missForest (MF), k-nearest
neighbor (kNN), expectation-maximization (EM), mean-imputation (Mean), zero-imputation (Zero)
and carrying-forward (CF). The imputed data from these baselines are fed into the prediction models
that share the same architecture as the ones trained by GIL, for fair comparison. For time-series, we
also compare to BRITS (Cao et al., 2018) which trains bidirectional LSTMs end-to-end by masking
out additional observed data and jointly imposing imputation and classification objectives.

4.2 MIMIC-III
We consider the early prediction of septic shock using the MIMIC-III dataset, following the frameworks provided in recent data-driven works for data pre-processing (Sheetrit et al., 2017; Fleuren
et al., 2020; Khoshnevisan et al., 2020). Specifically, we use 14 commonly-utilized vital signs and
lab results over a 2-hour observation window to predict whether a patient will develop septic shock in
the next 4-hour window. The resulting dataset used for training and testing contains 2,166 sequences
each of which has length between 3 and 125, and the overall missing rate is 71.63%; the missing
data in this dataset can be considered as a mixture of MCAR, MAR and MNAR. More details can be
found in Appendix D.1.

To evaluate performance, we consider the prediction model constituted by a 1024-unit LSTM layer
for encoding followed by a dense output layer for inference. Considering that some of the baseline
methods are not designed to capture the temporal correlations within the sequences which do not
follow a fixed horizon, besides testing with the varied-length sequences (Var-l.) obtained in above,
we also test with a fixed-length version (Fix-l.) where the maximum time step for each sequence is
set to be 8 (see Appendix D.1 for details). The accuracy[1] and AUC for the two tests are summarized
in Table 1. When varied-length sequences are considered, GIL(-D) slightly outperforms GAIN
and MICE while it significantly dominates the other baseline methods. However, when fixedlength sequences are considered, GAIN and MICE’s performance decreases dramatically and are
significantly dominated by that of GIL(-D). This could be caused by omitting the records beyond
the 8-th time step as both models are flexible enough to capture such information from that point onwards, while in contrast, the performance increases for kNN, MF and EM which are in general based
on much simpler models. On the other hand, the models trained by GIL(-D) was not significantly
affected by the information loss. In fact, it emphasizes that applying importance to the gradients
during training can enable the model to capture the information behind missing data. Moreover, the
dramatically increased performance from GIL-H to GIL(-D) in both tests underscores the significance
of training the models using the importance determined by the RL policy learned by GIL(-D), instead
of a pre-defined heuristic. Note that according to a recent survey of septic shock predictions (Fleuren
et al., 2020), the overall maximum AUC attained by existing data-driven methods, which uses domain
knowledge to design models specifically for ICU time-series, is 0.96 and it is comparable to that
attained by our method. Finally, GP-VAE and BRITS require the input sequences to have the same
horizon, so we only report for the Fix-l. case where they under-perform. These are possibly due to

1All accuracy values in this work are obtained using a decision threshold of 0.5.


-----

Table 2: Average Accuracy and AUC obtained from the Ophthalmic dataset over 3 different random

|masks. Subscripts are standard deviations.|Col2|
|---|---|
|M.R. GIL -D -H|GAIN GP-VAE MIWAE MICE Zero|
|Acc. 86.841.43 87.131.09 83.630.83 25% AUC 92.401.33 92.421.44 88.872.16|85.381.09 85.670.83 80.991.8 84.210.72 84.501.8 90.132.09 91.471.02 84.041.47 91.351.35 90.590.32|
|Acc. 35% 83.33 0.72 85.09 2.58 80.41 1.49 AUC 88.49 0.68 90.68 2.36 87.02 3.03|80.41 0.83 80.410.83 79.244.38 80.411.49 80.121.09 88.02 2.15 84.853.29 85.51 3.47 85.87 3.24 88.02 0.97|



Figure 3: t-SNE visualization of the feature space learned by (a) GIL, (b) MIWAE and (c) GAIN on
the MNIST dataset with 90% missing rate.

the high variation in GPVAE’s uncertainty estimation component, and the mechanism of imposing
additional missingness in BRITS[2].

4.3 OPHTHALMIC DATA
We consider identifying diabetic retinopathy (DR), an eye disease caused by diabetes, using the
de-identified data obtained from Duke Eye Center, constituted by a mixture of image feature vectors
and tabular information. Specifically, the data collected from each subject is formulated as a
vector with 4,101 dimensions containing two 2,048-length feature vectors from two retinal imaging
modalities, optic coherence tomography (OCT) and color fundus photography (CFP), followed by a
5-dimensional data characterizing demographic (e.g., age) and diabetic information (e.g., years of
diabetes). Additional details of the dataset and examples of the OCT and CFP images are provided
in the Appendix D.2. At most one of the two types of retinal images could be missing due to
MAR as sometimes the diagnosis can be performed with a single modality (Cui et al., 2021), while
the demographic/diabetic information can be missing due to MCAR. A total of 1,148 subjects are
included in this dataset with a 17% missing rate (M.R.)[3]. We apply additional random masks over
both the input image features (by assuming one or two entire image feature vectors are missing) and
demographic/diabetic information (by assuming some entries are missing) to increase the missing
rate to 25% and 35%, where for each case 3 different random masks are applied.

For this experiment we consider an MLP with 2 hidden layers with 1,000 nodes each for prediction.
The results are shown in Table 2 where the mean and standard deviations (in subscripts) of accuracy
and AUC are reported. Although GIL and GP-VAE result in similar performance in the 25% missing
rate case, GIL still achieves higher accuracy and AUC; these are significantly improved, with
significantly lower standard deviation over GP-VAE, when the missing rate increases to 35%. In
general, GIL-D outperforms all the other methods, with an exception of the higher standard deviation
of accuracy in the 35% M.R. case, which benefits from the flexibility of the framework we proposed.
With increased missingness, the performance of most baselines drops significantly and are close to
the zero-imputation baseline as the image feature vectors occupy most dimensions of the inputs while
imputing with zeros would be a reasonable heuristic.

4.4 MNIST
This work focuses mostly on tabular inputs or time-series; however, we test on MNIST images given
its simple structure, which could be classified using MLPs and more importantly, because results
from these data are easier to interpret. Specifically, we test on the MCAR version of MNIST where a
pre-determined portion of pixels (i.e., out of 50%, 70% and 90%) are masked off uniformly at random
from the original images, and the MAR version, which assumes part of the image is always observable
and the positions of missing pixels are determined by a distribution conditioned on features extracted
from the observable portion following Mattei & Frellsen (2019). For each test, the masks are applied

2The hyper-parameters used to train their models can be found in Appendix D.1
3When calculating missing rates, each feature vector is only counted as a single entry regardless of size.


-----

Table 3: Average Accuracy reported for the MNIST dataset over different missing rates. For each
missing rate, 5 random masks were used resulting in standard deviations shown as subscripts (×10[−][3]).

|M.R. GIL -H|GAIN MIWAE MICE Zero|GP-VAE|
|---|---|---|


|50% 96.29 96.08 0.7 1.2 MCAR 70% 93.35 93.26 0.9 0.9 90% 78.47 78.44 3.9 4.7|96.23 95.46 94.58 95.83 0.7 0.8 0.2 1.7 91.98 93.20 90.53 92.80 3.8 0.8 3 0.8 72.58 77.91 73.48 76.67 9 11.3 5 5.3|96.57 3.1 93.49 3 76.25 2|
|---|---|---|


|MAR - 93.23 93.15 0.3 0.5|92.97 93.18 92.62 82.50 0.8 0.4 1.1 2.9|92.62 1.3|
|---|---|---|



Table 4: Correlation between imputation MSE and prediction accuracy for different MRs.

|Col1|MNIST 50% M.R. 70% M.R. 90% M.R.|Ophthalmic 25% M.R. 35% M.R.|MIMIC Ground-truths Unknown|
|---|---|---|---|


|Pearson’s Coeff. p-value|-0.48 -0.88 -0.86 0.018 <0.01 <0.01|-0.12 -0.54 0.581 <0.01|- -|
|---|---|---|---|



with 5 different random seeds, resulting in the standard deviations reported in Table 3 which also
summarizes obtained results.

We consider using MLPs constituted by two hidden dense layers with 500 nodes each as the prediction
models[4]. Note that in the MCAR setting, GP-VAE achieves slightly higher average accuracy than our
method when 50% and 70% of the pixels missing, because it relies on convolutional encoding layers,
which intrinsically have advantages on this dataset. However, GP-VAE’s performance is associated
with higher standard deviation in both cases and is significantly outperformed by our method when
90% pixels are missing. The models trained by GIL also significantly outperforms the baselines
following the imputation-prediction framework. This suggests that most of the IAs fail to reconstruct
the input data when the missing rate is high (see Figure 1). Moreover, the errors produced during
the imputation are then propagated into the prediction step which results in the inferior performance.
In the MAR setting, the task becomes more challenging as entire rows of pixels could be missing.
Finally, the model trained by GIL outperforms most of the baselines depending on IAs with an
exception of MIWAE, which only slightly outperforms it. Note that GIL-H’s performance is close
to GIL in most of the settings due to the simple structure of the MNIST digits where the missing
indicator m could be a good estimation of the importance a.

**Feature Space Visualization In Figure 3, we use t-distributed stochastic neighbor embedding (t-**
SNE) to visualize the features learned by the prediction model trained by GIL compared to MIWAE
and GAIN. We observe that GIL results in more expressive features than the others as it generates
clusters with clearer boundaries.

4.5 CORRELATION BETWEEN IMPUTATION AND PREDICTION PERFORMANCE
In this section we study if imputation error is correlated with the downstream prediction performance,
under the imputation-prediction framework. Each column of Table 4 is obtained by calculating the
Pearson’s correlation coefficient (Freedman et al., 2007) and its p-value, between the imputation
mean squared error (MSE) and prediction accuracy, across different imputation methods under the
same dataset. It can be observed that the imputation MSE is negatively correlated with prediction
performance for most of the datasets and data with higher M.R. tends to have higher degree of
negative correlation, which indicates that imputation error could be propagated to downstream tasks.
5 CONCLUSION

We have developed the GIL method, training DL models (MLPs and LSTMs) to directly perform
inference with incomplete data, without the need for imputation. Existing methods addressing the
problem of missing data mostly follow the two-step (imputation-prediction) framework. However,
the error produced during imputation can be propagated into subsequent tasks, especially when the
data have high missingness rates or small sample sizes. GIL circumvents these issues by applying
importance weighting to the gradients to leverage the information underlying the missingness patterns
in the data. We have evaluated our method by comparing it to the state-of-the-art baselines using IAs
on two real-world datasets and one benchmark dataset.

4The performance of MIWAE reported here is different from Mattei & Frellsen (2019) because they used a
CNN classifier.


-----

ACKNOWLEDGMENTS

This research was supported in part by the ONR under agreements N00014-17-1-2504 and N00014-201-2745, AFOSR under award number FA9550-19-1-0169, NIH under NIH/NIDDK R01-DK123062
and NIH/NINDS 1R61NS120246 awards, DARPA, as well as the NSF under CNS-1652544 award
and the National AI Institute for Edge Computing Leveraging Next Generation Wireless Networks,
Grant CNS-2112562.

We would also like to thank Ruiyi Zhang (Adobe Research) for advice that leads to improved paper
presentations, and Ge Gao (North Carolina State University) for sharing insights toward processing
the MIMIC-III dataset as well as the use of some imputation baselines.

REFERENCES

Martin Abadi, Paul Barham, Jianmin Chen, Zhifeng Chen, Andy Davis, Jeffrey Dean, Matthieu
Devin, Sanjay Ghemawat, Geoffrey Irving, Michael Isard, Manjunath Kudlur, Josh Levenberg,
Rajat Monga, Sherry Moore, Derek G. Murray, Benoit Steiner, Paul Tucker, Vijay Vasudevan, Pete
Warden, Martin Wicke, Yuan Yu, and Xiaoqiang Zheng. Tensorflow: A system for large-scale
machine learning. 2016.

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

Jimmy Ba, Volodymyr Mnih, and Koray Kavukcuoglu. Multiple object recognition with visual
attention. arXiv preprint arXiv:1412.7755, 2014.

Dzmitry Bahdanau, Kyunghyun Cho, and Yoshua Bengio. Neural machine translation by jointly
learning to align and translate. arXiv preprint arXiv:1409.0473, 2014.

Faraj Bashir and Hua-Liang Wei. Handling missing data in multivariate time series using a vector
autoregressive model-imputation (var-im) algorithm. Neurocomputing, 276:23–30, 2018.

Christopher M Bishop. Pattern recognition and machine learning. Springer, 2006.

Justin Boyan and Andrew W Moore. Learning evaluation functions to improve optimization by local
search. Journal of Machine Learning Research, 1(Nov):77–112, 2000.

Jacob S Calvert, Daniel A Price, Uli K Chettipally, Christopher W Barton, Mitchell D Feldman, Jana L
Hoffman, Melissa Jay, and Ritankar Das. A computational approach to early sepsis detection.
_Computers in biology and medicine, 74:69–73, 2016._

Wei Cao, Dong Wang, Jian Li, Hao Zhou, Lei Li, and Yitan Li. Brits: Bidirectional recurrent
imputation for time series. arXiv preprint arXiv:1805.10572, 2018.

Ting Chen, Simon Kornblith, Mohammad Norouzi, and Geoffrey Hinton. A simple framework for
contrastive learning of visual representations. In International Conference on Machine Learning,
pp. 1597–1607. PMLR, 2020.

Jianpeng Cheng, Li Dong, and Mirella Lapata. Long short-term memory-networks for machine
reading. arXiv preprint arXiv:1601.06733, 2016.

Ying Cui, Ying Zhu, Jay C Wang, Yifan Lu, Rebecca Zeng, Raviv Katz, Filippos Vingopoulos,
Rongrong Le, Inês Laíns, David M Wu, et al. Comparison of widefield swept-source optical
coherence tomography angiography with ultra-widefield colour fundus photography and fluorescein
angiography for detection of lesions in diabetic retinopathy. British Journal of Ophthalmology,
105(4):577–581, 2021.

Aiman Darwiche and Sumitra Mukherjee. Machine learning methods for septic shock prediction. In
_Proceedings of the 2018 International Conference on Artificial Intelligence and Virtual Reality, pp._
104–110, 2018.


-----

Chao Dong, Chen Change Loy, Kaiming He, and Xiaoou Tang. Image super-resolution using deep
convolutional networks. IEEE transactions on pattern analysis and machine intelligence, 38(2):
295–307, 2015.

Lucas M Fleuren, Thomas LT Klausch, Charlotte L Zwager, Linda J Schoonmade, Tingjie Guo,
Luca F Roggeveen, Eleonora L Swart, Armand RJ Girbes, Patrick Thoral, Ari Ercole, et al.
Machine learning for the prediction of sepsis: a systematic review and meta-analysis of diagnostic
test accuracy. Intensive care medicine, 46(3):383–400, 2020.

Vincent Fortuin and Gunnar Rätsch. Deep mean functions for meta-learning in gaussian processes.
_arXiv preprint arXiv:1901.08098, 2019._

Vincent Fortuin, Gideon Dresdner, Heiko Strathmann, and Gunnar Rätsch. Scalable gaussian
processes on discrete domains. arXiv preprint arXiv:1810.10368, 2018.

Vincent Fortuin, Dmitry Baranchuk, Gunnar Rätsch, and Stephan Mandt. Gp-vae: Deep probabilistic
time series imputation. In International Conference on Artificial Intelligence and Statistics, pp.
1651–1661. PMLR, 2020.

David Freedman, Robert Pisani, and Roger Purves. Statistics (international student edition). Pisani,
_R. Purves, 4th edn. WW Norton & Company, New York, 2007._

Qitong Gao, Joshua Amason, Scott Cousins, Miroslav Pajic, and Majda Hadziahmetovic. Automated
identification of referable retinal pathology in teleophthalmology setting. Translational vision
_science & technology, 10(6):30–30, 2021._

Pedro J García-Laencina, José-Luis Sancho-Gómez, and Aníbal R Figueiras-Vidal. Pattern classification with missing data: a review. Neural Computing and Applications, 19(2):263–282,
2010.

Sepp Hochreiter and Jürgen Schmidhuber. Long short-term memory. Neural computation, 9(8):
1735–1780, 1997.

James Honaker and Gary King. What to do about missing values in time-series cross-section data.
_American journal of political science, 54(2):561–581, 2010._

Niels Ipsen, Pierre-Alexandre Mattei, and Jes Frellsen. How to deal with missing data in supervised
deep learning? In ICML Workshop on the Art of Learning with Missing Values (Artemiss), 2020.

Alistair EW Johnson, Tom J Pollard, Lu Shen, H Lehman Li-Wei, Mengling Feng, Mohammad
Ghassemi, Benjamin Moody, Peter Szolovits, Leo Anthony Celi, and Roger G Mark. Mimic-iii, a
freely accessible critical care database. Scientific data, 3(1):1–9, 2016.

[Kaggle. Kaggle diabetic retinopathy detection competition. https://www.kaggle.com/c/](https://www.kaggle.com/c/diabetic-retinopathy-detection/)
[diabetic-retinopathy-detection/, 2015. Accessed: 2021-05-08.](https://www.kaggle.com/c/diabetic-retinopathy-detection/)

Hye Jin Kam and Ha Young Kim. Learning representations for the early detection of sepsis with
deep neural networks. Computers in biology and medicine, 89:248–255, 2017.

Daniel S Kermany, Michael Goldbaum, Wenjia Cai, Carolina CS Valentim, Huiying Liang, Sally L
Baxter, Alex McKeown, Ge Yang, Xiaokang Wu, Fangbing Yan, et al. Identifying medical
diagnoses and treatable diseases by image-based deep learning. Cell, 172(5):1122–1131, 2018.

Farzaneh Khoshnevisan et al. A variational recurrent adversarial multi-source domain adaptation
framework for septic shock early prediction across medical systems. 2020.

Diederik P Kingma and Max Welling. Auto-encoding variational bayes. _arXiv preprint_
_arXiv:1312.6114, 2013._

Yann LeCun and Corinna Cortes. MNIST handwritten digit database.

Yann LeCun, Patrick Haffner, Léon Bottou, and Yoshua Bengio. Object recognition with gradientbased learning. In Shape, contour and grouping in computer vision, pp. 319–345. Springer,
1999.


-----

Ke Li and Jitendra Malik. Learning to optimize. arXiv preprint arXiv:1606.01885, 2016.

Steven Cheng-Xian Li, Bo Jiang, and Benjamin Marlin. Misgan: Learning from incomplete data
with generative adversarial networks. arXiv preprint arXiv:1902.09599, 2019.

Timothy P Lillicrap, Jonathan J Hunt, Alexander Pritzel, Nicolas Heess, Tom Erez, Yuval Tassa,
David Silver, and Daan Wierstra. Continuous control with deep reinforcement learning. arXiv
_preprint arXiv:1509.02971, 2015._

Zachary C Lipton, David Kale, and Randall Wetzel. Directly modeling missing data in sequences
with rnns: Improved classification of clinical time series. In Machine learning for healthcare
_conference, pp. 253–270. PMLR, 2016._

Roderick JA Little and Donald B Rubin. Statistical analysis with missing data, volume 793. John
Wiley & Sons, 2019.

Ran Liu, Joseph L Greenstein, Stephen J Granite, James C Fackler, Melania M Bembea, Sridevi V
Sarma, and Raimond L Winslow. Data-driven discovery of a novel sepsis pre-shock state predicts
impending septic shock in the ICU. Scientific reports, 9(1):1–9, 2019.

Chao Ma, Sebastian Tschiatschek, Konstantina Palla, José Miguel Hernández-Lobato, Sebastian
Nowozin, and Cheng Zhang. Eddi: Efficient dynamic discovery of high-value information with
partial VAE. arXiv preprint arXiv:1809.11142, 2018.

Behrooz Mamandipoor, Mahshid Majd, Monica Moz, and Venet Osmani. Blood lactate concentration
prediction in critical care patients: handling missing values. arXiv preprint arXiv:1910.01473,
2019.

Qingqing Mao, Melissa Jay, Jana L Hoffman, Jacob Calvert, Christopher Barton, David Shimabukuro,
Lisa Shieh, Uli Chettipally, Grant Fletcher, Yaniv Kerem, et al. Multicentre validation of a sepsis
prediction algorithm using only vital sign data in the emergency department, general ward and
ICU. BMJ open, 8(1), 2018.

Pierre-Alexandre Mattei and Jes Frellsen. Miwae: Deep generative modelling and imputation of
incomplete data sets. In International Conference on Machine Learning, pp. 4413–4423. PMLR,
2019.

Rahul Mazumder, Trevor Hastie, and Robert Tibshirani. Spectral regularization algorithms for
learning large incomplete matrices. The Journal of Machine Learning Research, 11:2287–2322,
2010.

Volodymyr Mnih, Nicolas Heess, Alex Graves, et al. Recurrent models of visual attention. Advances
_in Neural Information Processing Systems, 27:2204–2212, 2014._

Marine Le Morvan, Julie Josse, Thomas Moreau, Erwan Scornet, and Gaël Varoquaux. Neumiss
networks: differentiable programming for supervised learning with missing values. 2020.

Matthieu Scherpf, Felix Gräßer, Hagen Malberg, and Sebastian Zaunseder. Predicting sepsis with a
recurrent neural network using the MIMIC III database. Computers in biology and medicine, 113:
103395, 2019.

Tobias Schnabel, Adith Swaminathan, Ashudeep Singh, Navin Chandak, and Thorsten Joachims.
Recommendations as treatments: Debiasing learning and evaluation. In International Conference
_on Machine Learning, pp. 1670–1679. PMLR, 2016._

Eitam Sheetrit, Nir Nissim, Denis Klimov, Lior Fuchs, Yuval Elovici, and Yuval Shahar. Temporal
pattern discovery for accurate sepsis diagnosis in ICU patients. arXiv preprint arXiv:1709.01720,
2017.

Ikaro Silva, George Moody, Daniel J Scott, Leo A Celi, and Roger G Mark. Predicting in-hospital
mortality of ICU patients: The physionet/computing in cardiology challenge 2012. In 2012
_Computing in Cardiology, pp. 245–248. IEEE, 2012._


-----

David Silver, Guy Lever, Nicolas Heess, Thomas Degris, Daan Wierstra, and Martin Riedmiller.
Deterministic policy gradient algorithms. 2014.

Mervyn Singer, Clifford S Deutschman, Christopher Warren Seymour, Manu Shankar-Hari, Djillali
Annane, Michael Bauer, Rinaldo Bellomo, Gordon R Bernard, Jean-Daniel Chiche, Craig M
Coopersmith, et al. The third international consensus definitions for sepsis and septic shock
(sepsis-3). Jama, 315(8):801–810, 2016.

Aude Sportisse, Claire Boyer, Aymeric Dieuleveut, and Julie Josses. Debiasing averaged stochastic
gradient descent to handle missing values. Advances in Neural Information Processing Systems,
33, 2020.

Suhas Sreehari, SV Venkatakrishnan, Katherine L Bouman, Jeffrey P Simmons, Lawrence F Drummy,
and Charles A Bouman. Multi-resolution data fusion for super-resolution electron microscopy. In
_Proceedings of the IEEE conference on computer vision and pattern recognition workshops, pp._
88–96, 2017.

Daniel J Stekhoven and Peter Bühlmann. Missforest—non-parametric missing value imputation for
mixed-type data. Bioinformatics, 28(1):112–118, 2012.

Richard S Sutton and Andrew G Barto. Reinforcement learning: An introduction. MIT press, 2018.

Christian Szegedy, Vincent Vanhoucke, Sergey Ioffe, Jon Shlens, and Zbigniew Wojna. Rethinking
the inception architecture for computer vision. In Proceedings of the IEEE conference on computer
_vision and pattern recognition, pp. 2818–2826, 2016._

Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz
Kaiser, and Illia Polosukhin. Attention is all you need. In Advances in Neural Information
_Processing Systems, pp. 5998–6008, 2017._

Jun Wang, Arjen P De Vries, and Marcel JT Reinders. Unifying user-based and item-based collaborative filtering approaches by similarity fusion. In Proceedings of the 29th annual international
_ACM SIGIR conference on Research and development in information retrieval, pp. 501–508, 2006._

Kaixuan Wei, Angelica Aviles-Rivero, Jingwei Liang, Ying Fu, Carola-Bibiane Schönlieb, and
Hua Huang. Tuning-free plug-and-play proximal algorithm for inverse imaging problems. In
_International Conference on Machine Learning, pp. 10158–10169. PMLR, 2020._

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

Andrew Gordon Wilson, Zhiting Hu, Ruslan Salakhutdinov, and Eric P Xing. Deep kernel learning.
In Artificial intelligence and statistics, pp. 370–378. PMLR, 2016.

Yan Wu, Mihaela Rosca, and Timothy Lillicrap. Deep compressed sensing. In International
_Conference on Machine Learning, pp. 6850–6860. PMLR, 2019._

Kelvin Xu, Jimmy Ba, Ryan Kiros, Kyunghyun Cho, Aaron Courville, Ruslan Salakhudinov, Rich
Zemel, and Yoshua Bengio. Show, attend and tell: Neural image caption generation with visual
attention. In International conference on machine learning, pp. 2048–2057, 2015.

Xi Yang, Yuan Zhang, and Min Chi. Time-aware subgroup matrix decomposition: Imputing missing
data using forecasting events. In 2018 IEEE International Conference on Big Data (Big Data), pp.
1524–1533. IEEE, 2018.

Christopher R Yee, Niven R Narain, Viatcheslav R Akmaev, and Vijetha Vemulapalli. A data-driven
approach to predicting septic shock in the intensive care unit. Biomedical informatics insights, 11:
1178222619885147, 2019.

Jinsung Yoon, James Jordon, and Mihaela Schaar. Gain: Missing data imputation using generative
adversarial nets. In International Conference on Machine Learning, pp. 5689–5698. PMLR, 2018.

Hsiang-Fu Yu, Nikhil Rao, and Inderjit S Dhillon. Temporal regularized matrix factorization for
high-dimensional time series prediction. In Advances in Neural Information Processing Systems,
pp. 847–855, 2016.


-----

A EXTENDING GIL TO LSTMS

**LSTM Model.** We also consider the use of an LSTM as the encoder for sequential inputs with
varied lengths. Specifically, in this case, we can define Xj = (x[⊤]j,1[,][ x][⊤]j,2[, . . .,][ x][⊤]j,Tj [)][⊤] _[∈]_ [R][T][j] _[×][d][ along]_
with Mj = (m[⊤]j,1[,][ m][⊤]j,2[, . . .,][ m][⊤]j,Ti [)][ ∈{][0][,][ 1][}][T][j] _[×][d][ such that each][ x][j,t][ ∈]_ [R][d][ and][ m][j,t][ ∈{][0][,][ 1][}][d][,]
where i [1, N ] and t [1, Tj]. Recall that the forward pass of an LSTM follows
_∈_ _∈_

**ot = σ(Woxj,t + Uoht** 1), it = σ(Wixj,t + Uiht 1), ft = σ(Wf **xj,t + Uf** **ht** 1),
_−_ _−_ _−_ (6)
**gt = tanh(Wgxj,t + Ught** 1), **ct = ft** **ct** 1 + it **gt,** **ht = ot** tanh(ct),
_−_ _⊙_ _−_ _⊙_ _⊙_

whereσ(x) = 1 xt ∈/(1 +R[d] eis the input at time step[−][x]) is the sigmoid function, and t, {Wo, ⊙ Wis the element-wise product. The output layer isi, Wg, Wf _, Uo, Ui, Ug, Uf_ _} are the weights,_
appended after ht to constitute the inference layer, i.e.,

**yˆt = φout(Woutht),** (7)

where φout is the activation function, Winf = Wout are the weights and t ∈ [1, Tj]. The following
proposition shows that given a loss function E(yˆ, y), the gradients of E w.r.t the encoding weights
**Wenc = {Wo, Wi, Wg, Wf** _} can be written in the form of outer products. Note that the gradients_
for (autoregressive) parameters depending on previous hidden states ht−1, i.e.,, {Uo, Ui, Ug, Uf _},_
like the bias terms, cannot be written as outer products; thus, these weights are updated following the
regular SGD, without importance weighting. The proof is provided following the proposition.

**Proposition 2 Given an LSTM as described in (6)-(7) and a smooth loss function E(yˆ, y), the**
_gradients of E w.r.t. Wenc = {Wo, Wi, Wg, Wf_ _} at time t can be written in outer product forms_

_– i.e., it holds that_ _∂[∂E]Wo_ _t_ [=][ ∆]t[o][x]t[⊤][,] _∂∂EWi_ _t_ [=][ ∆]t[i][x][⊤]t _[,]_ _∂∂EWg_ _t_ [=][ ∆]t[g][x]t[⊤][,] _∂∂EWf_ _t_ [=][ ∆]t[f] **[x][⊤]t** _[,][ where]_

**∆[o]t** _[,][ ∆]t[i][,][ ∆]t[g][,][ ∆]t[f]_ _[are the gradients propagated from the inference layers defined in below, and][ x][t]_
_represents the observation obtained at time step t from any sequence Xj_ _in general._
_∈X_

Then it follows that if we definegardless of it index, (4) can be directly used for training the LSTM encoding weights x as the t-th observation from any sequence Xj ∈ R W[T][j] _[×]enc[d], re- =_
**Wo, Wi, Wg, Wf** .
_{_ _}_

Now we prove the Proposition above.

We first define


**∆[o]t** [=][∂E/∂][h][t] (8)

_[⊙]_ [tanh(][c][t][)][ ⊙] **[o][t]** _[⊙]_ [(1][ −] **[o][t][)][,]**

**∆[i]t** [=][∂E/∂][h][t] (9)

_[⊙]_ **[o][t]** _[⊙]_ [(1][ −] [tanh][2][(][c][t][))][ ⊙] **[g][t]** _[⊙]_ **[i][t]** _[⊙]_ [(1][ −] **[i][t][)][,]**

**∆[g]t** [=][∂E/∂][h][t] _[⊙]_ **[o][t]** _[⊙]_ [(1][ −] [tanh][2][(][c][t][))][ ⊙] **[i][t]** _[⊙]_ [(1][ −] **[g]t[2][)][,]** (10)

**∆[f]t** [=][∂E/∂][h][t] (11)

_[⊙]_ **[o][t]** _[⊙]_ [(1][ −] [tanh][2][(][c][t][))][ ⊙] **[c][t][−][1]** _[⊙]_ **[f][t]** _[⊙]_ [(1][ −] **[f][t][)][.]**


Now to prove the proposition we start from the derivative of the smooth loss function E w.r.t. x,
which can be derived as

_∂E_ _∂(Woutht)_ _⊤_

= (12)
_∂ht_ _∂ht_ _·_
 

_∂E_

_φ[′](Woutht)_ (13)
_∂yˆt_ _⊙_

 

_∂E_

=Wout[⊤] _∂yˆt_ _⊙_ _φ[′](Woutht)_ _._ (14)

 


-----

Then the derivatives of E w.r.t. the ot, ct, ft, ct 1, it and gt are
_−_


_∂E_

= _[∂E]_
_∂ot_ _∂ht_


_∂ht_
_∂ot_


= _[∂E]_ tanh(ct) (15)

_∂ht_ _⊙_


_∂E_
= _[∂E]_
_∂ct_ _∂ht_


_∂ht_

_∂ct_


= _[∂E]_ **ot** (1 tanh[2](ct)) (16)

_∂ht_ _⊙_ _⊙_ _−_


_∂E_
= _[∂E]_
_∂ft_ _∂ht_


_∂ht_

_∂ct_


_∂ct_

_∂ft_


= _∂[∂E]ht_ _⊙_ **ot ⊙** (1 − tanh[2](ct)) ⊙ **ct−1** (17)

_∂E_ _∂ht_ _∂ct_

= _[∂E]_
_∂ct−1_ _∂ht_ _∂ct_ _∂ct−1_

= _[∂E]_ **ot** (1 tanh[2](ct)) **ft** (18)

_∂ht_ _⊙_ _⊙_ _−_ _⊙_


_∂E_

= _[∂E]_
_∂it_ _∂ht_


_∂ht_

_∂ct_


_∂ct_

_∂it_


= _[∂E]_ **ot** (1 tanh[2](ct)) **gt** (19)

_∂ht_ _⊙_ _⊙_ _−_ _⊙_


_∂E_

= _[∂E]_
_∂gt_ _∂ht_


_∂ht_

_∂ct_


_∂ct_
_∂gt_


= _[∂E]_ **ot** (1 tanh[2](ct)) **it.** (20)

_∂ht_ _⊙_ _⊙_ _−_ _⊙_

Now, the derivatives of E w.r.t the weights Wo, Wi, Wg and Wf are


_∂E_ _∂ht_ _∂ot_ (21)

_∂Wo_ _t_ [=][ ∂E]∂ ∂Eht _∂ot_ _∂Wo_

= _∂ht_ _⊙_ tanh(ct) ⊙ **ot ⊙** (1 − **ot)** **x[⊤]t** (22)
 

=∆[o]t **[x]t[⊤]** (23)
_∂E_ _∂ht_ _∂ct_ _∂it_ (24)

_∂Wi_ _t_ [=][ ∂E]∂ht _∂ct_ _∂it_ _∂Wi_

_∂E_
= _∂ht_ _⊙_ **ot ⊙** (1 − tanh[2](ct)) ⊙ **gt ⊙** **it ⊙** (1 − **it)** **x[⊤]t** (25)
h i

=∆[i]t[x][⊤]t (26)
_∂E_ _∂ht_ _∂ct_ _∂gt_ (27)

_∂Wg_ _t_ [=][ ∂E]∂ht _∂ct_ _∂gt_ _∂Wg_

_∂E_
= **ot** (1 tanh[2](ct)) **it** (1 **gt[2][)]** **x[⊤]t** (28)

_∂ht_ _⊙_ _⊙_ _−_ _⊙_ _⊙_ _−_

h i

=∆[g]t **[x]t[⊤]** (29)
_∂E_ _∂ht_ _∂ct_ _∂ft_ (30)

_∂Wf_ _t_ [=][ ∂E]∂ht _∂ct_ _∂ft_ _∂Wf_

_∂E_
= _∂ht_ _⊙_ **ot ⊙** (1 − tanh[2](ct)) ⊙ **ct−1 ⊙** **ft ⊙** (1 − **ft)** **x[⊤]t** (31)
h i

=∆[f]t **[x][⊤]t** _[.]_ (32)


-----

Note that since σ[′](·) = σ(·)(1 − _σ(·)), in the transition between (21) and (22) it follows that_

_∂ot_

_∂Wo_ =σ[′](Woxt + Uoht−1)x[⊤]t (33)

= _σ(Woxt + Uoht−1) ⊙_ (1 − _σ(Woxt + Uoht−1))_ **x[⊤]t** (34)
h i

= [ot ⊙ (1 − **ot)] x[⊤]t** _[.]_ (35)

Similarly, the transition between (24) and (25) follows

_∂it_

= it (1 **it).** (36)
_∂Wi_ _⊙_ _−_

At last, the transition between (30) and (31) follows

_∂ft_

= ft (1 **ft).** (37)
_∂Wf_ _⊙_ _−_

Furthemore, since tanh[′](·) = 1 − tanh[2](·), the transition between (27) and (28) follows

_∂gt_

_∂Wg_ = tanh[′](Wgxt + Ught−1)x[⊤]t (38)


=(1 − tanh[2](Wgxt + Ught−1))x[⊤]t (39)

=(1 − **gt[2][)][x]t[⊤][,]** (40)


which concludes the proof.


-----

B DESCRIPTION OF ALGORITHM 1

The algorithm takes as input the training dataset X, the training targets Y, the missing indicators M,
the weights of the encoding layers Wenc, the weights for the inference layers Winf, the actor πθ, the
critic Qν, learning rates _α, αθ, αν_ and training loss function E. Our approach starts by initializing
_{_ _}_
all the parameters Wenc, Winf, πθ, Qν, sampling x with the corresponding m that
_∈X_ _∈M_
will be used for training in the first iteration, obtaining the feature ζ and the prediction ˆy, which
constitute the initial state as s = (x, m, ζ, ˆy). In each iteration, first the importance is generated
from a behavioral policy β that is conditioned on the target policy πθ, such as the noisy exploration
policy proposed in Lillicrap et al. (2015). Then the encoding layer is trained following (4), while the
inference layers are trained following the regular gradient descent. After training, the new prediction
is obtained following the updated weights and its value is assigned to ˆy, which is then used to generate
the reward following the reward function R. Then the training sample x[′] for the next iteration is
sampled and the corresponding m[′], ζ _[′], ˆy[′]_ are obtained to constitute the next state s[′]. Finally, the actor
_πθ and critic Qν are updated following (5). We refer to Silver et al. (2014); Lillicrap et al. (2015) for_
more details on actor-critic RL.


-----

C IMPORTANCE VS ATTENTIONS

We now illustrate the connections and distinctions between the importance in the GIL and visual
attentions Mnih et al. (2014); Ba et al. (2014); Xu et al. (2015), which are commonly used to train
CNNs to focus on specific dimensions of the inputs that are most helpful for making predictions in
the context of image captioning and multi-object recognition. In visual attentions, an input image is
first encoded into a set of vectors I = {l1, . . ., lL} where each li is a feature vector corresponding
to a specific region in the image as they are retrieved from lower convolutional layers. Then, the
attention αt,i for time step t and region i ∈ [1, L] is generated following Bahdanau et al. (2014) as
_αt,i = exp(et,i)/_ _j=1_ [exp (][e][t,k][)][, where][ e][t,i][ =][ f][att][(][l][i][,][ h][t][−][1][)][,][ f][att][ is usually an MLP (or the so-]
called attention model) and ht 1 is the hidden state of a recurrent neural network (RNN) Hochreiter
_−_
& Schmidhuber (1997) used for prediction. Then, a weighted average of the feature vectors[P][L] _i_ _[l][t,i][l][i]_

is fed into the prediction network for further inference where lt,i is a random variable parametrized
by αt,i following p(lt,i = 1 _lj<t,_ ) = αt,i, and lj<t represents the historical values of lt,i for all
_|_ _I_ [P]
_i ∈_ [1, L].

Given the definition of the multinoulli variable lt,i, the model is not smooth and thus cannot be
trained following the regular back-propagation. So the prediction network is trained to maximize the
evidence lower bound (ELBO) which updates the prediction network parameter θ using


_∂_ log p(y ˜li, ) _li_ )
_|_ _I_ + (log p(y _li,_ ) _b)_ _[∂]_ [log][ p][(˜] _|I_

_∂θ_ _|[˜]_ _I_ _−_ _∂θ_

h


_∂J_

_∂θ_ _n_

_[≈]_ [1]


(41)


_i=1_


where [˜]li are the Monte-Carlo samples drawn from p(lt,i|lj<t, I) and b represents the performance
from a baseline E[log p(y|θbaseline)]. It is notable that (41) is equivalent to the REINFOCE algorithm
in RL Williams (1992). Specifically, the search space can be interpreted as an MDP which is usually
used to model the environment in RL problems – i.e., the state space is constituted by I, the RNN
hidden state ht and the historical visitations lj<t; the actions are multinouli variables lt,i; the reward
function is defined as the marginal log-likelihood log p(y _lt,i,_ ); and the policy is characterized by
_|_ _I_
_p(lt,i_ _lj<t,_ ) Xu et al. (2015); Mnih et al. (2014); Ba et al. (2014).
_|_ _I_

However, these methods cannot be applied directly to our problem as it requires the features I
to exclusively associate with specific parts of the inputs spatially, which is attainable by using
convolutional encoders with image inputs but intractable with tabular inputs considered in our case
(see Section 2.1). Instead, our approach overcomes this issue by directly applying importance,
generated by RL, into the gradient space during back-propagation. Moreover, our method does not
require to formulate the learning objective as ELBOs and as a result GIL can adopt any state-of-the-art
RL algorithm – not limited to REINFORCE only as in Mnih et al. (2014); Ba et al. (2014); Xu et al.
(2015).


-----

D EXPERIMENTAL DETAILS AND ADDITIONAL EXPERIMENTS

The case studies are run on a work station with three Nvidia Quadro RTX 6000 GPUs with 24GB
of memory for each. We use Tensorflow to implement the models and training algorithms. To
train the imputation-free prediction models using GIL, we perform a grid search for the model
learning rate α ∈{0.001, 0.0007, 0.0005, 0.0003, 0.0001, 0.00005, 0.00001}, the exponential decay
step for α is selected from {1000, 750, 500} and the exponential decay rate for α is selected from
_{0.95, 0.9, 0.85, 0.8}. The actor πθ and critic Qν in the GIL (i.e., Alg. 1) are trained using deep_
deterministic policy gradient (DDPG) Lillicrap et al. (2015) where the discounting factor γ = 0.99.
We choose the behavioural policy β in GIL as


_πθ(s)_ with probability p1,
missing indicator m with probability p2, (42)
( random action with probability p3.


_β(s) =_


The learning rates of the actor απ and the critic αQ are selected by performing grid search
from {0.0005, 0.0001, 0.00005, 0.00001} and {0.001, 0.0005, 0.0001} respectively. To train the
imputation-free models using GIL-H/GIL-D, we follow the same grid search for α along with its
decay steps and decay rates. From the implementation perspective, we replace the missing entries,
in the inputs x, with a placeholder value before feeding them into the model. This can avoid value
errors thrown by Tensorflow if the input vectors contain NaN’s. However, note these values will not
be used to update the parameters. For the state-of-the-art baselines – GAIN[5], MIWAE[6], GPVAE[7] and
BRITS[8] – we use the implementations published on the Github by the authors. Adam optimizer is
used to train all the prediction models for baselines, or to compute (∂E/∂W)SGD for GIL, GIL-H
and GIL-D. All the models are trained using a batch size of 128. The details of selecting the other
hyper-parameters of each case study can be found in the corresponding sub-section below.

D.1 MIMIC-III

Figure 4: Graphical depiction of the 2-hour observation window and 4-hour early prediction window
(EPW) considered in this case study.

**Dataset Formulation** The MIMIC-III[9] contains EHRs obtained from roughly 40,000 patients who
stayed in the ICUs in the Beth Israel Deaconess Medical Center between 2001 and 2012 (Johnson
et al., 2016) and septic shock is a severe type of sepsis that results in above 40% mortality rate.
Figure 4 illustrates the 2-hour observation window, 4-hour early prediction window (EPW) and the
relation between inputs and targets for predictions. Note that each work related to septic shock
predictions adopts a slightly different strategy in terms of selecting lab/vital attributes and the lengths
of the observation and EPW (Mao et al., 2018; Darwiche & Mukherjee, 2018; Yee et al., 2019;
Sheetrit et al., 2017; Liu et al., 2019; Khoshnevisan et al., 2020; Fleuren et al., 2020). In this work, the
4-hour EPW allows including sufficient number of subjects which can ensure statistical significance
of the results, while the smaller observation window keeps the task challenging. To formulate the
training and testing dataset, we first selected 14 commonly used attributes for predictions as suggested
in previous works (Sheetrit et al., 2017; Fleuren et al., 2020; Khoshnevisan et al., 2020), which are
constituted by vital signs including temperature, respiratory rate, heart rate, systolic blood pressure,

[5https://github.com/jsyoon0823/GAIN](https://github.com/jsyoon0823/GAIN)
[6https://github.com/pamattei/miwae](https://github.com/pamattei/miwae)
[7https://github.com/ratschlab/GP-VAE](https://github.com/ratschlab/GP-VAE)
[8https://github.com/caow13/BRITS](https://github.com/caow13/BRITS)
[9Data acquired from https://mimic.mit.edu. Access to this dataset follows the MIT-CLP license.](https://mimic.mit.edu)


-----

mean arterial pressure, peripheral oxygen saturation (or SpO2), fraction of inspired oxygen (or
FIO2), and lab tests including white blood cell count, serum lactate level, platelet count, creatinine,
bilirubin, bandemia of white blood cells, blood urea nitrogen. To form the septic shock cohort,
we first identify the patients who were diagnosed with sepsis by indexing with the International
Classification of Diseases 9 (ICD-9) code 995.91 and 995.92. Then the patients with septic shock
history are selected using ICD-9 code 785.52 with their septic shock onset time determined following
the third international consensus for sepsis and septic shock (Sepsis-3) standard (Singer et al., 2016).
Finally, we remove the patients whose septic shock onset time was less than 6 hours after admission
to the ICUs since the data is not enough to be formulated as sequences pertaining to the 2-hour
observation and 4-hour prediction window. Consequently, a total of 1,083 septic shock patients are
identified. To form the non-shock (or control) cohort, first 1,083 patients are randomly sampled from
all admissions who have at least 2 hours of records excluding the ones who are in the septic shock
cohort, then a 2-hour time frame is randomly selected to be used as the observation window. All
patients selected following the above procedure are split into 8:2 to formulate the training and testing
datasets, which are referred to as the varied-length (Var-l.) sequences in Section 4.2. In the test set,
the number of subjects associated with positive and negative labels, respectively, are selected to be
equivalent.

Figure 5: Histogram of the sequences lengths of the MIMIC-III dataset (with maximum length
truncated at 20). It can be observed that most of the sequences have length ≤ 8.

**Formulation of Fix-l. Sequences** To formulate the fixed-length (Fix-l.) sequences, we follow
this protocol – if sequence length is less than 8, we pad it with the latest recording available until
the length reaches 8, otherwise we truncate the length to 8 by discarding the data from that point
on-wards. We specifically choose this threshold because most of the sequences have length ≤ 8, as
shown in Figure 5.

**Type of Missingness** The missing data in this dataset can be considered as a mixture of MCAR,
MAR and MNAR. Specifically, recordings from human mistakes or malfunctioned equipment
manifest as MCAR, the dependency across vital signs and lab results give rise to MAR (e.g., specific
lab tests are ordered only if abnormalities found in the related vital readings or other lab results),
and the mismatch among the sampling frequency of some periodically recorded attributes such as
temperature (obtained hourly) and blood test (conducted daily) can be considered as MNAR (Little &
Rubin, 2019; Mamandipoor et al., 2019).

**Training Details** To train the prediction models for both GIL and baselines, we use maximum
training steps of 2,000 and 4,000 for varied-length and fixed-length sequences respectively. For
imputation of the missing values, we train the imputation method MIWAE using the learning rates
_{0.001, 0.0001} with other hyper-parameters provided by the authors in the code. To train GAIN, we_
use the default learning rates provided by the authors to train the generator and discriminator. GAIN
contains one hyper-parameter that balances between the two losses LG and LM defined in Yoon et al.
(2018) which is selected from {0.1, 1, 10, 100}. To train GP-VAE on the Fix-l. case, the set of hyperparameters we use contain latent_dim = {6, 12, 35}, encoder_size = 256, 256, decoder_size =
256, 256, 256, window_size = {3, 4, 6}, beta = {0.2, 0.5, 0.8}, sigma = 1.005, length_scale =


-----

_{2, 4, 7}. To train BRITS we used the recommended impute_weight = 0.3, label_weight = 1.0._
The term D(ζ [+], ζ _[−]) included in the reward function for GIL-D is defined as in the follows. Assume_
that in each training epoch b/2 inputs with label 0 and b/2 inputs with label 1 are sampled from the
dataset, where b is the batch size. We use F [+] = (ζ1[+][, . . ., ζ]b/[+] 2[)][ ∈] [R][b/][2][×][e][ to denote the features]

associated with positive labels (i.e., 1) in the current batch and F _[−]_ = (ζ1[−][, . . ., ζ]b/[−]2[)][ ∈] [R][b/][2][×][e][ are,]
where e is the dimension of each individual feature ζ. Then we define D as

_MSE(F_ [+][0 : b/4], F _[−][0 : b/4]) + MSE(F_ [+][b/4 : b/2], F _[−][b/4 : b/2])_

_−MSE(F_ [+][0 : b/4], F [+][b/4 : b/2]) − _MSE(F_ _[−][0 : b/4], F_ _[−][b/4 : b/2]),_ (43)

where the slicing index follows the syntax from Python, e.g., F [+][0 : b/4] corresponds to the 0-th to
(b/4 − 1)-th rows of F [+].

D.2 OPHTHALMIC DATA

Figure 6: Example of a 3D OCT volume scan and a 2D OCT image slice from the volume scan.

**Dataset Introduction** This dataset is originally constituted by OCT images, CFP images and
patient EHRs including demographic information and test results related to diabetes. A total of 1,148
subjects are included in this dataset. OCT refers to the two- or three-dimensional images capturing
the retinal architectures of the eyes that are scanned by low-coherence lights. The 3D OCT image is
usually called the volume scan and it is constituted by 61 or 101 2D image slices depending on the
spec of the scanner. Examples of 3D and 2D OCT scans are shown in Figure 6. In this experiment
the volume scans we use contain the 61 slices, while we only consider the center slice (31th) 2D
image. The CFP images are the color images captured by a fundus camera showing the condition of
the interior surface of the eye, where an example is shown in Figure 7.

**Dataset Formulation** To formulate the data into tabular inputs, the 2D OCT and CFP retinal
images are first fed in to two inception-v3 (Szegedy et al., 2016) CNNs, pre-trained with similar
images provided by (Kermany et al., 2018) and (Kaggle, 2015) respectively, and the image feature
vectors with dimension of 2,048 each are output by the global average pooling layer. Then the patient
EHR information include age (integer), sex (boolean), length of diabetic history (integer), A1C result
(integer), which measures blood sugar levels, and if insulin has been used (boolean) constitute a
5-dimensional vector that is concatenated to the end of the OCT and CFP feature vectors. We split
all the subjects into a training cohort and a testing cohort following a ratio of 9:1. In the test set,
the number of subjects associated with positive and negative labels, respectively, are selected to be
equivalent. The raw images from the training cohort are augmented through cropping and rotation
before feeding into inception-v3. All the data from the testing cohort are not augmented. Figure 7
illustrates how the training and testing inputs are formulated, as well as the input-output relation of
the prediction models.

**Training Details** To train the prediction models for both GIL and baselines, we use maximum training steps of 2,000. For imputation of the missing values, MIWAE is trained using
the learning rates {0.001, 0.0001} with other hyper-parameters provided by the authors in the


-----

Figure 7: Diagram showing the pipline of this experiment.

code. GAIN is trained using hyper-parameters selected from {0.1, 1, 10, 100}. GP-VAE is trained
by applying the arguments latent_dim = 256, encoder_size = 256, 256, decoder_size =
256, 256, 256, window_size = 3, beta = 0.8, sigma = 1 to its code published in Github
[(https://github.com/ratschlab/GP-VAE). MF, EM and kNN are not included in this](https://github.com/ratschlab/GP-VAE)
experiment due to the high dimensionality of the inputs which makes these algorithms very computational inefficient as imputation results cannot be produced within days.

D.3 MNIST

**Training Details** To train the prediction models for both GIL and baselines, we use maximum training steps of 1,0000. For imputation of the missing values, MIWAE is trained using the learning rates {0.001, 0.0001} with other hyper-parameters provided by the authors in
the code. GAIN is trained using hyper-parameters selected from {0.1, 1, 10, 100}. GP-VAE is
trained by applying the arguments latent_dim = 256, encoder_size = 256, 256, decoder_size =
256, 256, 256, window_size = 3, beta = 0.8, sigma = 1 to its code published in Github
[(https://github.com/ratschlab/GP-VAE). MF, EM and kNN are not included in this](https://github.com/ratschlab/GP-VAE)
experiment due to the high dimensionality of the inputs which makes these algorithms very computational inefficient as imputation results cannot be produced within days.

**Additional Results** We also trained GIL-D on the MCAR version of MNIST dataset and the
performance are shown in Table 5 below.

Table 5: Accuracy for GIL-D on the MCAR version of MNIST dataset. Standard deviations are in
subscripts (×10[−][3]).

Missing Rate 50% 70% 90%

Acc. 96.290.09 93.350.9 78.473.9

D.4 PHYSIONET

Table 6: Accuracy, AUC and average precision (AP) obtained from the Physionet dataset

GIL GIL-H GAIN GP-VAE MIWAE BRITS Mean Zero

Acc. 87.45 87.38 86.23 86.85 87.12 86.20 86.87 87.15
AUC 82.20 82.00 78.02 82.21 83.47 81.30 82.57 81.14
AP 49.19 48.64 39.89 47.05 49.31 43.57 48.30 47.37


-----

Other than MIMIC-III, we also tested on a smaller scaled ICU time-series from 2012 Physionet
challenge Silva et al. (2012) which contain data obtained from 12,000 patients. We use the data
pre-processed and open-sourced by Fortuin et al. (2020).[10] For each patient the values of 35 different
attributes (e.g., blood pressure, temperature) are recorded over a 48-hour window. As a result, the data
for each patient can be formulated into a matrix in R[48][×][35], i.e., the sequence length across patients
are the same. This dataset has an overall 78.5% missing rate and a binary label is assigned to each
patient where 87% of them are 0’s and 13% are 1’s. Therefore we include average precision (AP)
which is calculated from the precision-recall curve as an additional metric to evaluate the performance
toward imbalanced labels. For this dataset, we consider a 1024-unit LSTM layer for encoding and a
dense output layer for inference. It can be observed from Table 6 that although our method achieves
the highest accuracy and 2-nd highest AP, while most methods perform very close to mean- and zeroimputation. This could be caused by the its simple structure, as all the sequences share the same time
horizon, which significantly reduces the difficulties for the classification task. Moreover, the highly
imbalanced labels and the small population result in the performance to be hardly distinguishable
across different methods. And this is the reason that we focus on the MIMIC-III dataset to study the
strengths and shortcomings of the methods.

[10https://github.com/ratschlab/GP-VAE](https://github.com/ratschlab/GP-VAE)


-----

E PROOF FOR PROPOSITION 1

Here we prove a generalized version of Proposition 1 as stated in the following proposition, which
shows that the gradient of the loss function E w.r.t. any layer in the MLP can be written in the outer
product format.

**Proposition 3 Given an MLP as in (1) and a smooth loss function E(yˆ|y), the gradients for any**
_hidden layer Wi, i_ [1, k] can be represented as
_∈_

_∂E_

= ∆iq[⊤]i 1[,] (44)
_∂Wi_ _−_


_where qi_ 1 = φi 1(Wi 1qi 2) is the output from the i 1-th layer with q0 = x, and ∆i =
_−_ _−_ _−_ _−_ _−_
**Wi[⊤]+1[∆][i][+1][ ⊙]** _[φ]i[′]_ [(][W][i][q][i][−][1][)][ with][ ∆][out][ =][ ∂E]∂yˆ _out[(][W][out][q][2][k][)][ and][ ⊙]_ _[denotes the element-wise]_

_(Hadamard) product._ _[⊙]_ _[φ][′]_

Now we prove Proposition 3.

We consider the MLP characterized by

**yˆ = φout(Woutφk(Wk . . . φ2(W2φ1(W1x))))** (45)

where x ∈ R[d] represents the input to the model, Wi is the weight matrix for the i-th hidden layer,
and φi is the activation function of the i-th hidden layer.

We start with deriving the derivatives for the output layer

_∂E_ _∂yˆ_

= _[∂E]_ (46)
_∂Wout_ _∂yˆ_ _∂Wout_


=[ _[∂E]_ _out[(][W][out][q][k][)]][ ·][ q]k[⊤]_ (47)

_∂yˆ_

_[⊙]_ _[φ][′]_

=∆out · q[⊤]k _[,]_ (48)

where qk is the output from the k-th hidden layer.

Now we show the derivative of E w.r.t. the k-th hidden layer


_∂E_ _∂yˆ_

= _[∂E]_
_∂Wk_ _∂yˆ_ _∂qk_


_∂qk_ (49)

_∂Wk_


=Wout[⊤] [[] _[∂E]_ _out[(][W][out][q][k][)]][ ·][ ∂][q][k]_ (50)

_∂yˆ_ _∂Wk_

_[⊙]_ _[φ][′]_

=Wout[⊤] **[∆][out]** (51)

_[·][ ∂]∂W[q][k]k_

=[(Wout[⊤] **[∆][out][)][ ⊙]** _[φ][′]k[(][W][k][q][k][−][1][)]][ ·][ q][⊤]k−1_ (52)

=∆k · q[⊤]k−1 (53)

We now prove Proposition 3 by induction. First assume that for j ∈ [2, k], the derivative of E w.r.t.
the j-th hidden layer is


_∂E_

_∂Wj_ =[(Wj[⊤]+1[∆][j][+1][)][ ⊙] _[φ][′]j[(][W][j][q][j][−][1][)]][ ·][ q][⊤]j−1_ (54)


=∆j **q[⊤]j** 1[.] (55)
_·_ _−_


Moreover, let us assume that


_∂E_

=Wj[⊤]+1[[][ ∂E] _φ[′]j+1[(][W][j][+1][q][j][)]]_ (56)
_∂qj_ _∂qj+1_ _⊙_

=Wj[⊤]+1[∆][j][+1][.] (57)


-----

Now we need to show that for j − 1, it holds that

_∂E_

=[(Wj[⊤][∆][j][)][ ⊙] _[φ]j[′]_ 1[(][W][j][−][1][q][j][−][2][)]][ ·][ q][⊤]j 2[,] (58)
_∂Wj−1_ _−_ _−_

_∂E_

=Wj[⊤][[][ ∂E] _φ[′]j[(][W][j][q][j][−][1][)]]_ (59)
_∂qj−1_ _∂qj_ _⊙_

=Wj[⊤][∆][j][.] (60)

We first prove that (60) holds, i.e.,


_∂E_ _∂qj_

= _[∂E]_ (61)
_∂qj−1_ _∂qj_ _∂qj−1_

_∂qj_
=Wj[⊤]+1[∆][j][+1] (62)

_∂qj−1_

=Wj[⊤][[][W]j[⊤]+1[∆][j][+1] _[⊙]_ _[φ]j[′]_ [(][W][j][q][j][−][1][)]] (63)

=Wj[⊤][∆][j][,] (64)

which is equivalent to (60).

To prove that (58) holds, we start with


_∂E_ = _[∂E]_ _∂qj_ _∂qj−1_ (65)

_∂Wj−1_ _∂qj_ _∂qj−1_ _∂Wj−1_

=Wj[⊤]+1[∆][j][+1] _∂qj_ _∂qj−1_ (66)

_∂qj−1_ _∂Wj−1_

=Wj[⊤][[][W]j[⊤]+1[∆][j][+1] _j[(][W][j][q][j][−][1][)]][ ∂][q][j][−][1]_ (67)

_[⊙]_ _[φ][′]_ _∂Wj_ 1

_−_

=Wj[⊤][∆][j] _∂qj−1_ (68)

_∂Wj−1_

=[(Wj[⊤][∆][j][)][ ⊙] _[φ]j[′]_ 1[(][W][j][−][1][q][j][−][2][)]][ ·][ q]j[⊤] 2[,] (69)
_−_ _−_

which is equivalent to (58).


-----

F ADDITIONAL RELATED WORKS

**RL for Optimization.** This work is also related to approaches that use RL to solve optimization
problems, or to improve existing solvers. For example, Boyan & Moore (2000) proposes a framework
to improve the efficiency of local search algorithms by predicting the search outcomes and biasing
toward directions that provide higher returns, which is approached by estimating value functions of
the MDP corresponding to the optimization problem. Recently, Li & Malik (2016) trains RL agents
to perform gradient descent steps over linear regression and neural network models. Similarly, Wei
et al. (2020) introduces a policy network to automatically determine parameters for the plug-and-play
frameworks (Sreehari et al., 2017) which can solve inverse imaging problems. Note that these
methods require complete inputs and cannot be adapted to the problem we consider.


-----

G ADDITIONAL ABLATION STUDY

In this section, we compare GIL against an ablation baseline that applies the importance directly to
the inputs x, instead of the gradient space. Specifically, suppose a differential parametric function
_hψ(x), parameterized by ψ, is multiplied with x element-wise. Then, as opposed to (4), the gradient_
updates w.r.t. the encoding layers can formulated as

**Wenc[′]** (70)

_[←]_ **[W][enc]** _[−]_ _[α][∆]_ _[·][ (][x][⊤]_ _[⊙]_ _[h][ψ][(][x][)][⊤][)][.]_

Consequently, all model elements (i.e., hψ, Wenc and Winf ) could be trained by gradient descent,
without introducing RL policies.

Table 7: Performance of the ablation model over the ophthalmic and MNIST datasets

|Ophthalmic|MNIST|
|---|---|


|25% M.R. 35% M.R.|70% M.R. 90% M.R.|
|---|---|


|Acc. 84.50 80.99 1.49 1.09 AUC 89.78 87.26 1.68 2.53|93.22 78.3 0.6e−03 1.8e−03 - -|
|---|---|



In Table 7 it shows the performance of the model described above over the ophthalmic dataset, as
well as MNIST digits with 70% and 90% missing rate (M.R.). Specifically, hψ(x) is captured by
a neural network with the same architecture as the policy πθ in GIL. The training and testing are
conducted following the same convention (e.g., random seeds that determine the missing entries and
hyper-parameter search) as introduced in Sections 4.3, 4.4 and Appendix D.2, D.3. It can be observed
that the ablation baseline attained similar performance as to GIL-H and zero imputation. The reason
that leads to these results would be that the gradients for training hψ(x) still follow the outer product
format ∆ _· x[⊤]_ which are left to be accounted for. Specifically, our method is built on top of the
idea, and the experiments in Section 4 also support that if part of the inputs are missing they may not
provide sufficient information to train the prediction models, given the outer product format of the
gradients. Similarly, the gradients for hψ(x) also follow this format and hψ(x) may not be properly
trained directly on incomplete inputs x, which could in general limit the overall performance of the
prediction model.

**Additional rationales for using RL to obtain gradient importance in GIL.** The ablation baseline above is closely related to visual attention (VA) models as discussed in Appendix C, and could
be seen as a preliminary version of them where hψ(x) is used to re-weight elements in the inputs
**x. VA techniques benefit from the fact that features learned by CNNs are spatially correlated to the**
inputs, so the attentions could be directly applied to learned features. To achieve this, VA formulates
the objective as an ELBO which is shown equivalent to the objective of a basic policy gradient RL
algorithm, REINFORCE (Mnih et al., 2014). However, it was not clear how such VA methods could
be adapted to the problem we consider, as the features learned by MLPs or LSTMs are not spatially
correlated with inputs.

**More on importance versus attentions.** The other baseline, BRITS (Cao et al., 2018), uses
bidirectional RNN to process time-series, with attentions applied to the hidden states of RNNs,
followed by optimizing over imputation and prediction objectives jointly during training. BRITS
is considered as the state-of-the-art method that could predict with incomplete time-series inputs
_end-to-end. However, we showed that BRITS is outperformed by our method on the MIMIC dataset;_
see Table 1. The reason would be that BRITS requires masking off part of the observed inputs to
constitute the imputation objective; thus, the information provided for model training is even more
limited given the intrinsic > 70% M.R. of MIMIC.


-----

