Year: 2020

Abstract:
Handling missing data is one of the most fundamental problems in machine learning. Among many approaches, the simplest and most intuitive way is zero imputation, which treats the value of a missing entry simply as zero. However, many studies have experimentally confirmed that zero imputation results in suboptimal performances in training neural networks. Yet, none of the existing work has explained what brings such performance degradations. In this paper, we introduce the variable sparsity problem (VSP), which describes a phenomenon where the output of a predictive model largely varies with respect to the rate of missingness in the given input, and show that it adversarially affects the model performance. We first theoretically analyze this phenomenon and propose a simple yet effective technique to handle missingness, which we refer to as Sparsity Normalization (SN), that directly targets and resolves the VSP. We further experimentally validate SN on diverse benchmark datasets, to show that debiasing the effect of input-level sparsity improves the performance and stabilizes the training of neural networks.

text: Handling missing data is one of the most fundamental problems in machine learning. Among many approaches, the simplest and most intuitive way is zero imputation, which treats the value of a missing entry simply as zero. However, many studies have experimentally confirmed that zero imputation results in suboptimal performances in training neural networks. Yet, none of the existing work has explained what brings such performance degradations. In this paper, we introduce the variable sparsity problem (VSP), which describes a phenomenon where the output of a predictive model largely varies with respect to the rate of missingness in the given input, and show that it adversarially affects the model performance. We first theoretically analyze this phenomenon and propose a simple yet effective technique to handle missingness, which we refer to as Sparsity Normalization (SN), that directly targets and resolves the VSP. We further experimentally validate SN on diverse benchmark datasets, to show that debiasing the effect of input-level sparsity improves the performance and stabilizes the training of neural networks.

heading: 1 INTRODUCTION
text: Many real-world datasets often contain data instances whose subset of input features is missing. While various imputing techniques, from imputing using global statistics such as mean, to individually imputing by learning auxiliary models such as GAN, can be applied with their own pros and cons, the most simple and natural way to do this is zero imputation, where we simply treat a missing feature as zero. In neural networks, at first glance, zero imputation can be thought of as a reasonable solution since it simply drops missing input nodes by preventing the weights associated with them from being updated. Some what surprisingly, however, many previous studies have reported that this intuitive approach has an adverse effect on model performances (Hazan et al., 2015; Luo et al., 2018; Śmieja et al., 2018), and none of them has investigated the reasons of such performance degradations.
In this work, we find that zero imputation causes the output of a neural network to largely vary with respect to the number of missing entries in the input. We name this phenomenon Variable Sparsity Problem (VSP), which should be avoided in many real-world tasks. Consider a movie recommender system, for instance. It is not desirable that users get different average of predicted ratings just because they have rated different number of movies (regardless of their actual rating values). One might argue that people with less ratings do not like movies in general and it is natural to give higher predicted values to people with more ratings. This might be partially true for users of some sparsity levels, but it is not a common case uniformly applicable for a wider range of sparsity levels. This can be verified in real collaborative filtering datasets as shown in Figure 1 (upper left corner) where users have a similar average rating for test data regardless of the number of known ratings (see also other two examples in Figure 1). However, in standard neural networks with zero imputation, we observe that the model’s inference correlates with the number of known entries of the data instance as shown in
the second row of Figure 11. It would be fatal in some safety-critical applications such as a medical domain: a patient’s probability of developing disease for example should not be evaluated differently depending on the number of medical tests they received (we do not want our model to predict the probability of death is high just because some patient has been screened a lot!).
In addition, we theoretically analyze the existence of VSP under several circumstances and propose a simple yet effective means to suppress VSP while retaining the intuitive advantages of zero imputation: normalizing with the number of non-zero entries for each data instance. We refer to this regularization as Sparsity Normalization, and show that it effectively deals away with the VSP, resulting in significant improvements in both the performance and the stability of training neural networks.
Our contribution in this paper is threefold:
• To best of our knowledge, we are the first in exploring the adverse effect of zero imputation, both theoretically and empirically.
1Note that, this tendency is very consistent with other test points and is observed throughout the entire learning process (even before the training).
• We identify the cause of adverse effect of zero imputation, which we refer to as variable sparsity problem, and formally describe how this problem actually affects training and inference of neural networks (Section 2). We further provide new perspectives using VSP to understand phenomena that have not been clearly explained or that we have misunderstood (Section 4 and 5).
• We present Sparsity Normalization (SN) and theoretically show that SN can solve the VSP under certain conditions (Section 3). We also experimentally reaffirm that simply applying SN can effectively alleviate or solve the VSP yielding significant performance gains (Section 4).

heading: 2 VARIABLE SPARSITY PROBLEM
text: We formally define the Variable Sparsity Problem (VSP) as follows: a phenomenon in which the expected value of the output layer of a neural network (over the weight and input distributions) depends on the sparsity (the number of zero values) of the input data (Figure 2a). With VSP, the activation values of neural networks could become largely different for exactly the same input instance, depending on the number of zero entries; this makes training more difficult and may mislead the model into incorrect predictions.
While zero imputation is intuitive in the sense that it drops the missing input features, we will show that it causes variable sparsity problem for several example cases. Specifically, we show the VSP under assumptions with increasing generality: (Case 1) where activation function is an identity mapping with no bias, (Case 2) where activation function is an affine function, and (Case 3) where activation function is a non-decreasing convex function such as ReLU (Glorot et al., 2011), leaky ReLU (Maas et al., 2013), ELU (Clevert et al., 2016), or Softplus (Dugas et al., 2001).
Here, we summarize the notation for clarity. For a L-layer deep network with non-linearity σ, we use W i ∈ Rni×ni−1 to denote the weight matrix of i-th layer, bi ∈ Rni to denote the bias, hi ∈ Rni to denote the activation vector. For simplicity, we use h0 ∈ Rn0 and hL ∈ RnL to denote input and output layer, respectively. Then, we have
hi = σ(W ihi−1 + bi), for i = 1, · · · , L.
Our goal in this section is to observe the change in hL as the sparsity of h0 (input x) changes. To simplify the discussion, we consider the following assumption:
Assumption 1. (i) Every coordinate of input vector, h0l , is generated by the element-wise multiplication of two random variables h̃0l and ml where ml is binary mask indicating missing value and h̃0l is a (possibly unobserved) feature value. Here, missing mask ml is MCAR (missing completely at random), with no dependency with other mask variables or their values h̃0. All ml follow some identical distribution with mean µm. (ii) The elements of matrix W i are mutually independent and follow the identical distribution with mean µiw. Similarly, b
i and h̃0 consist of i.i.d. coordinates with mean µib and µx, respectively. (iii) µ i w is not zero uniformly over all i.
(i) assumes the simplest missing mechanism. (ii) is similarly defined in Glorot & Bengio (2010) and He et al. (2015) in studying weight initialization techniques. (iii) may not hold under some initialization strategies, but as the learning progresses, it is very likely to hold.
(Case 1) For simplicity, let us first consider networks without the non-linearity nor the bias term. Theorem 1 shows that the average value of the output layer E[hLl ] is directly proportional to the expectation of the mask vector µm:
Theorem 1. Suppose that activation σ is an identity function and that bil is uniformly fixed as zero under Assumption 1. Then, we have E[hLl ] = ∏L i=1 ni−1µ i wµxµm.
(Case 2) When the activation function is affine but now with a possibly nonzero bias, E[hLl ] is influenced by µm in the following way:
Theorem 2. Suppose that activation σ is an affine function under Assumption 1. Suppose further that fi(x) is defined as σ(ni−1µiwx+ µ i b). Then, E[h L l ] = fL ◦ · · · ◦ f1(µxµm).
(Case 3) Finally, when the activation function is non-linear but non-decreasing and convex, we can show that E[hLl ] is lower-bounded by some quantity involving µm:
Theorem 3. Suppose that σ is a non-decreasing convex function under Assumption 1. Suppose further that fi(x) is defined as σ(ni−1µiwx+ µ i b) and µ i w > 0. Then, E[h L l ] ≥ fL ◦ · · · ◦ f1(µxµm).
If the expected value of the output layer (or the lower bound of it) depends on the level of sparsity/missingness as in Theorem 1-3, even similar data instances may have different output values depending on their sparsity levels, which would hinder fair and correct inference of the model. As shown in Figure 1 (second row), the VSP can easily occur even in practical settings of training neural networks where the above conditions do not hold.

heading: 3 SPARSITY NORMALIZATION
text: Algorithm 1 Sparsity Normalization (SN) Input: Dataset D, constant K. Output: Sparsity Normalized Dataset DSN.
Empty set S = φ for each (h0,m) ∈ D do
h0SN ← K · h0/ ‖m‖1 S ← S ∪ { h0SN }
end for DSN ← S
In this section, we propose a simple yet surprisingly effective method to resolve the VSP. We first revisit (Case 2) to find a way of making expected output independent of input sparsity level since the linearity in activation simplifies the correction. Recalling the notation of h0 = h̃0 m ( represents the element-wise product), we find that simply normalizing via h0SN = (h̃
0 m) ·K1/µm for any fixed constant K1, can debias the dependency on the input sparsity level. We name this simple normalizing technique Sparsity Normalization (SN) and describe it in Algorithm 1. Conceptually, this method scales the size of each input value according to its sparsity level so that the change in output size is less sensitive to the sparsity level (Figure 2b). The formal description on correcting sparsity bias by SN is as follows in this particular case:
Theorem 4. (With Sparsity Normalization) Suppose that activation σ is an affine function under Assumption 1. Suppose further that fi(x) = σ(ni−1µiwx+ µ i b) and replace the input layer using SN, i.e. h0SN = (h̃ 0 m)·K1/µm for any fixed constantK1. Then, we haveE[hLl ] = fL◦· · ·◦f1(µx·K1).
Unlike in Theorem 2, SN in Theorem 4 makes average activation to be independent of µm, which determines the sparsity levels of input. It is not trivial to show the counterpart of (Case 3) using SN since E[σ(x)] = σ(E[x]) does not hold in general. However, we show through extensive experiments in Section 4 that SN is practically effective even in more general cases.
While Theorem 4 assumes that µm is known and fixed across all data instances, we relax this assumption in practice and consider varying µm across data instances. Specifically, by a maximum likelihood principle, we can estimate µm for each instance by ∥∥h0∥∥ 0 /n0 = ‖m‖1/n0. Thus, we have h0SN = K · h0/ ‖m‖1 where K = n0 · K1 (see Algorithm 1)2. In practice, we recommend using K as the average of ‖m‖1 over all instances in the training set. We could encounter the dying ReLU phenomenon (He et al., 2015) if K is too small (e.g., K = 1). Since the hyper-parameter K can bring in a regularization effect via controlling the magnitude of gradient (Salimans & Kingma, 2016), we define K = E(h0,m)∈D[‖m‖1] so that the average scales remain constant before and after the normalization, minimizing such side effects caused by SN.

heading: 4 EXPERIMENTS
text: In this section, we empirically show that VSP occurs in various machine learning tasks and it can be alleviated by SN. In addition, we also show that resolving VSP leads to improved model performance on diverse scenarios.

heading: 4.1 COLLABORATIVE FILTERING (RECOMMENDATION) DATASETS
text: We identify VSP and the effect of SN on several popular benchmark datasets for collaborative filtering with extremely high missing rates. We train an AutoRec (Sedhain et al., 2015) using user vector on Movielens (Harper & Konstan, 2016) 100K dataset for validating VSP and SN. Going back to the first column of Figure 1, the prediction with SN is almost constant regardless of ‖m‖1. Another perceived phenomenon in Figure 1 is that the higher the ‖m‖1, the smaller the variation in the prediction of
2When ‖m‖1 is 0, calculation is impossible. Hence, in this case, the ‖m‖1 is assumed to be 1.
the model with SN. Note that the same tendency has been observed regardless of test instances or datasets. This implies that models with SN yield more calibrated predictions; as more features are known for a particular instance, the variance of prediction for that instance should decrease (since we generated independent masks in Figure 1). It is also worthwhile to note that an AutoRec is a sigmoid-based network and Movielens datasets are known to have no MCAR hypothesis (Wang et al., 2018), in which Assumption 1 does not hold at all.
It is also validated that performance gains can be obtained by putting SN into AutoRec (Sedhain et al., 2015), CF-NADE (Zheng et al., 2016), and CF-UIcA (Du et al., 2018), which are the state-ofthe-art among neural network based collaborative filtering models on several Movielens datasets (see Appendix B for detailed settings). In Table 13, we consider three different sized Movielens datasets. Note that AutoRec and CF-NADE allow two types of models according to data encoding (user- or item-rating vector based) and we consider both types. While we obtain performance improvements with SN in most cases, it is more prominent in user-rating based model in case of AutoRec and CF-NADE.
Furthermore, Table 2 compares our simply modification using SN on AutoRec and CF-UIcA with other state-of-the-art collaborative filtering models beyond neural networks. Unlike experiments of AutoRec in Table 1, which use the same network architectures proposed in original papers, here we could successfully learn more expressive network due to the stability obtained by using SN4. For Movielens 100K and 1M datasets, applying SN yields better or similar performance compared to other state-of-the-art collaborative filtering methods beyond neural network based models. It is important to note that all models outperforming AutoRec with SN on Movielens 10M are ensemble models while AutoRec with SN is a single model and it shows consistently competitive results across all datasets.

heading: 4.2 ELECTRONIC MEDICAL RECORDS (EMR) DATASETS
text: We further test VSP and SN for clinical time-series prediction with two Electronic Medical Records (EMR), namely PhysioNet Challenge 2012 (Silva et al., 2012) and the National Health Insurance Service (NHIS) datasets, which have intrinsic missingness since patients will only receive medical examinations that are considered to be necessary. We identify whether the VSP exists with the PhysioNet Challenge 2012 dataset (Silva et al., 2012). We randomly select one test point and plot in-hospital death probabilities as the number of examinations varies (Second column of Figure 1). Without SN, the in-hospital death probability increases as the number of examinations increases,
3We consider a CF-NADE without weight sharing and re-run experiments for fair comparisons because applying SN with weight sharing is not trivial. We also exclude averaging possible choices because it does not make big differences given unnecessary extra computational costs.
4Because overfitting is less with SN in AutoRec, we use twice the capacity than the existing AutoRec model.
even though there is no such tendency in the dataset statistics. However, SN corrects this bias so that in-hospital death probability is consistent regardless of the number of examinations. We observe a similar tendency for examples from the NHIS dataset as well.
Although SN corrects the VSP in both datasets, we perceive different behaviors in both cases in terms of actual performance changes. While SN significantly outperforms its counterpart without SN on NHIS dataset as shown in Table 3, it just performs similarly on PhysioNet dataset (results and detailed settings are deferred to Appendix C). However, SN is still valuable for its ability to prevent biased predictions in this mission-critical area.
In addition, we compare SN with other missing handling techniques to show the efficacy of SN although our main purpose is to provide a deeper understanding and the corresponding solution about the issue that the zero imputation, which is the simplest and most intuitive way of handling missing data, degrades the performance in training neural networks. Even with its simplicity, SN exhibits better or similar performances compared to other more complex techniques. Detailed descriptions of other missing handling techniques are deferred to Appendix H.

heading: 4.3 SINGLE-CELL RNA SEQUENCE DATASETS
text: Single-cell RNA sequence datasets contain expression levels of specific RNA for each single cells. AutoImpute (Talwar et al., 2018) is one of the state-of-the-art methods that imputes missing data on single-cell RNA sequence datasets. We reproduce their experiments using authors’ official implementation, and follow most of their experimental settings (see Appendix D for details).
As before, we first check whether VSP occurs in AutoImpute model. The third column in Figure 1 shows how the prediction of a AutoImpute model changes as the number of known entries changes. Although the number of RNAs found in the specific cell is less related to cell characteristics (upper right corner in Figure 1), the prediction increases as the number of RNAs found in the cell increases. This tendency is dramatically reduced with SN.
Figure 3 shows how imputation performance changes by being worked with SN to several single cell RNA sequence datasets with respect to the portion of train set (see Appendix D for more results). As we can see in Figure 3, SN significantly increases the imputation performance of AutoImpute model. In particular, the smaller the train data, the better the effect of SN in all datasets consistently. AutoImpute model is a sigmoid-based function and single cell RNA datasets (Talwar et al., 2018) do not have the MCAR hypothesis, unlike Assumption 1. Nevertheless, VSP occurs even here and it can be successfully alleviated by SN with huge performance gain. It implies that SN would work for other neural network based imputation techniques.

heading: 4.4 DROPOUT ON UCI DATASETS
text: While SN primarily targets to fix the VSP in the input layer, it could be also applied to any layers of deep neural networks to resolve VSP. The typical example of having heterogeneous sparsity in hidden layers is when we use dropout (Srivastava et al., 2014), which can be understood as another form of zero imputation but at the hidden layers; with Bernoulli dropout, the variance of the number of zero units across instances is np(1− p) (n: the dimension of hidden layer, p: drop rate). While dropout partially handles this VSP issue by scaling 1/(1− p) in the training phase5, SN can exactly correct VSP of hidden layers by considering individual level sparsity (Note that the scaling of dropout can be viewed as applying SN in an average sense: E [‖m‖1] = n(1− p) and K = n).
5In almost all of the deep learning frameworks such as PyTorch, TensorFlow, Theano and Caffe, this inverted dropout is supported.
Figure 4 shows how the RMSE changes as the drop rate changes with and without SN on three popular UCI regression datasets (Boston Housing, Diabetes, and California Housing)6. As illustrated in Figure 4, the larger the drop rate, the greater the difference of RMSE between with and without SN. To explain this phenomenon, we define the degree of VSP as the inverse of signal-to-noise ratio with respect to the number of active units in hidden layers: √ p/(n(1− p)) (expected number of active units over its standard deviation). As can be seen from the figure, the larger drop rate p is, the more severe degree of VSP is and thus the greater protection by SN against performance degradation.

heading: 4.5 DENSITY ESTIMATION
text: In the current literature of estimating density based on deep models, inputs with missing features are in general not largely considered. However, we still may experience the VSP since the proportion of zeros in the data itself can vary greatly from instance to instance. In this experiment, we apply SN to MADE (Germain et al., 2015), which is one of the modern architectures in neural network-based density estimation. We reproduce binarized MNIST (LeCun, 1998) experiments of Germain et al. (2015) measuring negative log likelihood (the lower the better) of test dataset while increasing the number of masks. Figure 5 illustrates the effect of using SN. Note that MADE uses masks in the hidden layer that are designed to produce proper autoregressive conditional probabilities and variable sparsity arises across hidden nodes. SN can be trivially extended to handle this case as well and the corresponding result is given in the figure denoted as w/SN(all)7. We reaffirm that SN is effective even when MCAR assumption is not established.

heading: 5 RELATED WORKS
text: Missing handling techniques Missing imputation can be understood as a technique to increase the generalization performance by injecting plausible noise into data. Noise injection using global statistics like mean or median values is the simplest way to do this (Lipton et al., 2016; Śmieja et al., 2018). However, it could lead to highly incorrect estimation since they do not take into consideration the characteristics of each data instance (Tresp et al., 1994; Che et al., 2018). To overcome this limitation, researchers have proposed various ways to model individualized noise using autoencoders (Pathak et al., 2016; Gondara & Wang, 2018), or GANs (Yoon et al., 2018; Li et al., 2019). However, those model based imputation techniques have not properly worked for high dimensional datasets with the large number of features and/or extremely high missing rates (Yoon et al., 2018) because excessive noise can ruin the training of neural networks rather increasing generalization performance.
For this reason, in the case of high dimensional datasets such as collaborative filtering or single cell RNA sequences, different methods of handling missing data have been proposed. A line of work simply used zero imputation by minimizing noise level and achieved state-of-the-art performance on their target datasets (Sedhain et al., 2015; Zheng et al., 2016; Talwar et al., 2018). In addition, methods using low-rank matrix factorization have been proposed to reduce the input dimension, but these methods not only cause lots of information loss but also fail to capture non-linearity of the input data (Hazan et al., 2015; Bachman et al., 2017; He et al., 2017). Vinyals et al. (2016); Monti et al. (2017) proposed recurrent neural network (RNN) based methods but computational costs for these methods are outrageous for high dimensional datasets. Also, it is not natural to use RNN-based models for non-sequential datasets.
6The most experimental settings are adopted from Klambauer et al. (2017); Littwin & Wolf (2018)’s UCI experiments (see Appendix E for detail).
7The detailed description of the extension is deferred to Appendix F.
Other forms of Sparsity Normalization We discuss other forms of SN to alleviate VSP, already in use unwittingly due to empirical performance improvements. DropBlock (Ghiasi et al., 2018) compensates activation from dropped features by exactly counting mask vector similar to SN (similar approach discussed in Section 4.4.) It is remarkable that we can find models using SN-like normalization even in handling datasets without missing features. For example, in CBOW model (Mikolov et al., 2013) where the number of words used as an input depends on the position in the sentence, it was later revealed that SN like normalization has a practical performance improvement. As an another example, Kipf & Welling (2017) applied Laplacian normalization which is the standard way of representing a graph in graph theory, can naturally handle heterogeneous node degrees and precisely matches the SN operation. In this paper, we explicitly extend SN, which was limited and unintentionally applied to only few settings, to a model agnostic technique.

heading: 6 CONCLUSION
text: We identified variable sparsity problem (VSP) caused by zero imputation that has not been explicitly studied before. To best of our knowledge, this paper provided the first theoretical analysis on why zero imputation is harmful to inference of neural networks. We showed that variable sparsity problem actually exists in diverse real-world datasets. We also confirmed that theoretically inspired normalizing method, Sparsity Normalization, not only reduces the VSP but also improves the generalization performance and stability of feed-forwarding of neural networks with missing values, even in areas where existing missing imputation techniques do not cover well (e.g., collaborative filtering, single cell RNA datasets).

heading: ACKNOWLEDGMENTS
text: This work was supported by the Institute of Information & Communications Technology Planning & Evaluation (IITP) grants (No.2016-0-00563, No.2017-0-01779, and No.2019-0-01371), the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT) grants (No.2018R1A5A1059921 and No.2019R1C1C1009192), the Samsung Research Funding & Incubation Center of Samsung Electronics via SRFC-IT1702-15, the National IT Industry Promotion Agency grant funded by the Ministry of Science and ICT, and the Ministry of Health and Welfare (NO. S0310-19-1001, Development Project of The Precision Medicine Hospital Information System (P-HIS)).

heading: A PROOFS

heading: A.1 PROOF OF THEOREM 1
text: Proof. From the definition of hl, w1l , h 0 l , h̃ 0 l ,ml, the following equation holds.
E[h1l ] = n0E[w 1 l h 0 1] = n0E[w 1 l h̃ 0 lml]
From the Assumption 1, w1l , h̃ 0 l , and ml are independent of each other. Thus,
E[h1l ] = n0E[w 1 l ]E[h̃ 0 l ]E[ml]
Similarly, the following holds.
E[hil] = ni−1E[w i lh i−1 l ] for i = 1, · · · , L
Since hi−1l and w i l are independent of each other by the Assumption 1 and the definition of h i−1 l , E[hil] = ni−1E[w i l ]E[h i−1 l ]. Therefore,
E[hLl ] = L∏ i=1 ni−1E[w i l ]E[h̃ 0 l ]E[ml] = L∏ i=1 ni−1µ i wµxµm

heading: A.2 PROOF OF THEOREM 2
text: Proof. From the definition of hl, w1l , h 0 l , h̃ 0 l ,ml and the property of an affine function σ(E[·]) = E[σ(·)], the following equation holds.
E[h1l ] = σ ( n0E[w 1 l h 0 l ] + E[b 1 l ] ) = σ ( n0E[w 1 l h̃ 0 lml] + E[b 1 l ] )
From the Assumption 1, w1l , h̃ 0 l , and ml are independent of each other. Thus, E[h1l ] = σ ( n0E[w 1 l ]E[h̃ 0 l ]E[ml] + E[b 1 l ] )
= σ ( n0µ 1 wµxµm + µ 1 b ) = f1(µxµm)
Similarly, the following holds. E[hil] = σ ( ni−1E[w i lh i−1 l ] + E[b i l] ) for i = 1, · · · , L
Since hi−1l and w i l are independent of each other by the Assumption 1 and the definition of h i−1 l ,
E[hil] = σ ( ni−1E[w i l ]E[h i−1 l ] + E[b i l] )
= σ ( ni−1µ i wE[h i−1 l ] + µ i b ) = fi(E[h i l])
Therefore,
E[hLl ] = fL ◦ · · · ◦ f1(µxµm)

heading: A.3 PROOF OF THEOREM 3
text: Proof. From the definition of hl, w1l , h 0 l , h̃ 0 l ,ml and the property of a convex function E[σ(·)] ≥ σ(E[·]), the following equation holds.
E[h1l ] ≥ σ ( n0E[w 1 l h 0 l ] + E[b 1 l ] ) = σ ( n0E[w 1 l h̃ 0 lml] + E[b 1 l ] )
From the Assumption 1, w1l , h̃ 0 l , and ml are independent of each other. Thus, E[h1l ] ≥ σ ( n0E[w 1 l ]E[h̃ 0 l ]E[ml] + E[b 1 l ] )
= σ ( n0µ 1 wµxµm + µ 1 b ) = f1(µxµm)
Similarly, the following holds. E[hil] ≥ σ ( ni−1E[w i lh i−1 l ] + E[b i l] ) for i = 1, · · · , L
Since hi−1l and w i l are independent of each other by the Assumption 1 and the definition of h i−1 l ,
E[hil] ≥ σ ( ni−1E[w i l ]E[h i−1 l ] + E[b i l] )
= σ ( ni−1µ i wE[h i−1 l ] + µ i b ) = fi(E[h i−1 l ])
Since we assume that σ is non-decreasing, we finally get
E[hLl ] ≥ fL ◦ · · · ◦ f1(µxµm)

heading: A.4 PROOF OF THEOREM 4
text: Proof. By theorem 2, E[hLl ] = fL ◦ · · · ◦ f1(E[h0SN]). Since E[h0SN] = E[h0] ·K/µm,
E[hLl ] = fL ◦ · · · ◦ f1(µxµm ·K/µm) = fL ◦ · · · ◦ f1(µx ·K)

heading: B COLLABORATIVE FILTERING (RECOMMENDATION) DATASETS

heading: B.1 DETAILED EXPERIMENTAL SETTINGS OF TABLE 1
text: This subsection describes the experimental settings of detailed collaborative filtering tasks in Section 4. As already mentioned, we follow the settings of AutoRec, CF-NADE, and CF-UIcA as much as possible. We perform each experiments by 5 times, and report mean and 95% confidence intervals. We randomly select 10% of the ratings of each datasets for the test set (Harper & Konstan, 2016). As the dataset is too small for Movielens 100K and 1M datasets, the confidence interval tends to be large by changing the dataset split. Hence, the same dataset split is used in each 5 experiments in Table 1.
AutoRec (Sedhain et al., 2015) We use two layer AutoRec model with 500 hidden units. For fair comparisons, we tune the hyper-parameters for weight decay in all experiments to have only one significant digit, and use a learning rate of 10−3 except on Movielens 10M where we use a learning rate of 10−4. We use full batch on Movielens 100K and 1M, mini-batch (1000) on Movielens 10M. Besides, we use Adam optimizer instead of Resilient Propagation (RProp) unlike the AutoRec paper. The RProp optimizer shows fast convergence but can only be used in full batch scenario. It is not possible for 12GB of GPU memory to use full batch in training a large dataset such as Movielens 10M. Thus, we decide to use Adam optimizer rather than RProp optimizer. Fortunately, although the optimizer is changed to Adam, the prediction performance is not degraded in most cases. The experimental results of comparing both optimizers are summarized in Table 4.
CF-NADE (Zheng et al., 2016) We use two layer CF-NADE model with 500 hidden units. For fair comparisons, we tune the hyper-parameters for weight decay in all experiments to have only one significant digit, and use a learning rate8 of 0.001. Also, we use mini-batch (512) just following CF-NADE. Although the CF-NADE used weight sharing and averaging possible choices in addition to weight decay, we report the results without weight sharing and averaging possible choices in Table 1 because it is not clear how to apply SN with weight sharing, and there is almost no performance gain with averaging possible choices despite of its high computational costs. Furthermore, we do not experiment on Movielens 10M with item vector encoding because the authors of CF-NADE did not provide the results for it due to the complexity of the model.
CF-UIcA (Du et al., 2018) We use the authors’ official code and the train/test dataset splits9 for these experiments. Since CF-UIcA is a model that accepts both user and item vector as input, it is not necessary to consider two types of encoding as in AutoRec or CF-NADE. On the other hand, it is reasonable to take different K values for user and item vector with SN. We treat K as 66 for the user vector and 110 for the item vector. For models without SN, 0.0001 is used as the parameter λ for weight decay as suggested by the CF-UIcA. It is natural to use different λ values for SN because the optimal λ must be changed along with SN. Hence, we use λ = 0.0005 for Movielens 100K and λ = 0.00006 for Movielens 1M in the case of SN. As in CF-NADE, we do not test Movielens 10M because the authors of CF-UIcA also did not report for it owing to its high computational cost.
8The CF-NADE uses learning rate 5 × 10−4 for Movielens 10M, but we use 10−3 for fast convergence. Therefore, the results can be somewhat different from original paper.
9https://github.com/thu-ml/CF-UIcA.

heading: B.2 DETAILED EXPERIMENTAL SETTINGS OF TABLE 2
text: In comparison with other state-of-the-art models, we used 1000 hidden units in AutoRec (Sedhain et al., 2015) with SN. While Sedhain et al. (2015) claimed that they were able to achieve enough performance only with 500 hidden units, 500 hidden units did not achieve sufficient performance when applying SN. The Figure 6 plots the test RMSE, changing the number of hidden units for Movielens 100K and 1M. We can see that 600 units for Movielens 100K and 900 units for Movielens 1M are necessary for getting better performance. Obviously, as datasets become more complex and larger, we need more network capacity. Therefore, we decide to use two times larger network capacity (1000 hidden units) to get better performance. The number of hidden units can also be viewed and tuned as a hyper-parameter, and we have not tuned much for the number of hidden units. Note that, unlike Table 1, we report the results with five random splits for all datasets in Table 2 to compare fairly with other state-of-the-art methods.

heading: C ELECTRONIC MEDICAL RECORDS (EMR) DATASETS

heading: C.1 NHIS DATASET
text: The NHIS dataset, which is from National Health Insurance Service (NHIS), consists of medical diagnosis of around 300,000 people. The goal is to predict the occurrence of 5 diseases. Each patients takes 34 examinations over 5 years. We split dataset into two set (train and test), where the ratio of train and test split is 3:1. We pre-process input data with min-max normalization, which makes min and max values of each features be zero and one following GAIN (Yoon et al., 2018). We train 2 layer neural networks which have 50 and 30 hidden units each, and evaluate the model with AUROC. We use ReLU activation, and dropout rate as 0.8 (if applied). Since these dataset is too imbalance, we apply class weight to the loss function for handling label imbalance. Besides, we use Adam optimizer with learning rate 10−2 without weight decay, and full batch. We evaluate the model on 5 times and report mean and 95% confidence interval. We also observe that Sparsity Normalization makes performance gain even when the dropout is integrated in the networks (See Table 5).
Data source This study used the National Health Insurance System-National Health Screening Cohort (NHIS-HEALS)* data derived from a national health screening program and the national health insurance claim database in the National Health Insurance System (NHIS) of South Korea. Data from the NHIS-HEALS10 was fully anonymized for all analyses and informed consent was not specifically obtained from each participant. This study was approved and exempt from informed consent by the Institutional Review Board of Yonsei University, Severance Hospital in Seoul, South Korea (IRB no.4-2016-0383).
Data Availability Data cannot be shared publicly because of the provisions of the National Health Insurance Service (NHIS). Korean legal restrictions prohibit authors from making the data publicly available, and the authority implemented the restrictions is NHIS (National Health Insurance Service), one of the government agency of Republic of Korea. NHIS provides limited portion of anonymized data to the researchers for the purpose of the public interest. However, they exclusively provide data to whom made direct contact of the NHIS and agreed to policies of NHIS. Redistribution of the data is not permitted for the researchers. The contact name and the information to which the data request can be sent: Haeryoung Park Information analysis department Big data operation room NHISS Tel: +82-33-736-2430. E-mail: lumen77@nhis.or.kr.
10Seong SC, Kim YY, Park SK, et al. Cohort profile: the National Health Insurance Service-National Health Screening Cohort (NHIS-HEALS) in Korea. BMJ Open 2017;7:e016640. pmid:28947447.

heading: C.2 PHYSIONET CHALLENGE 2012 DATASET
text: PhysioNet Challenge 2012 dataset (Silva et al., 2012) consists of 48 hours-long multivariate clinical time series data from intensive care unit (ICU). Most of the experimental settings are followed by BRITS (Cao et al., 2018). We divide 48 hours into 288 timesteps, which contain 35 examinations each. The goal of this task is to predict in-hospital death. We use dataset split given by PhysioNet Challenge 2012 (Silva et al., 2012) where each of them contains 4000 data points. In the preprocessing phase, we standardize (make mean as zero and standard deviation as one for each features) the input features and fill zero for missing values (zero imputation). We train single layer LSTM network which has 108 hidden units, and evaluate the model with AUROC. We apply class weight to handle imbalance problem in the dataset like in setting above. We use Adam optimizer with learning rate 2 × 10−4, 512 batch size, and early stopping method based on AUROC of validation set. We evaluate the model on 5 times and report mean and 95% confidence interval. Note that input of LSTM model is the results for 35 medical examinations at a specific timestamp. Hence, we apply SN separately for each timestamp. As the aforementioned state of Section 4.2, SN could not make significant performance gain in PhysioNet Challenge 2012 dataset (Silva et al., 2012) though Sparsity Normalization eases the VSP (See Table 6). Nevertheless, SN is still valuable for its ability to prevent biased predictions in this mission-critical area.

heading: D SINGLE-CELL RNA SEQUENCE DATASETS
text: In the AutoImpute experiments, we run the experiments using the author’s public code11.Talwar et al. (2018) reported experimental results on eight datasets (Blakeley, Jurkat, Kolodziejczyk, Preimplantation, Quake, Usoskin, Zeisel, and PBMC). Since we can obtain preprocessed datasets only for seven datasets except PBMC, we run experiments on seven single cell RNA sequence datasets12. All experimental settings without SN are exactly followed by the author’s code and the hyper-parameter settings published in Table 2 of their paper. For the model integrated with SN, all the experimental settings are followed by the original model except that the smaller threshold for early stopping is taken because the models with SN tends to be underfitted by using the threshold as suggested by the authors13. In addition, Talwar et al. (2018) conducted experiments only on cases with test set ratios of {0.1, 0.2, 0.3, 0.4, 0.5}, but we explored more test set ratios {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9} to show that the SN works well under more sparsity (extremely high test set ratio). Like the authors’ settings, we perform 10 experiments and report mean and 95% confidence intervals. We change the test set split on every trials. It is remarkable that although hyper-parameters are not favorable for SN, SN performs similarly or better on all seven datasets that we concerned (See Figure 3 and 7).

heading: E DROPOUT ON UCI DATASETS
text: Most experimental settings are adopted from Klambauer et al. (2017); Littwin & Wolf (2018)’s UCI experiments. We use ReLU (Glorot et al., 2011) networks of 4 hidden layers with 256 units. We use Mean Square Error (MSE) for loss function and Adam optimizer without weight decay, = 10−8 and learning rate 10−4. The batch size used for training is 128. We split 10% from the dataset for test and use 20% of training set as the validation set. For all inputs, we applied min-max normalization, which makes min and max values of each features be 0 and 1. Note that, we use dense datasets without any missing attributes to focus on the effects of dropout. All the datasets are used as provided in the package in sklearn.datasets.
11https://github.com/divyanshu-talwar/AutoImpute 12https://drive.google.com/drive/folders/1q2ho_cNfsQJNbdCt9j0nwlZv-Roj_
yK1 13We slightly tune the hyper-parameter λ of weight decay only for the Preimplantation dataset with SN (λ = 20).

heading: F DENSITY ESTIMATION
text: To reproduce binarized MNIST experiments of MADE (Germain et al., 2015), we adopt and slightly modify a public implementation of MADE14. Figure 5 shows our reproduction results of Figure 2 in the original paper. We follow most settings of the MADE. We use a single hidden layer MADE networks with 500 units, learning rate 0.005, and test on authors’ binarized MNIST dataset15. The only difference from the original authors’ implementation is that we used Adam ( = 10−8, and no weight decay), rather using Adagrad. Because it is observed that the model is undefitted with Adagrad while the learning rates of each element of weight matrix is dropped so rapidly.
Since there is no missingness in the Binarized MNIST dataset, it cannot be divided by ‖M‖1 as suggested by Algorithm 1. In this experiment, we regard ‖M‖1 as ‖h0‖0 so as not to lose the generality. That is, all pixels that are 0 in the binarized MNIST dataset are regarded as missing. We label results of this with w/SN and plot them in Figure 5. As the aforementioned state in Section 4.5, MADE can also cause variable sparsity by mask matrices of each weight. In MADE, the connection between specific units is forcibly controlled through a mechanism of element-wise product of a specific mask matrix M i in weight W i ∈ Rni×ni−1 . These mask matrices also cause variation in sparsity. We plot the results of using the new mask matrix M iSN in place of the mask matrix M
i used by MADE with w/SN(all) in Figure 5. The method of calculating the new mask matrix M iSN using the existing mask matrix M i is as follows:
M iSN ← (1TM i1/ni) ·M i (M i11T )
where denotes element-wise division and 1 is a column vector where all elements are 1. (1TM i1/ni) corresponds to K and (M i11T ) does to ‖m‖1 in Algorithm 1.
We demonstrate the effectiveness of SN even in situations where the learning rate is smaller. Smaller learning rate (0.001) shows the effect of SN more clearly despite of their longer training time as show in Figure 8.

heading: G MACHINE DESCRIPTION
text: We perform all the experiments on a Titan X with 12GB of VRAM. 12 GB of VRAM is not always necessary, and most experiments require smaller VRAM.
14https://github.com/karpathy/pytorch-made 15https://github.com/mgermain/MADE/releases/download/ICML2015/binarized_
mnist.npz

heading: H COMPARISON TO OTHER MISSING HANDLING TECHNIQUES
text: In this section, we compare Sparsity Normalization with other missing handling techniques. Although the main contribution of our paper is to provide a deeper understanding and the corresponding solution about the issue that the zero imputation, the simplest and most intuitive way of handling missing data, degrades the performance in training neural networks. Nonetheless, we show that Sparsity Normalization is effective for high dimensional datasets (with a large number of features and high missing rates) via a collaborative filtering dataset, while showing that Sparsity Normalization has competitive results with other modern missing handling techniques for datasets with non-high dimensional setting (electronic medical records datasets, UCI datasets w/ and w/o MCAR assumption). For fair comparisons, we consider the tasks in Section 4, as well as the tasks used by modern missing handling techniques such as GAIN (Yoon et al., 2018) and Śmieja et al. (2018).
As baselines, we consider modern missing handling techniques such as GAIN, Śmieja et al. (2018) and their baselines as well.
• Zero Imputation w/o Sparsity Normalization: Missing values are replaced with zero. • Zero Imputation w/ Sparsity Normalization (ours.): Based on zero imputation, apply Sparsity
Normalization (SN).
• Zero Imputation with Batch Normalization (Ioffe & Szegedy, 2015): Based on zero imputation, apply Batch Normalization (BN) only on the first layer.
• Zero Imputation with Layer Normalization (Lei Ba et al., 2016): Based on zero imputation, apply Layer Normalization (LN) only on the first layer. • Dropout16: Missing values are replaced with zero and other values are divided by( E(h0,m)∈D[‖m‖1]/n0 ) like standard dropout (Srivastava et al., 2014). Dropout uses a sin-
gle missing (drop) probability uniformly across all instances of the dataset while SN normalizes each data instance with its own missing rate.
• Mean Imputation: Missing values are replaced with the mean of those features. • Median Imputation: Missing values are replaced with the median of those features. • k-Nearest Neighbors (k-NN): Missing values are replaced with the mean of those features from k nearest neighbor samples. We use k = 5 following Śmieja et al. (2018)’s experimental setting.
• Multivariate Imputation by Chained Equations (MICE): Proposed by Buuren & GroothuisOudshoorn (2010).
• SoftImpute (Mazumder et al., 2010)17
• Gaussian Mixture Model Compensator (GMMC): Proposed by Śmieja et al. (2018). In the case of GMMC, any activation functions except ReLU and RBF (Radial Basis Function) is prohibitive on the first hidden layer. Thus, the activation function of the first hidden layer is replaced by ReLU in all base architectures without ReLU.
• GAIN (Yoon et al., 2018)
We implement Mean Imputation, Median Imputation, MICE, k-NN and SoftImpute by python package fancyimpute. We use authors’ official codes for GMMC18 and GAIN19. Layer Normalization and Batch Normalization are not commonly considered in studies of handling missing data. However, we additionally take these as baselines because they have similarities to Sparsity Normalization in terms of stabilizing the statistics of a hidden layer20 (See Appendix H.1.2 for deeper analysis).
16The GMMC (Śmieja et al., 2018) used this method as their baseline. 17In Yoon et al. (2018), this method was named Matrix. 18https://github.com/lstruski/Processing-of-missing-data-by-neural-networks 19https://github.com/jsyoon0823/GAIN 20We only consider LN and BN in the case of applying the first hidden layer. Because we find that the
prediction performance is more worse when LN or BN applied to all the hidden layers, and it is difficult to fairly compare with the Sparsity Normalization.

heading: H.1 COLLABORATIVE FILTERING (RECOMMENDATION) DATASET
text: In this section, we compare the SN with other missing handling techniques using the collaborative filtering dataset. Appendix H.1.1 compares the prediction performance of the baseline methods and the SN, while Appendix H.1.2 deeply analyzes the characteristics of the SN in comparison with Layer Normalization (Lei Ba et al., 2016) and Batch Normalization (Ioffe & Szegedy, 2015).

heading: H.1.1 COMPARISON OF PREDICTION PERFORMANCE
text: It is considered training an AutoRec (Sedhain et al., 2015) model on the Movielens 100K dataset. Most experimental settings are adopted from Section 4.121. We evaluate each missing handling techniques on both data encoding (user- or item-rating vector) as shown in Table 7. In both encoding, Sparsity Normalization performs better or similar to other missing handling techniques. While some missing handling techniques perform poorly rather than zero imputation depending on the encoding, Sparsity Normalization improves performance consistently for both data encodings. It is worth mentioning that Sparsity Normalization performs statistically significantly better than all other baselines with item vector encoding, which is considered a better encoding scheme in most collaborative filtering models (Salakhutdinov et al., 2007; Sedhain et al., 2015; Zheng et al., 2016).
H.1.2 IS BATCH NORMALIZATION OR LAYER NORMALIZATION ABLE TO SOLVE VSP?
Someone might wonder if Batch Normalization (Ioffe & Szegedy, 2015) or Layer Normalization (Lei Ba et al., 2016) have a similar effect to Sparsity Normalization by alleviating VSP. However, BN or LN can not solve VSP even though these three methods have something in common in terms of stabilizing the statistics of the hidden layer. To validate this, we compare SN with LN and BN by controlling the strength of weight decay regularization. We use the AutoRec model with Movielens 100K for these experiments as in Appendix H.1.1.
Figure 9 shows that VSP occurs in all cases except SN with weak regularization (left column). The model’s prediction highly correlates the number of known entries in all methods except SN. On the other hand, strong regularization might seem to solve the VSP while the model shows relatively constant inference regardless of the number of known entries. However, the strong regularization is not an acceptable solution because it gives the less freedom of the model’s inference making constant prediction (right column). It must not be natural that the predicted values of the model are almost constant regardless of input sample (we do not want a model that recommends the same movie no
21Only when applying Batch Normalization and Layer Normalization, we set the number of early stop iteration as 10, 000 (10 times of the other models) in order to prevent underfitting.
matter which movies a user like/dislike!). This trend can also be seen in the process of tuning the hyper-parameter λ for each model through that the optimal λ value of each model except LN and BN is determined at around 500, whereas that of BN and LN is determined above 500000 (inordinate regularization). In other words, unlike SN, the VSP is not solved with LN or BN. Rather, strong regularization is able to solve the VSP, but this is not a direct solution to VSP forcing the model to choice constant predicted values irrespectively of the input. It is instructive note that the trend of Figure 9 is extremely consistent with the test points as Figure 1.

heading: H.2 ELECTRONIC MEDICAL RECORDS (EMR) DATASETS
text: We also compare Sparsity Normalization and baselines for the five disease identification tasks in the NHIS dataset used in Section 4.2. The results are described in Table 3. Sparsity Normalization shows better or similar performance compared to other baseline methods as well.

heading: H.3 UCI DATASETS
text: We further compare Sparsity Normalization with other missing handling techniques on UCI datasets which have relatively low missing rates and small feature dimension (non-high dimensional datasets). We consider the UCI datasets used in GAIN (Yoon et al., 2018) and GMMC (Śmieja et al., 2018). The datasets used in the both papers can be divided into two categories: missing features are intentionally injected (w/ MCAR assumption) or missing features exist inherently (w/o MCAR assumption).
We adopt the settings of Klambauer et al. (2017); Littwin & Wolf (2018)’s UCI exeperiments to use the same Multi Layer Perceptron (MLP) architecture as in Section 4.4: ReLU (Glorot et al., 2011) networks of 4 hidden layers with 256 units. The main purpose of imputation should be to improve prediction performance rather than imputation performance. In this reason, we just focus on the prediction performance of each missing handling techniques for UCI datasets. Because all UCI datasets used in this section are for imbalanced binary classification tasks, prediction performance is reported with AUROC rather than accuracy, and the class weight is considered in loss function. On top of that, we use Adam Optimizer in all experiments for fair comparison with baselines.
Though we adopt datasets used in GAIN (Yoon et al., 2018) and GMMC (Śmieja et al., 2018), we report quite different performance from that of the papers. Several possible reasons are as follows. First, GAIN and GMMC did not publish the train/test dataset split, thus we use our own split which is made under the similar settings of both papers. Second, MLP is used rather than logistic regression or Radial Basis Function Network (RBFN) which are used in GAIN and GMMC respectively. It is because we think that MLP is more reasonable and widely acceptable architecture nowadays22 than the others. Furthermore, we use AUROC and class weights, unlike the GAIN and GMMC. The final possible reason is that we use Adam Optimizer for all models because SGD with learning rate decay is difficult for fair comparison when hyper-parameters are set in favor to a particular model. In these ways, we do our best to compare Sparsity Normalization and other missing handling techniques including GMMC and GAIN in the most fair and reasonable setting.

heading: H.3.1 UCI DATASETS WITH MCAR ASSUMPTION
text: We deliberately inject missing values into the datasets which don’t have any missing attributes internally (w/ MCAR assumption) to perform binary classification tasks. We consider Breast, Spam, and Credit datasets from GAIN (Yoon et al., 2018) and Crashes and Heart datasets from GMMC (Śmieja et al., 2018). We make 20% of all features be missing for the Breast, Spam, and Credit datasets following GAIN, and 50% for the Crashes and Heart datasets following GMMC. The summary of the datasets are described in Table 8. For Breast, Spam, and Credit datasets taken by GAIN, we perform min-max normalization which makes min and max values of each features be 0 and 1 following GAIN paper, and for Crashes and Heart datasets taken by GMMC, we perform another kind of min-max normalization which makes min and max values of each features be -1 and 1 following GMMC paper.
The experimental results are summarized in Table 9. It is difficult to find a significant difference in prediction performance among each missing handling techniques for datasets with small feature dimensions. The results of these experiments are also consistent with experiments of GAIN (Yoon et al., 2018). Though the GAIN showed significantly better imputation performance compared to their baseline methods, prediction performances were not statistically significant (See Table 3 of the GAIN paper, and note that they didn’t report 95% confidence interval but standard deviation). From these overall results, we conclude that SN is quite comparable for the datasets of low dimension/missing rate with MCAR assumption.

heading: H.3.2 UCI DATASETS WITHOUT MCAR ASSUMPTION
text: We compare SN and each missing handling techniques on the datasets that have internal missingness (w/o MCAR assumption). We consider Bands, Hepartitis, Horse, Mammographics, and Pima datasets experimented in the GMMC (Śmieja et al., 2018) paper (See Table 8 for statistics of the datasets). Following the GMMC paper, the min-max normalization is performed, which makes min and max values of each features be -1 and 1. As shown in Table 10, it is concluded that even without MCAR assumption, SN shows comparable results for the datasets of low dimension/missing rate.

heading: H.4 CONCLUSION
text: In conclusion, Sparsity Normalization is significantly superior to other missing handling techniques for the high dimensional/missing rate datasets. Sparsity Normalization performs well compared to other modern missing handling techniques even on non-high dimensional datasets. Sparsity Normalization is valuable in that it performs better than other models and does not require additional training or parameters. Moreover, Sparsity Normalization is computationally inexpensive compared to Mean or Median Imputation because sparse tensors can be used to save computational costs to calculate first hidden layer not with mean or median imputation but with zero imputation (w/ or w/o SN). The reduced computational cost by using sparse tensors is relatively higher when we deal with high-dimensional datasets.

Review summaries:
This paper provides a novel solution to the variable sparsity problem, where the output of neural networks biased with respect to the number of missing inputs. The authors proposed a sparsity normalization algorithm to process the input vectors to encounter the bias. In experiments, the authors evaluated the proposed sparsity normalization model on multiple datasets: collaborative filtering datasets, electric medical records datasets, single-cell RNA sequence datasets and UCI datasets. Results show that the proposed normalization method improves the prediction performance and the predicted values of the neural network is more uniformly distributed according to the number of missing entries.

This paper studies a very interesting phenomena in machine learning called VSP, that is the output of the model is highly affected via the level of missing values in its input.  The authors demonstrate the existence of such phenomena empirically, analyze the root cause for it theoretically, and propose a simple yet effective normalization method to tackle the problem. Several experiments demonstrate the effectiveness of this method.