Year: 2021

Abstract:
Irregularly sampled time series commonly occur in several domains where they present a significant challenge to standard deep learning models. In this paper, we propose a new deep learning framework for probabilistic interpolation of irregularly sampled time series that we call the Heteroscedastic Temporal Variational Autoencoder (HeTVAE). HeTVAE includes a novel input layer to encode information about input observation sparsity, a temporal VAE architecture to propagate uncertainty due to input sparsity, and a heteroscedastic output layer to enable variable uncertainty in output interpolations. Our results show that the proposed architecture is better able to reflect variable uncertainty through time due to sparse and irregular sampling than a range of baseline and traditional models, as well as recently proposed deep latent variable models that use homoscedastic output layers.

heading: 1 INTRODUCTION
text: In this paper, we propose a novel deep learning framework for probabilistic interpolation of irregularly sampled time series. Irregularly sampled time series data occur in multiple scientific and industrial domains including finance (Manimaran et al., 2006), climate science (Schulz & Stattegger, 1997) and healthcare (Marlin et al., 2012; Yadav et al., 2018). In some domains including electronic health records and mobile health studies (Cheng et al., 2017), there can be significant variation in inter-observation intervals through time. This is due to the complexity of the underlying observation processes that can include “normal” variation in observation times combined with extended, blockstructured periods of missingness. For example, in the case of ICU EHR data, this can occur due to patients being moved between different locations for procedures or tests, resulting in missing physiological sensor data for extended periods of time. In mobile health studies, the same problem can occur due to mobile sensor batteries running out, or participants forgetting to wear or carry devices.
In such situations, it is of critical importance for interpolation models to be able to correctly reflect the variable input uncertainty that results from variable observation sparsity so as not to provide overly confident inferences. However, modeling time series data subject to irregular sampling poses a significant challenge to machine learning models that assume fully-observed, fixed-size feature representations (Marlin et al., 2012; Yadav et al., 2018; Shukla & Marlin, 2021b). The main challenges in dealing with such data include the presence of variable time gaps between the observation time points, partially observed feature vectors caused by the lack of temporal alignment across different dimensions, as well as different data cases, and variable numbers of observations across dimensions and data cases. Significant recent work has focused on developing specialized models and architectures to address these challenges in modeling irregularly sampled multivariate time series (Li & Marlin, 2015; 2016; Lipton et al., 2016; Futoma et al., 2017; Che et al., 2018a; Shukla & Marlin, 2019; Rubanova et al., 2019; Horn et al., 2020; Li & Marlin, 2020; Shukla & Marlin, 2021a; De Brouwer et al., 2019; Tan et al., 2020; Kidger et al., 2020).
Recently, Shukla & Marlin (2021a) introduced the Multi-Time Attention Network (mTAN) model, a variational autoencoder (VAE) architecture for continuous-time interpolation of irregularly sampled time series. This model was shown to provide state-of-the-art classification and deterministic interpolation performance. However, like many VAEs, the mTAN architecture produces a homoscedastic output distribution conditioned on the latent state. This means that the model can only reflect uncertainty due to variable input sparsity through variations in the VAE latent state. As we will show,
this mechanism is insufficient to capture differences in uncertainty over time. On the other hand, Gaussian Process Regression-based (GPR) methods (Rasmussen & Williams, 2006) have the ability to reflect variable uncertainty through the posterior inference process. The main drawbacks of GPR-based methods are their significantly higher run times during both training and inference, and the added restriction to define positive definite covariance functions for multivariate time series.
In this work, we propose a novel encoder-decoder architecture for multivariate probabilistic time series interpolation that we refer to as the Heteroscedastic Temporal Variational Autoencoder or HeTVAE. HeTVAE aims to address the challenges described above by encoding information about input sparsity using an uncertainty-aware multi-time attention network (UnTAN), flexibly capturing relationships between dimensions and time points using both probabilistic and deterministic latent pathways, and directly representing variable output uncertainty via a heteroscedastic output layer.
The proposed UnTAN layer generalizes the previously introduced mTAN layer with an additional intensity network that can more directly encode information about input uncertainty due to variable sparsity. The proposed UnTAN layer uses an attention mechanism to produce a distributed latent representation of irregularly sampled time series at a set of reference time points. The UnTAN module thus provides an interface between input multivariate, sparse and irregularly sampled time series data and more traditional deep learning components that expect fixed-dimensional or regularly spaced inputs. We combat the presence of additional local optima that arises from the use of a heteroscedastic output layer by leveraging an augmented training objective where we combine the ELBO loss with an uncertainty agnostic loss component. The uncertainty agnostic component helps to prevent learning from converging to local optima where the structure in data is explained as noise.
We evaluate the proposed architecture on both synthetic and real data sets. Our approach outperforms a variety of baseline models and recent approaches in terms of log likelihood, which is our primary metric of interest in the case of probabilistic interpolation. Finally, we perform ablation testing of different components of the architecture to assess their impact on interpolation performance.

heading: 2 RELATED WORK
text: Keeping in mind the focus of this work, we concentrate our discussion of related work on deterministic and probabilistic approaches applicable to the interpolation and imputation tasks.
Deterministic Interpolation Methods: Deterministic interpolation methods can be divided into filtering and smoothing-based approaches. Filtering-based approaches infer the values at a given time by conditioning only on past observations. For example, Han-Gyu Kim et al. (2017) use a unidirectional RNN for missing data imputation that conditions only on data from the relative past of the missing observations. On the other hand, smoothing-based methods condition on all possible observations (past and future) to infer any unobserved value. For example, Yoon et al. (2018) and Cao et al. (2018) present missing data imputation approach based on multi-directional and bi-directional RNNs. These models typically use the gated recurrent unit with decay (GRU-D) model (Che et al., 2018a) as a base architecture for dealing with irregular sampling. Interpolation-prediction networks take a different approach to interfacing with irregularly sampled data that is based on the use of temporal kernel smoother-based layers (Shukla & Marlin, 2019). Shan & Oliva (2021) propose hierarchical imputation strategy based on set-based architectures for imputation in irregularly sampled time series. Of course, the major disadvantage of deterministic interpolation approaches is that they do not express uncertainty over output interpolations and thus can not be applied to the problem of probabilistic interpolation without modifications.
Probabilistic Interpolation Methods: The two primary building blocks for probabilistic interpolation and imputation of multivariate irregularly sampled time series are Gaussian processes regression (GPR) (Rasmussen & Williams, 2006) and variational autoencoders (VAEs) (Rezende et al., 2014; Kingma & Welling, 2014). GPR models have the advantage of providing an analytically tractable full joint posterior distribution over interpolation outputs when conditioned on irregularly sampled input data. Commonly used covariance functions have the ability to translate variable input observation density into variable interpolation uncertainty. GPR-based models have been used as the core of several approaches for supervised learning and forecasting with irregularly sampled data (Ghassemi et al., 2015; Li & Marlin, 2015; 2016; Futoma et al., 2017). However, GPR-based models can become somewhat cumbersome in the multivariate setting due to the positive definiteness con-
straint on the covariance function (Rasmussen & Williams, 2006). The use of separable covariance functions is one common approach to the construction of GPR models over multiple dimensions (Bonilla et al., 2008), but this construction requires all dimensions to share the same temporal kernel parameters. A further drawback of GP-based methods is their significantly higher run times relative to deep learning-based models when applied to larger-scale data (Shukla & Marlin, 2019).
Variational autoencoders (VAEs) combine probabilistic latent states with deterministic encoder and decoder networks to define a flexible and computationally efficient class of probabilistic models that generalize classical factor analysis (Kingma & Welling, 2014; Rezende et al., 2014). Recent research has seen the proposal of several new VAE-based models for irregularly sampled time series. Chen et al. (2018) proposed a latent ordinary differential equation (ODE) model for continuous-time data using an RNN encoder and a neural ODE decoder. Building on the prior work of Chen et al. (2018), Rubanova et al. (2019) proposed a latent ODE model that replaces the RNN with an ODE-RNN model as the encoder. Li et al. (2020) replace the deterministic ODEs with stochastic differential equations(SDEs). Norcliffe et al. (2021) extends the prior work on neural ode by combining it with neural processes (Garnelo et al., 2018). Shukla & Marlin (2021a) proposed the Multi-Time Attention Network (mTAN) model, a VAE-based architecture that uses a multi-head temporal cross attention encoder and decoder module (the mTAND module) to provide the interface to multivariate irregularly sampled time series data. Fortuin et al. (2020) proposed a VAE-based approach for the task of smoothing in multivariate time series with a Gaussian process prior in the latent space to capture temporal dynamics. Garnelo et al. (2018); Kim et al. (2019) used heteroscedastic output layers to represent uncertainty in case of fixed dimensional inputs but these approaches are not applicable to irregularly sampled time series.
Similar to the mTAN model, the Heteroscedastic Temporal Variational Autoencoder (HeTVAE) model proposed in this work is an attention-based VAE architecture. The primary differences are that mTAN uses a homoscedastic output distribution that assumes constant uncertainty and that the mTAN model’s cross attention operation normalizes away information about input sparsity. These limitations are problematic in cases where there is variable input density through time resulting in the need for encoding, propagating, and reflecting that uncertainty in the output distribution. As we describe in the next section, HeTVAE addresses these issues by combining a novel sparsity-sensitive encoder module with a heteroscedastic output distribution and parallel probabilistic and deterministic pathways for propagating information through the model. Another important difference relative to these previous methods is that HeTVAE uses an augmented learning objective to address the underfitting of predictive variance caused by the use of the heteroscedastic layer.

heading: 3 PROBABILISTIC INTERPOLATION WITH THE HETVAE
text: In this section, we present the proposed architecture for probabilistic interpolation of irregularly sampled time series, the Heteroscedastic Temporal Variational Autoencoder (HeTVAE). HeTVAE leverages a sparsity-aware layer as the encoder and decoder in order to represent input uncertainty and propagate it to output interpolations. We begin by introducing notation. We then describe the architecture of the encoder/decoder network followed by the complete HeTVAE architecture.

heading: 3.1 NOTATION
text: We letD = {sn|n = 1, ..., N} represent a data set containing N data cases. An individual data case consists of a D-dimensional, sparse and irregularly sampled multivariate time series sn. Different dimensions d of the multivariate time series can have observations at different times, as well as different total numbers of observations Ldn. We follow the series-based representation of irregularly sampled time series (Shukla & Marlin, 2021b) and represent time series d for data case n as a tuple sdn = (tdn,xdn) where tdn = [t1dn, ..., tLdndn] is the list of time points at which observations are defined and xdn = [x1dn, ..., xLdndn] is the corresponding list of observed values. We drop the data case index n for brevity when the context is clear.

heading: 3.2 REPRESENTING INPUT SPARSITY
text: As noted in the previous section, the mTAN encoder module does not represent information about input sparsity due to the normalization of the attention weights. To address this issue, we propose
an augmented module that we refer to as an Uncertainty Aware Multi-Time Attention Network (UnTAN). The UnTAN module is shown in Figure 1a. This module includes two encoding pathways that leverage a shared time embedding function and a shared attention function. The first encoding pathway (the intensity pathway, INT) focuses on representing information about the sparsity of observations while the second encoding pathway (the value pathway, VAL) focuses on representing information about values of observations. The outputs of these two pathways are concatenated and mixed via a linear layer to define the final output of the module. The mathematical description of the module is given in Equations 1 to 3 and is explained in detail below.
inth(rk, td) = pool({exp(αh(rk, tid)) | tid ∈ td}) pool({exp(αh(rk, ti′u)) | ti′u ∈ tu})
(1)
valh(rk, td,xd) = pool({exp(αh(rk, tid)) · xid | tid ∈ td, xid ∈ xd})
pool({exp(αh(rk, ti′d)) | ti′d ∈ td}) (2)
αh(t, t ′) =
( φh(t)wv
Tφh(t ′)T√
de
) (3)
Time Embeddings and Attention Weights: Similar to the mTAN module, the UnTAN module uses time embedding functions φh(t) to project univariate time values into a higher dimensional space. Each time embedding function is a one-layer fully connected network with a sine function non-linearity φh(t) = sin(ω · t+ β). We learn H time embeddings each of dimension de. w and v are the parameters of the scaled dot product attention function αh(t, t′) shown in Equation 3. The scaling factor 1/ √ de is used to normalize the dot product to counteract the growth in the dot product magnitude with increase in the time embedding dimension de.
Intensity Encoding: The intensity encoding pathway is defined by the function inth(rk, td) shown in Equation 1. The inputs to the intensity function are a query time point rk and a vector td containing all the time points at which observations are available for dimension d. The numerator of the intensity function exponentiates the attention weights between rk and each time point in td to ensure positivity, then pools over the observed time points. The denominator of this computation is identical to the numerator, but the set of time points tu that is pooled over is the union over all observed time points for dimension d from all data cases.
Intuitively, if the largest attention weight between rk and any element of td is small relative to attention weights between rk and the time points in tu, then the output of the intensity function will be low. Importantly, due to the use of the non-linear time embedding function, pairs of time points with high attention weights do not necessarily have to be close together in time meaning the notion of intensity that the network expresses is significantly generalized.
We also note that different sets could be used for tu including a regularly spaced set of reference time points. One advantage of using the union of all observed time points is that it fixes the maximum value of the intensity function at 1. The two pooling functions applicable in the computation of the intensity function are max and sum. If the time series is sparse, max works well because using sum in the sparse case can lead to very low output values. In a more densely observed time series, either sum or max can be used.
Value Encoding: The value encoding function valh(rk, td,xd) is presented in Equation 2 in a form that highlights the symmetry with the intensity encoding function. The primary differences are that valh(rk, td,xd) takes as input both observed time points td and their corresponding values xd, and the denominator of the function pools over td itself. While different pooling options could be used for this function, in practice we use sum-based pooling. These choices lead to a function valh(rk, td,xd) that interpolates the observed values at the query time points using softmax weights derived from the attention function. The values of observed points with higher attention weights contribute more to the output value. This structure is equivalent to that used in the mTAN module when sum-based pooling is used. We can also clearly see that this function on its own can not represent information about input sparsity due to the normalization over td. Indeed, the function is completely invariant to an additive decrease in all of the attention weights α′h(rk, tid) = αh(rk, tid)− δ.
Module Output: The last stage of the UnTAN module concatenates the value and intensity pathway representations and then linearly weights them together to form the final J-dimensional repre-
sentation that is output by the module. The parameters of this linear stage of the model are U inthdj and Uvalhdj . The value of the j th dimension of the output at a query time point rk is given by Equation 4.
UnTAN(rk, t,x)[j] = H∑ h=1 D∑ d=1
[ inth(rk, td)
valh(rk, td,xd)
]T [ U inthdj
Uvalhdj
] (4)
Finally, we note that the UnTAN module defines a continuous function of t given an input time series and hence cannot be directly incorporated into standard neural network architectures. We adapt the UnTAN module to produce fully observed fixed-dimensional discrete sequences by materializing its output at a set of reference time points. Reference time points can be fixed set of regularly spaced time points or may need to depend on the input time series. For a given set of reference time points r = [r1, · · · , rK ], the discretized UnTAN module UnTAND(r, t,x) is defined as UnTAND(r, t,x)[i] = UnTAN(ri, t,x). This module takes as input the time series s = (t,x) and the set of reference time points r and outputs a sequence of K UnTAN embeddings, each of dimension J corresponding to each reference point. As described in the next section, we use the UnTAND module to provide an interface between sparse and irregularly sampled data and fully connected MLP network structures.

heading: 3.3 THE HETVAE MODEL
text: In this section, we describe the overall architecture of the HeTVAE model, as shown in Figure 1b.
Model Architecture: The HeTVAE consists of parallel deterministic and probabilistic pathways for propagating input information to the output distribution, including information about input sparsity. We begin by mapping the input time series s = (t,x) through the UnTAND module along with a collection of K reference time points r. In the probabilistic path, we construct a distribution over latent variables at each reference time point using a diagonal Gaussian distribution q with mean and variance output by fully connected layers applied to the UnTAND output embeddings henc = [henc1 , · · · ,hencK ] as shown in Equation 6. In the deterministic path, the UnTAND output embeddings henc are passed through a feed-forward network g to produce a deterministic temporal representation (at each reference point) of the same dimension as the probabilistic latent state.
The decoder takes as input the representation from both pathways along with the reference time points and a set of query points t′ (Eq 8). The UnTAND module produces a sequence of embeddings hdec = [hdec1 , · · · ,hdec|t′| ] corresponding to each time point in t
′. The UnTAND embeddings are then independently decoded using a fully connected decoder fdec and the result is used to parameterize the output distribution. We use a diagonal covariance Gaussian distribution where both the mean µ = [µ1, · · · ,µ|t′|],µi ∈ RD and variance σ2 = [σ21 , · · · ,σ2|t′|],σ 2 i ∈ RD are predicted for each time point by the final decoded representation as shown in Eq 9. The generated time series is sampled from this distribution and is given by ŝ = (t′,x′) with all data dimensions observed.
The complete model is described below. We define qγ(z|r, s) to be the distribution over the probabilistic latent variables z = [z1, · · · , zK ] induced by the input time series s = (t,x) at the reference time points r. We define the prior p(zi) over the latent states to be a standard multivariate normal distribution. We let phetθ (x ′ id | zcat, t′id) define the final probability distribution over the value of time point t′id on dimension d given the concatenated latent state z cat = [zcat1 , · · · , zcatK ]. γ and θ represent the parameters of all components of the encoder and decoder respectively.
henc = UnTANDenc(r, t,x) (5)
zk ∼ qγ(zk |µk,σ2k), µk = fencµ (henck ), σ2k = fencσ (henck ) (6) zcatk = concat(zk, g(h enc k )) (7)
hdec = UnTANDdec(t′, r, zcat) (8)
phetθ (x ′ id | zcat, t′id) = N (x′id; µi [d], σ2i [d]), µi = fdecµ (hdeci ), σ2i = fdecσ (hdeci ) (9)
x′id ∼ phetθ (x′id | zcat, t′id) (10)
Compared to the constant output variance used to train the mTAN-based VAE model proposed in prior work (Shukla & Marlin, 2021a), our proposed model produces a heteroscedastic output distribution that we will show provides improved modeling for the probabilistic interpolation task. However, the increased complexity of the model’s output representation results in an increased space of local optima. We address this issue using an augmented learning objective, as described in the next section. Finally, we note that we can easily obtain a simplified homoscedastic version of the model with constant output variance σ2c using the alternate final output distribution pcθ(x ′ id | z, t′id) = N (x′id; µi [d], σ2c ).
Augmented Learning Objective: To learn the parameters of the HeTVAE framework given a data set of sparse and irregularly sampled time series, we propose an augmented learning objective based on a normalized version of the evidence lower bound (ELBO) combined with an uncertainty agnostic scaled squared loss. We normalize the contribution from each data case by the total number of observations so that the effective weight of each data case in the objective function is independent of the total number of observed values. The augmented learning objective is defined below. µn is the predicted mean over the test time points as defined in Equation 9. Also recall that the concatenated latent state zcat depends directly on the probabilistic latent state z.
LNVAE(θ, γ) = N∑ n=1 1∑ d Ldn ( Eqγ(z|r,sn)[log p het θ (xn|zcatn , tn)]−DKL(qγ(z|r, sn)||p(z)) (11)
− λEqγ(z|r,sn)‖xn − µn‖ 2 2] )
DKL(qγ(z|r, sn)||p(z)) = K∑ i=1 DKL(qγ(zi|r, sn)||p(zi)) (12)
log phetθ (xn|zcatn , tn) = D∑ d=1 Ldn∑ j=1 log phetθ (xjdn|zcatn , tjdn) (13)
We include the uncertainty agnostic scaled squared loss term to counteract the propensity of the heteroscedastic model to become stuck in poor local optima where the mean is essentially flat and all of the structure in the data is explained as noise. This happens because the model has the ability to learn larger variances at the output, which allows the mean to underfit the data. The extra component (scaled squared loss) helps to push the optimization process to find more informative
Observed data Ground truth Reconstructions
n = 3 n = 10 n = 20
Figure 2: In this figure, we show example interpolations on the synthetic dataset. The set of 3 columns correspond to interpolation results with increasing numbers of observed points: 3, 10 and 20 respectively. The first, second and third rows correspond to STGP, HeTVAE and HTVAE mTAN respectively. The shaded region corresponds to ± one standard deviation. STGP and HetVAE exhibit variable output uncertainty in response to input sparsity while mTAN does not.
parameters by introducing a fixed penalty for the mean deviating from the data. As we will show in the experiments, the use of this augmented training procedure has a strong positive impact on final model performance. Since, we are focusing on the interpolation task, we train the HeTVAE by maximizing the augmented learning objective (Equation 11) on the interpolated time points (more details on training has been provided in the experimental protocols in Section 4).

heading: 4 EXPERIMENTS
text: In this section, we present interpolation experiments using a range of models on three real-world data sets. PhysioNet Challenge 2012 (Silva et al., 2012) and MIMIC-III (Johnson et al., 2016) consist of multivariate, sparse and irregularly sampled time series data. We also perform experiments on the Climate dataset (Menne et al., 2016), consisting of multi-rate time series. We also show qualitative results on a synthetic dataset. Details of each dataset can be found in the Appendix A.6.1.
Experimental Protocols: We randomly divide the real data sets into a training set containing 80% of the instances, and a test set containing the remaining 20% of instances. We use 20% of the training data for validation. In the interpolation task, we condition on a subset of available points and produce distributions over the rest of the time points. On the real-world datasets, we perform interpolation experiments by conditioning on 50% of the available points. At test time, the values of observed points are conditioned on and each model is used to infer single time point marginal distributions over values at the rest of the available time points in the test instance. In the case of methods that do not produce probabilistic outputs, we make mean predictions. In the case of the synthetic dataset where we have access to all true values, we use the observed points to infer the values at the rest of the available points. We repeat each real data experiment five times using different random seeds to initialize the model parameters. We assess performance using the negative log likelihood, which is our primary metric of interest. We also report mean squared and mean absolute error. For all experiments, we select hyper-parameters on the held-out validation set using grid search and then apply the best trained model to the test set. The hyper-parameter ranges searched for each model and dataset are fully described in Appendix A.5.
Models: We compare our proposed model HeTVAE to several probabilistic and deterministic interpolation methods. We compare to two Gaussian processes regression (GPR) approaches. The most basic GP model for multivariate time series fits one GPR model per dimension. This approach is known as a single task GP model (STGP) (Rasmussen & Williams, 2006). A potentially better option is to model data using a Multi Task GP (MTGP) (Bonilla et al., 2008). This approach models the correlations both across different dimensions and across time by defining a kernel expressed as the Hadamard product of a temporal kernel (as used in the STGP) and a task kernel. We also compare to several VAE-based approaches. These approaches use a homoscedastic output distribution with different encoder and decoder architectures. HVAE RNN employs a gated recurrent unit network (Chung et al., 2014) as encoder and decoder, HVAE RNN-ODE (Chen et al., 2018) replaces the RNN decoder with a neural ODE, HVAE ODE-RNN-ODE (Rubanova et al., 2019) employs a ODE-RNN encoder and neural ODE decoder. Finally, we compare to HTVAE mTAN (Shukla & Marlin, 2021a), a temporal VAE model consisting of multi-time attention networks producing
homoscedastic output. For VAE models with homoscedastic output, we treat the output variance term as a hyperparameter and select the variance using log likelihood on the validation set. Architecture details for these methods can be found in Appendix A.4. As baselines, we also consider deterministic mean and forward imputation-based methods. Forward imputation always predicts the last observed value on each dimension, while mean imputation predicts the mean of all the observations for each dimension.
Synthetic Data Results: Figure 2 shows sample visualization output for the synthetic dataset. For this experiment, we compare HTVAE mTAN, the single task Gaussian process STGP, and the proposed HeTVAE model. We vary the number of observed points (3, 10, 20) and each model is used to infer the distribution over the remaining time points. We draw multiple samples from the VAE latent state for HeTVAE and HTVAE mTAN and visualize the distribution of the resulting mixture. As we can see, the interpolations produced by HTVAE mTAN have approximately constant uncertainty across time and this uncertainty level does not change even when the number of points conditioned on increases. On the other hand, both HeTVAE and STGP show variable uncertainty across time. Their uncertainty reduces in the vicinity of input observations and increases in gaps between observations. Even though the STGP has an advantage in this experiment (the synthetic data were generated with an RBF kernel smoother and STGP uses RBF kernel as the covariance function), the proposed model HeTVAE shows comparable interpolation performance. We show more qualitative results in Appendix A.3.
Real Data Results: Tables 1, 2 and 3 compare the interpolation performance of all the approaches on PhysioNet, MIMIC-III and Climate dataset respectively. HeTVAE outperforms the prior approaches with respect to the negative log likelihood score on all three datasets. Gaussian Process based methods − STGP and MTGP achieve second and third best performance respectively. We emphasize that while the MAE and MSE values for some of the prior approaches are close to those obtained by the HeTVAE model, the primary metric of interest for comparing probabilistic interpolation approaches is log likelihood, where the HeTVAE performs much better than the other methods.
We note that the MAE/MSE of the VAE-based models with homoscedastic output can be improved by using a small fixed variance during training. However, this produces even worse log likelihood values. Further, we note that the current implementation of MTGP is not scalable to the Climate dataset (270 dimensions). We provide experiments on an additional dataset in Appendix A.1.
Ablation Results: Table 4 shows the results of ablating several different components of the HeTVAE model and training procedure. The first row shows the results for the full proposed approach. The HeTVAE - ALO ablation shows the result of removing the augmented learning objective and training the model only using the ELBO. This results in an immediate drop in performance on PhysioNet. HeTVAE - DET removes the deterministic pathway from the model, resulting in a performance drop on both MIMIC-III and PhysioNet. HeTVAE - INT removes the intensity encoding pathway from the UnTAND module. It results in a large drop in performance on both datasets. HeTVAE - HET- ALO removes the heteroscedastic layer and the augmented learning objective (since the augmented learning objective is introduced to improve the learning in the presence of heteroscedastic layer), resulting in a highly significant drop on both datasets. These results show that all of the components included in the proposed model contribute to improved model performance. We provide more ablation results in Appendix A.2 and discuss hyperparameter selection in Appendix A.5.

heading: 5 DISCUSSION AND CONCLUSIONS
text: In this paper, we have proposed the Heteroscedastic Temporal Variational Autoencoder (HeTVAE) for probabilistic interpolation of irregularly sampled time series data. HeTVAE consists of an input sparsity-aware encoder, parallel deterministic and probabilistic pathways for propagating input uncertainty to the output, and a heteroscedastic output distribution to represent variable uncertainty in the output interpolations. Furthermore, we propose an augmented training objective to combat the presence of additional local optima that arise from the use of the heteroscedastic output structure. Our results show that the proposed model significantly improves uncertainty quantification in the output interpolations as evidenced by significantly improved log likelihood scores compared to several baselines and state-of-the-art methods. While the HeTVAE model can produce a probability distribution over an arbitrary collection of output time points, it is currently restricted to producing marginal distributions. As a result, sampling from the model does not necessarily produce smooth trajectories as would be the case with GPR-based models. Augmenting the HeTVAE model to account for residual correlations in the output layer is an interesting direction for future work.

heading: 6 REPRODUCIBILITY STATEMENT
text: We have submitted the source code as part of the supplementary materials which contains the instructions to reproduce the results in the paper including the hyperparameters. The hyperparameter ranges searched for each model are fully described in Appendix A.5. The source code also includes the synthetic dataset generation process as well as one of the real-world dataset. The other datasets can be downloaded and prepared following the preprocessing steps notes in Appendix A.6.1.

heading: A APPENDIX

heading: A.1 ADDITIONAL RESULTS
text: We also perform experiments on the UCI electricity dataset (described in Appendix A.6.1). We follow the same experiment protocols described in Section 4. As we can see from Table 5, the proposed model HeTVAE outperforms the prior approaches across all three metrics.

heading: A.2 ABLATION STUDY
text: Tables 6 and 7 show the complete results of ablating several different components of the HeTVAE model and training procedure with respect to all three evaluation metrics on PhysioNet and MIMICIII respectively. We denote different components of the HeTVAE model as − HET: heteroscedastic output layer, ALO: augmented learning objective, INT: intensity encoding, DET: deterministic pathway. The results show selected individual and compound ablations of these components and indicate that all of these components contribute significantly to the model’s performance in terms of the negative log likelihood score. We provide detailed comments below.
Effect of Heteroscedastic Layer: Since the augmented learning objective is introduced to improve the learning in the presence of heteroscedastic layer, we remove the augmented learning objective (ALO) with the heteroscedastic layer (HET). This ablation corresponds to HeTVAE - HET - ALO. As we can see from both Table 6 and 7, this results in a highly significant drop in the log likelihood performance as compared to the full HeTVAE model on both datasets. However, it results in only a slight drop in performance with respect to MAE and MSE, which is sensible as the HET component only affects uncertainty sensitive performance metrics.
Effect of Intensity Encoding: HeTVAE - INT removes the intensity encoding pathway from the UnTAND module. It results in an immediate drop in performance on both datasets. We also compare the effect of intensity encoding after removing the deterministic pathway and the augmented learning objective. These ablations are shown in HeTVAE - DET - ALO and HeTVAE - INT - DET - ALO. The performance drop is less severe in this case because of the propensity of the heteroscedastic output layer to get stuck in poor local optima in the absence of the augmented learning objective (ALO).
Effect of Augmented Learning Objective: The HeTVAE - ALO ablation shows the result of removing the augmented learning objective and training the model only using only the ELBO. This results in an immediate drop in performance on PhysioNet. The performance drop is less severe on MIMIC-III. We further perform this ablation without the DET component and observe severe drops in performance across all metrics on both datasets. These ablations correspond to HeTVAE - DET and HeTVAE - DET - ALO. This shows that along with ALO component, the DET component also constrains the model from getting stuck in local optima where all of the structure in the data is explained as noise. We show interpolations corresponding to these ablations in Appendix A.3.1.
Effect of Deterministic Pathway: HeTVAE - DET removes the deterministic pathway from the model, resulting in a performance drop on both MIMIC-III and PhysioNet across all metrics. We further compare the performance of both the probabilistic and deterministic pathways in isolation as shown by ablation HeTVAE - DET - ALO and HeTVAE - PROB - ALO. We observe that the deterministic pathway HeTVAE - PROB - ALO outperforms the probabilistic pathway HeTVAE -
DET - ALO in terms of log likelihood on MIMIC-III while the opposite is true in case of PhysioNet. However, on both datasets using only the deterministic pathway (HeTVAE - PROB - ALO) achieves better MAE and MSE scores as compared to using only the probabilistic pathway (HeTVAE - DET - ALO).
A.3 VISUALIZATIONS
A.3.1 INTERPOLATIONS ON PHYSIONET
Figure 3 shows example interpolations on the PhysioNet dataset. Following the experimental setting mentioned in Section 4, the models were trained using all dimensions and the inference uses all dimensions. We only show interpolations corresponding to Heart Rate as an illustration. As we can see, the STGP and HeTVAE models exhibit good fit and variable uncertainty on the edges where there are no observations. We can also see that mTAN trained with homoscedastic output is not able to produce as good a fit because of the fixed variance at the output (discussed in Section 4).
The most interesting observation is the performance of HeTVAE - DET - ALO, an ablation of HeTVAE model that retains heteroscedastic output, but removes the deterministic pathways and the augmented learning objective. This ablation significantly underfits the data and performs similar to mTAN. This is an example of local optima that arises from the use of a heteroscedastic output layer where the mean is excessively smooth and all of the structure in the data is explained as noise. We address this with the use of augmented learning objective described in Section 3.3. As seen in the Figure 3, adding the augmented learning objective (HeTVAE - DET) clearly improves performance.

heading: A.3.2 SYNTHETIC DATA VISUALIZATIONS: SPARSITY
text: In this section, we show supplemental interpolation results on the synthetic dataset. The setting here is same as in Section 4. Figure 4 compares HTVAE mTAN, the single task Gaussian process STGP, the proposed HeTVAE model and an ablation of proposed model without intensity encoding HeTVAE - INT. We vary the number of observed points (3, 10, 20) and each model is used to infer the distribution over the remaining time points. We draw multiple samples from the VAE latent state for HeTVAE, HeTVAE - INT and HTVAE mTAN, and visualize the distribution of the resulting mixture. Figure 4 illustrates the interpolation performance of each of the models. As we can see, the interpolations produced by HTVAE mTAN have approximately constant uncertainty across time and this uncertainty level does not change even when the number of points conditioned on increases. On the other hand, both HeTVAE and STGP show variable uncertainty across time. Their uncertainty reduces in the vicinity of input observations and increases in gaps between observations. The HeTVAE-INT model performs slightly better than HTVAE mTAN model but it does not show variable uncertainty due to input sparsity like HeTVAE.

heading: A.3.3 SYNTHETIC DATA VISUALIZATIONS: INTER-OBSERVATION GAP
text: To demonstrate the effectiveness of intensity encoder (INT), we perform another experiment on synthetic dataset where we increase the maximum inter-observation gap between the observations. We follow the same training protocol as described in Section 4. At test time, we condition on 10 observed points with increasing maximum inter-observation gap. We vary the maximum interobservation gap from 20% to 80% of the length of the original time series. Each model is used to infer single time point marginal distributions over values at the rest of the available time points in the test instance.
Figure 5 shows the interpolations with increasing maximum inter-observation gap. STGP and HeTVAE show variable uncertainty with time and the uncertainty increases with increasing maximum inter-observation gap. On the other hand, HTVAE mTAN with homoscedastic output shows approximately constant uncertainty with time and also across different maximum inter-observation gaps. These results clearly show that HTVAE mTAN produces over-confident probabilistic interpolations over large gaps.
Furthermore, we show an ablation of the proposed model HeTVAE - INT, where we remove the intensity encoder and perform the interpolations. As we see from the figure, this leads to approximately constant uncertainty across time as well as different maximum inter-observation gaps. This shows that the HeTVAE model is not able to capture uncertainty due to input sparsity as effectively without the intensity encoder.

heading: A.4 ARCHITECTURE DETAILS
text: HeTVAE: Learnable parameters in the UnTAND architecture shown in Figure 1a include the weights of the three linear layers and the parameters of the shared time embedding functions. Each time embedding function is a one layer fully connected network with a sine function non-linearity. The two linear layers on top of embedding function are linear projections from time embedding dimension de to de/H where H is the number of time embeddings. Note that these linear layers do not share parameters. The third linear layer performs a linear projection from 2 ∗ D ∗ H to J . It takes as input the concatenation of the VAL encoder output and INT encoder output and produces an output of dimension J . de, H and J are all hyperparameters of the architecture. The ranges considered are described in the next section.
The HeTVAE model shown in the Figure 1b consists of three MLP blocks apart from the UnTAND modules. The MLP in the deterministic path is a one layer fully connected layer that projects the UnTAND output to match the dimension of the latent state. The remaining MLP blocks are twolayer fully connected networks with matching width and ReLU activations. The MLP in the decoder takes the output of UnTAND module and outputs the mean and variance of dimension D and sequence length t′. We use a softplus transformation on the decoder output to get the variance σi = 0.01 + softplus(fdecσ (h dec i )). Similarly, in the probabilistic path, we apply an exponential transformation to get the variance of the q distribution σ2k = exp(f enc σ (h enc k )). We use K reference time points regularly spaced between 0 and 1. K is considered to be a hyperparameter of the architecture. The ranges considered are described in the next section.
Baselines: For the HTVAE mTAN, we use a similar architecture as HeTVAE where we remove the deterministic path, heteroscedastic output layer and use the mTAND module instead of the UnTAND module (Shukla & Marlin, 2021a). We use the same architectures for the ODE and RNN-based VAEs as Rubanova et al. (2019).

heading: A.5 HYPERPARAMETERS
text: HeTVAE: We fix the time embedding dimension to de = 128. The number of embeddings H is searched over the range {1, 2, 4}. We search the number of reference points K over the range {4, 8, 16, 32}, the latent dimension over the range {8, 16, 32, 64, 128}, the output dimension of UnTAND J over the range {16, 32, 64, 128}, and the width of the two-layer fully connected layers over {128, 256, 512}. In augmented learning objective, we search for λ over the range {1.0, 5.0, 10.0}. We use the Adam Optimizer for training the models. Experiments are run for 2, 000 iterations with a learning rate of 0.0001 and a batch size of 128. The best hyperparameters are reported in the code. We use 100 samples from the probabilistic latent state to compute the evaluation metrics.
Ablations: We note that the ablations were not performed with a fixed architecture. For all the ablation models, we tuned the hyperparameters and reported the results with the best hyperparameter setting. We also made sure that the hyperparameter ranges for ablated models with just deterministic/probabilistic path were wide enough that the optimal ablated models did not saturate the end of the ranges for architectural hyper-parameter values including the dimensionality of the latent representations.
VAE Baselines: For VAE models with homoscedastic output, we treat the output variance term as a hyperparameter and select the variance over the range {0.01, 0.1, 0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0}. For HTVAE mTAN, we search the corresponding hyperparameters over the same range as HeTVAE. For ODE and RNN based VAEs, we search for GRU hidden units, latent dimension, the number of hidden units in the fully connected network for the ODE function in the encoder and decoder over the range {20, 32, 64, 128, 256}. For ODEs, we also search the number of layers in fully connected network in the range {1, 2, 3}. We use a batch size of 50 and a learning rate of 0.001. We use 100 samples from the latent state to compute the evaluation metrics.
Gaussian Processes: For single task GP, we use a squared exponential kernel. In case of multitask GP, we experimented with the Matern kernel with different smoothness parameters, and the
squared exponential kernel. We found that Matern kernel performs better. We use maximum marginal likelihood to train the GP hyperparameters. We search for learning rate over the range {0.1, 0.01, 0.001} and run for 100 iterations. We search for smoothness parameter over the range {0.5, 1.5, 2.5}. We search for the batch size over the range {32, 64, 128, 256}.

heading: A.6 TRAINING DETAILS

heading: A.6.1 DATA GENERATION AND PREPROCESSING
text: Synthetic Data Generation: We generate a synthetic dataset consisting of 2000 trajectories each consisting of 50 time points with values between 0 and 1. We fix 10 reference time points and draw values for each from a standard normal distribution. We then use an RBF kernel smoother with a fixed bandwidth of α = 120.0 to construct local interpolations over the 50 time points. The data generating process is shown below:
zk ∼ N (0, 1), k ∈ [1, · · · , 10] rk = 0.1 ∗ k ti = 0.02 ∗ i, i ∈ [1, · · · , 50]
xi = ∑ k exp(−α(ti − rk)2) · zk∑ k′ exp(−α(ti − rk′)2) +N (0, 0.12)
We randomly sample 3 − 10 observations from each trajectory to simulate a sparse and irregularly sampled univariate time series.
PhysioNet: The PhysioNet Challenge 2012 dataset (Silva et al., 2012) consists of multivariate time series data with 37 physiological variables from intensive care unit (ICU) records. Each record contains measurements from the first 48 hours after admission. We use the protocols described in Rubanova et al. (2019) and round the observation times to the nearest minute resulting in 2880 possible measurement times per time series. The data set consists includes 8000 instances that can be used for interpolation experiments. PhysioNet is freely available for research use and can be downloaded from https://physionet.org/content/challenge-2012/.
MIMIC-III: The MIMIC-III data set (Johnson et al., 2016) is a multivariate time series dataset containing sparse and irregularly sampled physiological signals collected at Beth Israel Deaconess Medical Center. We use the procedures proposed by Shukla & Marlin (2019) to process the data set. This results in 53, 211 records each containing 12 physiological variables. We use all 53, 211 instances to perform interpolation experiments. MIMIC-III is available through a permissive data use agreement which can be requested at https://mimic.mit.edu/iii/gettingstarted/. Once the request is approved, the dataset can be downloaded from https://mimic.mit.edu/ iii/gettingstarted/dbsetup/. The instructions and code to extract the MIMIC-III dataset is given at https://github.com/mlds-lab/interp-net.
Climate Dataset: The U.S. Historical Climatology Network Monthly (USHCN) dataset (Menne et al., 2016) is a publicly available dataset consisting of daily measurements of 5 climate variables − daily maximum temperature, daily minimum temperature, whether it was a snowy day or not, total daily precipitation, and daily snow precipitation. It contains data from the last 150 years for 1, 218 meteorological stations scattered over the United States. Following the preprocessing steps of Che et al. (2018b), we extract daily climate data for 100 consecutive years starting from 1910 to 2009 from 54 stations in California. To get multi-rate time series data, we split the stations into 3 groups with sampling rates of 2 days, 1 week, and 1 month respectively. We divide the data into smaller time series consisting of yearly data and end up with a dataset of 100 examples each consisting of 270 features. We perform the interpolation task on this dataset where we compute the feature values every day using the multi-rate time series data. The dataset is available for download at https://cdiac.ess-dive.lbl.gov/ftp/ushcn_daily/.
Electricity Dataset: The UCI household electricity dataset contains measurements of seven different quantities related to electricity consumption in a household. The data are recorded every minute for 47 months between December 2006 and November 2010, yielding over 2 million observations. To simulate irregular sampling, we keep observations only at durations sampled from
an exponential distribution with λ = 20. Following the preprocessing step of Binkowski et al. (2018), we also do random feature sampling where we choose one out of seven features at each time step. We divide the data into smaller time series consisting of monthly data and end up with a dataset of 1431 examples each consisting of 7 features. We perform interpolation experiments on this dataset where we compute feature values every minute using the irregularly sampled data. The dataset is available for download at https://archive.ics.uci.edu/ml/datasets/ individual+household+electric+power+consumption.
Dataset Preprocessing: We rescale time to be in [0, 1] for all datasets. We also re-scale all dimensions. In case of PhysioNet and MIMIC-III, for each dimensions we first remove outliers in the outer 0.1% percentile region. We then compute the mean and standard deviation of all observations on that dimension. The outlier detection step is used to mitigate the effect of rare large values in the data set from affecting the normalization statistics. Finally, we z-transform all of the available data (including the points identified as outliers). No data points are discarded from the data sets during the normalization process.

heading: A.6.2 SOURCE CODE
text: The source code for reproducing the results in this paper is attached with the supplementary materials. It will be made publicly available upon presentation of the paper.

heading: A.6.3 COMPUTING INFRASTRUCTURE
text: All experiments were run on a Nvidia Titan X and 1080 Ti GPUs. The time required to run all the experiments in this paper including hyperparameter tuning was approximately eight days using eight GPUs.

Review summaries:
This paper introduces a novel model, HeTVAE, for probabilistic interpolation of time series that are irregularly sampled. HeTVAE builds on prior work by complementing it with a learned time-dependent output variance in the VAE and architectural improvements. The latter include a new branch accounting for the distribution of the sampled timestamps in the series and the addition of a deterministic branch bypassing the stochastic variational latent variable. The performance of HeTVAE is evaluated on multiple datasets against various baselines and via ablation studies.

This paper introduces a VAE-based model for interpolation of irregularly sampled time series. The temporal input data is mapped to a latent representation over fixed reference points with an attention mechanism, using an intensity network that allows to encode data sparsity information. This latent representation can then be used to interpolate points at new time steps. Thanks to the intensity network and the heteroscedastic output layer, the proposed HeTVAE model can capture uncertainty estimates over the interpolated points.  The model is tested on a number of datasets containing irregularly sampled points, and outperforms competing methods in the interpolation task.

The paper introduces a novel model of Variational Autoencoder that deals with irregularly sampled time series with a probabilistic approach to do time series interpolation.  The main contribution is the architecture by itself, its components, and the training process.  The model was evaluated on both real-world data sets from the medical and climate domain and synthetic data. 

This paper introduces several improvements over the previous work mTAN to better support probabilistic interpolation. Specifically, intensity encoding is introduced to make the model be aware of information about input sparsity. Also, the homoscedastic output distribution used by previous work is replaced by a heteroscedastic distribution. Experiments results show that this improved model (HeTVAE) achieves both better likelihood estimation and mean prediction compared to previous works.