Title: SET FUNCTIONS FOR TIME SERIES

Abstract: Despite the eminent successes of deep neural networks, many architectures are often hard to transfer to irregularly-sampled and asynchronous time series that occur in many real-world datasets, such as healthcare applications. This paper proposes a novel framework for classifying irregularly sampled time series with unaligned measurements, focusing on high scalability and data efficiency. Our method SEFT (Set Functions for Time Series) is based on recent advances in differentiable set function learning, extremely parallelizable, and scales well to very large datasets and online monitoring scenarios. We extensively compare our method to competitors on multiple healthcare time series datasets and show that it performs competitively whilst significantly reducing runtime.

1. INTRODUCTION:

With the increasing digitalization, measurements over extensive time periods are becoming ubiquitous. Nevertheless, in many application domains, in particular healthcare (Yadav et al., 2018) , measurements might not necessarily be observed at a regular rate or could be misaligned. Moreover, the presence or absence of a measurement and its observation frequency may carry information of its own (Little & Rubin, 2014) , such that imputing the missing values is not always desired.

While some algorithms can be readily applied to datasets with varying length, these methods usually assume regular sampling of the data and/or require the measurements across modalities to be aligned/synchronized, preventing their application to the aforementioned settings. Existing approaches for unaligned measurements, by contrast, typically rely on imputation to obtain a regularlysampled version of a dataset for classification. Learning a suitable imputation scheme, however, requires understanding the underlying dynamics of a system; this task is significantly more complicated and not necessarily required when classification is the main goal. Furthermore, even though a decoupled imputation scheme followed by classification is generally more scalable, it may lose information (in terms of "missingness patterns") that could be crucial for prediction tasks. In addition, the fact that decoupled schemes perform worse than methods that are trained end-to-end has been has been empirically demonstrated by Li & Marlin (2016) . Approaches that jointly optimize both tasks also add a large computational overhead, thus suffering from poor scalability or high memory requirements.

Our method is motivated by the understanding that, while RNNs and similar architectures are well suited for capturing and modelling the dynamics of a time series and thus excel at tasks such as forecasting, retaining the order of an input sequence can even be a disadvantage in classification scenarios (Vinyals et al., 2015) . We show that by relaxing the condition that a sequence must be processed in order, we can naturally derive an architecture that directly accounts for (i) irregular sampling, and (ii) unsynchronized measurements. Our method SEFT: Set Functions for Time Series, extends recent advances in set function learning to irregular sampled time series classification tasks, yields state-of-the-art performance, is highly scalable and improves over current approaches by almost an order of magnitude in terms of runtime.

With SEFT, we propose to rephrase the problem of classifying time series as classifying a set of observations. We show how set functions can be exploited to learn classifiers that are naturally applicable to unaligned and irregularly sampled time series, leading to state-of-the-art performance in irregularly-sampled time series classification tasks. Our approach can be interpreted as learning dataset-specific summary statistics of time series which are optimized to separate instances by class. Furthermore, our method is highly parallelizable and can be readily extended to an online monitoring setup with up to thousands of patients.

2. RELATED WORK:

This paper focuses on classifying time series with irregular sampling and potentially unaligned measurements. We briefly discuss recent work in this field; all approaches can be broadly grouped into the following three categories.

Irregular sampling as missing data While the problem of supervised classification in the presence of missing data is closely related to irregular sampling on time series, there are some core differences. Missing data is usually defined with respect to a number of features that could be observed, whereas time series themselves can have different lengths and a "typical" number of observed values might not exist. Generally, an irregularly-sampled time series can be converted into a missing data problem by discretizing the time axis into non-overlapping intervals, and declaring intervals in which no data was sampled as missing. This approach is followed by Marlin et al. (2012) , where a Gaussian Mixture Model was used to do semi-supervised clustering on electronic health records. Similarly, Lipton et al. (2016) discretize the time series into intervals, aggregate multiple measurements within an interval, and add missingness indicators to the input of a Recurrent Neural Network. By contrast, Che et al. (2018) present several variants of the Gated Recurrent Unit (GRU) combined with imputation schemes. Most prominently, the GRU-model was extended to include a decay term (GRU-D), such that the last observed value is decayed to the empirical mean of the time series via a learnable decay term. While these approaches are applicable to irregularly-sampled data, they either rely on imputation schemes or empirical global estimates on the data distribution (our method, by contrast, requires neither), without directly exploiting the global structure of the time series.

Frameworks supporting irregular sampling Some frameworks support missing data. For example, Lu et al. (2008) directly defined a kernel on irregularly-sampled time series, permitting subsequent classification and regression with kernel-based classifiers or regression schemes. Furthermore, Gaussian Processes (Williams & Rasmussen, 2006) constitute a common probabilistic model for time series; they directly permit modelling of continuous time data using mean and covariance functions. Along these lines, Li & Marlin (2015) derived a kernel on Gaussian Process Posteriors, allowing the comparison and classification of irregularly-sampled time series using kernel-based classifiers. Nevertheless, all of these approaches still rely on separate tuning/training of the imputation method and the classifier so that structures supporting the classification could be potentially missed in the imputation step. An emerging line of research employs Hawkes processes (Hawkes, 1971; Liniger, 2009) , i.e. a specific class of self-exciting point processes, for time series modelling and forecasting (Mei & Eisner, 2017; Yang et al., 2017; Xiao et al., 2017) . While Hawkes processes exhibit extraordinary performance in these domains, there is no standardised way of using them for classification. Previous work (Lukasik et al., 2016) trains multiple Hawkes processes (one for each label) and classifies a time series by assigning it the label that maximises the respective likelihood function. Since this approach does not scale to our datasets, we were unable to perform a fair comparison. We conjecture that further research will be required to make Hawkes processes applicable to general time series classification scenarios.

End-to-end learning of imputation schemes Methods of this type are composed of two modules with separate responsibilities, namely an imputation scheme and a classifier, where both components are trained discriminatively and end-to-end using gradient-based training. Recently, Li & Marlin (2016) proposed the Gaussian Process Adapters (GP Adapters) framework, where the parameters of a Gaussian Process Kernel are trained alongside a classifier. The Gaussian Process gives rise to a fixed-size representation of the irregularly-sampled time series, making it possible to apply any differentiable classification architecture. This approach was further extended to multivariate time series by Futoma et al. (2017) using Multi-task Gaussian Processes (MGPs) (Bonilla et al., 2008) , which allow correlations between the imputed channels. Moreover, Futoma et al. (2017) made the approach more compatible with time series of different lengths by applying a Long Short Term Memory (LSTM) (Hochreiter & Schmidhuber, 1997) . We treat the j th observation as a tuple (t j , z j , m j ), comprising a time t j , a value z j , and a modality indicator m j .

All observations are summarized as a set of such tuples. Each set of tuples belonging to the same modality is then separately embedded (f ) and subsequently aggregated (agg). An attention mechanism (attn) as described in Section 3.3 is then applied to learn the importance of individual and consecutive observations. Respective query matrices for 2 attentions head are illustrated in purple and orange blocks. The results of each attention head are then concatenated and used as the input for final classification layers.

scheme, the interpolation prediction networks. It applies multiple semi-parametric interpolation schemes to obtain a regularly-sampled time series representation. The parameters of the interpolation network are trained with the classifier in an end-to-end setup.

3. PROPOSED METHOD:

Our paper focuses on the problem of time series classification of irregularly sampled and unaligned time series. We first define the required terms before describing our models

3.1. NOTATION & REQUIREMENTS:

Definition 1 (Time series). We describe a time series of an instance i as a set S i of M := len(S i ) observations s j such that S i := {s 1 , . . . , s M }. We assume each observation s j to be represented as a tuple (t j , z j , m j ), consisting of a time t j ∈ R + , an observed value z j ∈ R, and a modality indicator m j ∈ {1 . . . D}, where D represents the dimensionality of the time series. We write Ω ⊆ R + × R × N + to denote the domain of observations. An entire time series can thus be represented as

EQUATION

where for notational convenience we omitted the index i.

We leave this definition very general on purpose, allowing the length of each time series (comprising all channels, such as "heart rate", "respiratory rate", etc. of one instance) to differ, since our models are capable of handling this. Likewise, we neither enforce nor expect all time series to be synchronized, i.e. being sampled at the same time, but rather we permit unaligned or non-synchronized observations in the sense of not having to observe all modalities at each time point. Time series are collected in a dataset D.

Definition 2 (Dataset). We consider a dataset D to contain n time series. Elements of D are tuples, i.e. D := {(S 1 , y 1 ), . . . , (S N , y N )}, where S i denotes the i th time series and y i ∈ {1, . . . , C} its associated class label.

Figure 1 gives a high-level overview of our method, including the individual steps required to perform classification. To get a more intuitive grasp of these definitions, we briefly illustrate our time series notation with an example. Let instance i be an in-hospital patient, while the time series represent measurements of two channels of vital parameters during a hospital stay, namely heart rate (HR) and mean arterial blood pressure (MAP). We enumerate those channels as modalities 1 and 2. Counting from admission time, a HR of 60 and 65 beats per minute was measured after 0.5 h and 3.0 h, respectively, whereas MAP values of 80, 85, and 87 mmHg were observed after 0.5 h, 1.7 h, and 2.5 h. According to Definition 1, the time series is thus represented as S i = {(0.5, 60, 1) , (3, 65, 1) , (0.5, 80, 2) , (1.7, 85, 2) , (3, 87, 2)}. In this example, observations are ordered by modality to increase readability; in practice, we are dealing with unordered sets.

Definition 3 (Non-synchronized time series). We call a D-dimensional time series nonsynchronized if there is at least one time point t j ∈ R + at which at least one modality is not observed, i.e. if there exists

t j ∈ R + such that |{(t k , z k , m k ) | t k = t j }| = D.

Furthermore, we assume that no two measurements of the same modality m k occur at the same time, i.e. t i = t j for i = j has to be satisfied for all measurements in m k . This assumption is not required for technical reasons but for consistency. It also makes it possible to interpret the results later on.

To summarize our generic setup, we do not require M , the number of observations per time series, to be the same, i.e. len(S i ) = len(S j ) for i = j is permitted, nor do we assume that the time points and modalities of the observations are the same across time series. This setting is common in biomedical time series, for example. Since typical machine learning algorithms are designed to operate on data of a fixed dimension, novel approaches to this non-trivial problem are required.

3.2. OUR MODEL:

In the following, we describe an approach inspired by differentiable learning of functions that operate on sets (Zaheer et al., 2017; Wagstaff et al., 2019) . We phrase the problem of classifying time series on irregular grids as learning a function f on a set of arbitrarily many time series observations following Definition 1, i.e.

S = {(t 1 , z 1 , m 1 ), . . . , (t M , z M , m M )}, such that f : S → R C

, where S represents a generic time series of arbitrary cardinality and R C corresponds to the logits of the C classes in the dataset. As we previously discussed, we interpret each time series as an unordered set of measurements, where all information is conserved because the observation time is included for each set element. Specifically, we define f to be a set function, i.e. a function that operates on a set and thus has to be invariant to the ordering of the elements in the set. Multiple architectures are applicable to constructing set functions such as Transformers (Lee et al., 2019; Vaswani et al., 2017) , or Deep Sets (Zaheer et al., 2017) . Due to preliminary experiments, where Transformers suffered from lower generalization performance in our setting1 , we base this work on the framework of Zaheer et al. (2017) . Intuitively, this can be seen as computing multivariate dataset-specific summary statistics, which are optimized to maximize classification performance. Thus, we sum-decompose the set function f into the form

EQUATION

where h : Ω → R d and g : R d → R C are neural networks, d ∈ N + determines the dimensionality of the latent representation, and s j represents a single observation of the time series S. We can view the averaged representations 1/|S| sj ∈S h(s j ) in general as a dataset-specific summary statistic learned to best distinguish the class labels. Equation 2 also implies the beneficial scalability properties of our approach: each embedding can be calculated independently of the others; hence, the constant computational cost of passing a single observation through the function h is scaled by the number of observations, resulting in a runtime of O(M ) for a time series of length M .

Recently, Wagstaff et al. (2019) derived requirements for a practical universal function representation of sum-decomposable set functions, i.e the requirements necessary for a sum-decomposable function to represent an arbitrary set-function given that h and g are arbitrarily expressive. In particular, they show that a universal function representation can only be guaranteed provided that d ≥ max i len(S i ) is satisfied. During hyperparameter search we thus independently sample the dimensionality of the aggregation space, and allow it to be in the order of the number of observations that are to be expected in the dataset. Further, we explored the utilization of max, sum, and mean as alternative aggregation functions inspired by Zaheer et al. (2017) ; Garnelo et al. (2018) .

Intuition Our method can be connected to Takens's embedding theorem (Takens, 1981) for dynamical systems: we also observe a set of samples from some unknown (but deterministic) dynamical process; provided the dimensionality of our architecture is sufficiently large2 , we are capable of reconstructing the system up to diffeomorphism. The crucial difference is that we do not have to construct a time-delay embedding but rather, we let the network learn an embedding that is suitable for classification.

Time encoding In order to represent the time point of an observation on a normalized scale, we employ variant of positional encodings, as introduced by Vaswani et al. (2017) . Preliminary results indicated that this encoding scheme reduces the sensitivity towards initialization and training hyperparameters of a model. Specifically, the time encoding converts the one-dimensional time axis into a multi-dimensional input by passing the time t of each observation through multiple sine and cosine functions of varying frequencies. Given a dimensionality τ ∈ N + of the time encoding, we refer to the encoded position as x ∈ R τ , where

x 2k (t) := sin t max ts 2k/τ

(3)

x 2k+1 (t) := cos t max ts 2k/τ (4)

with k ∈ {0, . . . , τ /2} and max ts representing the maximal time scale that is expected in the data. Intuitively, we select the wavelengths using a geometric progression from 2π to max ts • 2π, and treat the number of steps and the maximum timescale max ts as hyperparameters of the model. For all experiments time encodings were used, such that an observation is represented as s j = (x (t j ) , z j , m j ).

Loss function:

If not mentioned otherwise, we choose h and g in Equation 2to be multilayer perceptron deep neural networks, parametrized by weights θ and ψ, respectively. We thus denote these neural networks by h θ and g ψ ; their parameters are shared across all instances per dataset. In our training setup, we follow Zaheer et al. (2017) and apply the devised set function to the complete time series, i.e. to the set of all observations for each time series. Overall, we optimize a loss function that is defined as

EQUATION

where (•) represents a task-specific loss function. In out setup, we either utilize the binary crossentropy in combination with a sigmoid activation function in the last layer for binary classification or multi-label classification tasks and categorical cross-entropy in combination with a softmax activation function in the last layer for multi-class classification tasks.

3.3. ATTENTION-BASED AGGREGATION:

So far, our method permits encoding sets of arbitrary sizes into a fixed-size representation. For increasingly large set sizes, however, many irrelevant observations could influence the result of the set function. The mean aggregation function is particularly susceptible to this because the influence of an observation to the embedding shrinks proportionally to the size of the set. We thus suggest to use a weighted mean in order to allow the model to decide which observations are relevant and which should be considered irrelevant. This is equivalent to computing an attention a(S, s j ) over the set input elements, and subsequently, computing the sum over all elements in the set.

Our approach is based on scaled dot-product attention with multiple heads i ∈ {1, . . . , m} in order to be able to cover different aspects of the aggregated set 3 . We define a(•), i.e. the attention weight function of an individual time series, to depend on the overall set of observations. This is achieved by computing an embedding of the set elements using a smaller set function f , and projecting the concatenation of the set representation and the individual set elements into a d-dimensional space. Specifically, we have K j,i = [f (S), s j ] T W i where W i ∈ R (im(f )+|sj |)×d and K ∈ R |S|×d . Furthermore, we define a matrix of query points Q ∈ R m×d , which allow the model to summarize different aspects of the dataset via

e j,i = K j,i • Q i √ d

and a j,i = exp(e j,i ) j exp(e j,i ) where a j,i represents the amount of attention that head i gives to set element j. The head-specific row Q i of the query matrix Q allows a head to focus on individual aspects (such as the distribution of one or multiple modalities) of a time series. For each head, we multiply the set element embeddings computed via the set function f with the attentions derived for the individual instances, i.e. r i = j a j,i f (s j ). The computed representation is concatenated and passed to the aggregation network h θ as in a regular set function, i.e. r * = [r 1 . . . r m ]. In our setup, we initialize Q with zeros, such that at the beginning of training, the attention mechanism is equivalent to computing the unweighted mean over the set elements.

Overall, this aggregation function is similar to Transformers (Vaswani et al., 2017) , but differs from them in a few key aspects. Standard Transformer blocks would use the information from all set elements in order to compute the embedding of an individual set element, leading to a runtime and space complexity of O(n 2 ). In contrast, our approach computes the embeddings of set elements independently, leading lower runtime and memory complexity of O(n). Further, we observed that computing embeddings with information from other set elements (as the Transformer does) actually decreases generalization performance (see Table 1 for details).

4. EXPERIMENTS:

We executed all experiments and implementations in a unified code base, which we also make available4 to the community. While some of the datasets used subsequently have access restrictions, anybody can gain access after satisfying the defined requirements. This ensures the reproducibility of our results. Please consult Appendix A.2 for further details.

4.1. DATASETS:

In order to benchmark the proposed method we selected 4 datasets with irregularly-sampled and non-synchronized measurements.

Healing MNIST The H-MNIST dataset was introduced by Krishnan et al. (2015) in order to simulate characteristics which typically occur in medical time series. In our setup, we use a variant of this dataset. Every instance of the dataset contains 10 frames, derived from a single instance of MNIST dataset, where the digit is rotated according to an angle uniformly sampled between -90 • to 90 • . Furthermore, 3 randomly-selected consecutive frames are augmented by a square artefact in the top left corner of the image in order to indicate seasonality in the time series. Finally, 60 % of the data points are randomly discarded in order to yield a final high-dimensional irregularly-sampled time series with non-synchronized measurements. Using these settings each instance has on average 3, 136 observations. MIMIC-III Tasks MIMIC-III (Johnson et al., 2016) is a widely-used, freely-accessible dataset containing around 50, 000 distinct ICU stays. The median length of stay is 2.1 d and a wide range of physiological measurements (e.g. arterial blood pressure, respiration rate, heart rate) are recorded with a resolution of 1 h. Furthermore, laboratory test results, collected at irregular time intervals are available. Recently, Harutyunyan et al. (2019) defined a set of machine learning tasks, labels, and benchmarks using a subset of the MIMIC-III dataset. We trained and evaluated our method and competing methods on the binary mortality prediction task (M3-Mortality) and on the multiclass problem of phenotype classification (M3-Phenotyping), while applying additional filtering described in Appendix A.1. The goal of the mortality prediction task is to predict whether a patient will die during his/her hospital stay using only data from the first 48 hours of the ICU stay. This dataset contains around 21, 000 stays of which approximately 10 % result in death. The phenotype classification task consists of 40, 000 patients, each of which can suffer from a multitude of 25 acute care conditions.

Physionet Mortality Prediction Challenge:

The 2012 Physionet challenge dataset (Goldberger et al., 2000) , which we abbreviate P-Mortality, contains 12, 000 ICU stays each of which lasts at least 48 h. For each stay, a set of general descriptors (such as gender, age, height, weight) were collected at admission time. Depending on the course of the stay and patient status, up to 37 time series variables were measured (e.g. blood pressure, lactate, respiration rate, temperature). While some modalities might be measured in regular time intervals (e.g. hourly or daily), some are only collected when required. Not all variables are available for each stay. The goal of the challenge was to predict if-and with which certainty -a patient will die during the hospital stay. The training set consists of 8, 000 stays while the testing set comprises 4, 000 ICU visits. Both datasets are similarly imbalanced, with a prevalence of around 14 %. For simplicity, the general descriptors (such as age and weight), were included as time points with a single observation at the beginning of the stay. This treatment is similar to the approach by Harutyunyan et al. (2019) in the MIMIC-III benchmarking datasets. Please refer to Table A .1, Table A.2, and Table A.3 in the appendix for a more detailed enumeration of samples sizes and label distributions. The total number of samples may slightly deviate from the originally published splits, as time series of excessive length prevented fitting some methods in reasonable time, and were therefore excluded.

4.2. COMPETITOR METHODS:

GRU-simple GRU-SIMPLE (Che et al., 2018) augments the input at time t of a Gated-Recurrent-Unit RNN with a measurement mask m d t and a δ t matrix, which contains the time since the last measurement of the corresponding modality d, such that

δ t =    s t -s t-1 + δ d t-1 t > 1, m d t-1 = 0 s t -s t-1 t > 1, m d t-1 = 1 0 t = 0

where s t represents the time associated with time step t.

Phased-LSTM The PHASED-LSTM (Neil et al., 2016) introduced a biologically inspired time dependent gating mechanism which regulates access to the hidden and cell state of a Long short-term RNN cell (Hochreiter & Schmidhuber, 1997) . While this allows the network to handle event-based sequences with irregularly spaced observations, the approach does not support unaligned measurements. In order to still provide the architecture with all relevant information, we augment the input in a similar fashion as described for the GRU-SIMPLE approach.

GRU-D GRU-D or GRU-Decay (Che et al., 2018) contains modifications to the GRU RNN cell, allowing it to decay past observations to the mean imputation of a modality using a learnable decay rate. By additionally providing the measurement masks as an input the recurrent neural network the last feed in value. Learns how fast to decay back to a mean imputation of the missing data modality.

Interpolation Prediction Networks IP-NETWORKS (Shukla & Marlin, 2019) apply multiple semiparametric interpolation schemes to irregularly-sampled time series to obtain regularly-sampled representations that cover long-term trends, transients, and also sampling information. The method combines a univariate interpolation step with a subsequent multivariate interpolation; the parameters of the interpolation network are trained with the classifier in an end-to-end fashion.

Transformer In the TRANSFORMER architecture (Vaswani et al., 2017) the elements of a sequence are encoded simultaneously and information between sequence elements is captured using Multi-Head-Attention blocks. In our case, an individual sequence element corresponds to all measurements available at a given time point, augmented with a measurement indicator. Transformers are normally used for sequence-to-sequence modelling tasks and in our setup were adapted to classification tasks by mean-aggregating the final representation. This representation is then fed into a one-layer MLP to predict logits for the individual classes. 

4.3. EXPERIMENTAL SETUP:

To permit a fair comparison between the methods, we executed hyperparameter searches for each model on each dataset, composed of uniformly sampling 20 parameters according to Appendix A.3. Training was stopped after 20 epochs without improvement of the validation loss, the hyperparameters with the best overall validation performance were selected for quantifying the performance on the test set. The train, validation, and test splits were the same for all models and all evaluations. Final performance on the test set was calculated by 3 independent runs of the models; evaluation took place after the model was restored to the state with the best validation loss. In all subsequent benchmarks, we use the standard deviation of the test performance of these runs as generalization performance estimates.

4.4. RESULTS:

The results are shown in Table 1 . Overall, our proposed method exhibits the lowest per-epoch runtime on most datasets, while either yielding competitive or state-the-art performance. Further, the trade-off between runtime and performance of the proposed method is very good on all datasets (see Figure A.1 and Figure A.2 in the appendix for a visualization of this argument). In order to elucidate the contribution of individual model components, we also provide an ablation study in Table A .4. Here we see that the attention mechanism contributes more to the model performance, while the positional encoding seems to be beneficial for datasets with highly-varying time series lengths, in particular M3-Phenotyping.

Opening the black box In the medical domain, it is of particular interest to understand the decisions a model makes based on the input it is provided with. The formulation of our model and its per observation perspective on time series gives it the unique property of being able to quantify to which extent an individual observation contributed to the output of the model. We exemplify this in Figure 2 with a patient time series that was combined with our models attention values, displayed for a set of clinically relevant variables. After reviewing these records with our medical expert, we find that our model is able to pick up regions with drastic changes in individual modalities. Moreover, it is able to inspect other modalities at the same associated time (for instance, at hour 20). This is behaviour similar to what one would expect from an alerted clinician reviewing the logged medical records. Interestingly, we observe that the model attends to known trends (that are consisting with domain knowledge about patient deterioration ultimately resulting in death) such as increase in lactate or hemodynamic instability, as indicated by drops in blood pressure. Furthermore, the model appears to be alerted by persisting low urine output. After several hours, this can be indicative of kidney failure.

Table 1 : Performance comparison of methods on benchmarking datasets. Performance metrics have been rescaled to 100 for readability reasons. "AUC" denotes the area under the Receiver Operating Characteristic (ROC) curve; "PR AUC" denotes the area under the precision recall curve. "MICRO" refers to evaluating the metric globally by treating each entry of the label indicator matrix as a label. For "MACRO", the metric is computed for each class and then averaged, whereas in "WEIGHTED" the class-wise metrics are weighted by class imbalance. Values denoted with "OOM" were not obtainable due to restrictions in GPU memory. 

5. CONCLUSION:

In this work, we presented a novel approach for classifying time series with irregularly-sampled and unaligned, that is non-synchronized, observations. Our approach yields state-of-the-art to strongly competitive performance on numerous simulated and real-world datasets, while reducing runtime by almost half. Moreover, we demonstrated that combining the perspective of individual observations with an attention mechanism permits increasing the interpretability of the model. This is particularly relevant for the medical and healthcare applications.

For future work, we reserve a more extensive exploration of the learned latent representation to evaluate its utility for clustering of time series or visualization of their similarity. Further we use a constant architecture for the attention network f with 2 layers, 64 nodes per layer, 4 heads and a dimensionality of the dot product space d of 128. We solely sample the amount of attention dropout uniformly from the values (0.0, 0.1, 0.25, 0.5).

Transformer We utilize the same model architecture as defined in Vaswani et al. (2017) , where we use a one hidden layer MLP as a feed-forward network, with dimensionality of the hidden layer selected to be twice the model dimensionality. The parameters for the Transformer network were sampled according to the following criteria. The dimensionality of the model was sampled uniformly from the values (64, 128, 256, 512, 1024) , the number of attention heads per layer from the values (2, 4, 8) and the number of layers from the range [1, 6] ∈ N. Further, we sampled the amount of dropout of the residual connections and the amount of attention dropout uniformly from the values (0.0, 0.1, 0.2, 0.3, 0.5), and the maximal timescale for the time embedding from the values (10, 100, 1000) (similar to the SEFT approach). 1 for more details.

:

Please see Section 4.4 for a quantification.

In Takens's embedding theorem, d > dB is required, where dB refers to the fractal box counting dimension(Liebovitch & Toth, 1989), which is typically well below the size of typical neural network architectures.

3 Since we are dealing only with a single instance (i.e. time series) in this section, we use i and j to denote a head and an observation, respectively.

https://osf.io/2hg74/?view_only=8d45fdf237954948a02f1e2bf701cdf1

Figures and Tables:

Figure 2: Visualizations of a single attention head on an instance of the P-Mortality dataset. We display a set of variables relevant for assessing patient stability and organ failure: Serum Potassium (K), Lactate, Systolic Arterial Blood Pressure (SysABP), and Urine output. Darker colors represent higher attention values.

Figure A.1: A visualisation of the runtime of all methods and their AUC for datasets with a binary classification scenario.

Due to the high dimensionality of H-MNIST and associated memory issues, the set elements were constructed by concatenating the observation time with all values associated with the time point and measurement indicators. Furthermore, as this dataset features only 10 time steps and missingness is induced randomly, we refrained from applying the attention-based aggregation.

Table content: <table><tr><td>DATASET</td><td>MODEL</td><td colspan="3">MICRO AUC MACRO AUC WEIGHTED AUC</td><td>RUNTIME</td></tr><tr><td/><td>GRU-SIMPLE</td><td colspan="3">99.09 ± 0.05 99.01 ± 0.05 99.03 ± 0.05</td><td>11.43 ± 0.47</td></tr><tr><td/><td colspan="4">PHASED-LSTM 98.63 ± 0.13 98.50 ± 0.15 98.52 ± 0.14</td><td>33.93 ± 1.11</td></tr><tr><td>H-MNIST</td><td>GRU-D IP-NETS</td><td colspan="3">99.42 ± 0.01 99.37 ± 0.02 99.38 ± 0.02 99.06 ± 0.05 98.96 ± 0.03 98.98 ± 0.03</td><td>11.81 ± 0.44 127.76 ± 0.95</td></tr><tr><td/><td colspan="4">TRANSFORMER 99.59 ± 0.05 99.55 ± 0.05 99.56 ± 0.05</td><td>21.62 ± 0.90</td></tr><tr><td/><td>SEFT*</td><td colspan="3">99.76 ± 0.01 99.75 ± 0.01 99.75 ± 0.01</td><td>4.05 ± 0.35</td></tr><tr><td/><td>GRU-SIMPLE</td><td colspan="3">79.89 ± 0.14 73.91 ± 0.19 72.55 ± 0.16</td><td>112.58 ± 2.03</td></tr><tr><td/><td colspan="4">PHASED-LSTM 80.00 ± 0.06 73.91 ± 0.09 72.65 ± 0.08</td><td>400.41 ± 14.14</td></tr><tr><td/><td>GRU-D</td><td colspan="3">82.16 ± 0.04 77.14 ± 0.03 76.08 ± 0.01</td><td>288.70 ± 16.66</td></tr><tr><td>M3-Phenotyping</td><td>IP-NETS</td><td>-OOM-</td><td>-OOM-</td><td>-OOM-</td><td>-OOM-</td></tr><tr><td/><td>TRANSFORMER</td><td>-OOM-</td><td>-OOM-</td><td>-OOM-</td><td>-OOM-</td></tr><tr><td/><td>SEFT</td><td colspan="3">81.22 ± 0.12 75.95 ± 0.09 74.90 ± 0.11</td><td>56.27 ± 2.14</td></tr><tr><td/><td>SEFT-ATTN</td><td colspan="3">82.00 ± 0.06 76.95 ± 0.09 75.88 ± 0.09</td><td>52.32 ± 0.74</td></tr><tr><td/><td/><td>ACCURACY</td><td>PR AUC</td><td>AUC</td><td>RUNTIME</td></tr><tr><td/><td>GRU-SIMPLE</td><td colspan="3">88.24 ± 0.38 36.36 ± 1.31 79.36 ± 0.26</td><td>22.80 ± 0.56</td></tr><tr><td/><td colspan="4">PHASED-LSTM 88.32 ± 0.31 35.30 ± 1.38 80.16 ± 0.22</td><td>25.54 ± 0.26</td></tr><tr><td/><td>GRU-D</td><td colspan="3">89.56 ± 0.38 46.76 ± 0.65 83.73 ± 0.21</td><td>31.85 ± 0.86</td></tr><tr><td>M3-Mortality</td><td>IP-NETS</td><td colspan="3">89.73 ± 0.16 45.88 ± 0.87 83.30 ± 0.56</td><td>101.12 ± 4.52</td></tr><tr><td/><td colspan="4">TRANSFORMER 89.14 ± 0.15 42.32 ± 0.41 82.60 ± 0.55</td><td>4.79 ± 0.02</td></tr><tr><td/><td>SEFT</td><td colspan="3">88.65 ± 0.49 36.18 ± 5.07 79.15 ± 3.00</td><td>3.72 ± 0.11</td></tr><tr><td/><td>SEFT-ATTN</td><td colspan="3">89.48 ± 0.16 45.25 ± 0.96 83.79 ± 0.59</td><td>16.64 ± 0.20</td></tr><tr><td/><td>GRU-SIMPLE</td><td colspan="3">85.66 ± 0.14 39.43 ± 0.71 79.79 ± 0.16</td><td>5.16 ± 0.06</td></tr><tr><td/><td colspan="4">PHASED-LSTM 85.57 ± 0.11 39.55 ± 0.62 78.71 ± 0.76</td><td>18.59 ± 1.15</td></tr><tr><td/><td>GRU-D</td><td colspan="3">87.19 ± 0.30 54.95 ± 0.54 86.58 ± 0.32</td><td>14.08 ± 0.38</td></tr><tr><td>P-Mortality</td><td>IP-NETS</td><td colspan="3">87.23 ± 0.18 54.87 ± 0.41 86.42 ± 0.18</td><td>7.21 ± 0.46</td></tr><tr><td/><td colspan="4">TRANSFORMER 86.47 ± 0.08 48.72 ± 0.61 83.49 ± 0.46</td><td>2.69 ± 0.43</td></tr><tr><td/><td>SEFT</td><td colspan="3">87.11 ± 0.32 52.07 ± 0.41 84.12 ± 0.32</td><td>3.07 ± 0.03</td></tr><tr><td/><td>SEFT-ATTN</td><td colspan="3">87.62 ± 0.16 54.05 ± 0.27 85.50 ± 0.13</td><td>7.54 ± 0.08</td></tr></table>

3: M3-Phenotyping prevalence of labels for the multi label classification task

Table content: <table><tr><td>Phenotype</td><td colspan="3">Training Validation Testing</td></tr><tr><td>Acute and unspecified renal failure</td><td>0.216</td><td>0.207</td><td>0.211</td></tr><tr><td>Acute cerebrovascular disease</td><td>0.0746</td><td>0.0753</td><td>0.0662</td></tr><tr><td>Acute myocardial infarction</td><td>0.103</td><td>0.103</td><td>0.108</td></tr><tr><td>Cardiac dysrhythmias</td><td>0.322</td><td>0.317</td><td>0.323</td></tr><tr><td>Chronic kidney disease</td><td>0.135</td><td>0.131</td><td>0.132</td></tr><tr><td>Chronic obstructive pulmonary disease</td><td>0.132</td><td>0.128</td><td>0.126</td></tr><tr><td>and bronchiectasis</td><td/><td/><td/></tr><tr><td>Complications of surgical procedures or</td><td>0.207</td><td>0.201</td><td>0.213</td></tr><tr><td>medical care</td><td/><td/><td/></tr><tr><td>Conduction disorders</td><td>0.0726</td><td>0.07</td><td>0.0704</td></tr><tr><td>Congestive heart failure; nonhyperten-</td><td>0.268</td><td>0.264</td><td>0.268</td></tr><tr><td>sive</td><td/><td/><td/></tr><tr><td>Coronary atherosclerosis and other</td><td>0.323</td><td>0.317</td><td>0.331</td></tr><tr><td>heart disease</td><td/><td/><td/></tr><tr><td>Diabetes mellitus with complications</td><td>0.0955</td><td>0.0945</td><td>0.094</td></tr><tr><td>Diabetes mellitus without complication</td><td>0.194</td><td>0.187</td><td>0.192</td></tr><tr><td>Disorders of lipid metabolism</td><td>0.291</td><td>0.287</td><td>0.289</td></tr><tr><td>Essential hypertension</td><td>0.421</td><td>0.41</td><td>0.424</td></tr><tr><td>Fluid and electrolyte disorders</td><td>0.267</td><td>0.276</td><td>0.265</td></tr><tr><td>Gastrointestinal hemorrhage</td><td>0.0715</td><td>0.0747</td><td>0.0788</td></tr><tr><td>Hypertension with complications and</td><td>0.133</td><td>0.131</td><td>0.13</td></tr><tr><td>secondary hypertension</td><td/><td/><td/></tr><tr><td>Other liver diseases</td><td>0.0884</td><td>0.0904</td><td>0.0883</td></tr><tr><td>Other lower respiratory disease</td><td>0.0514</td><td>0.0484</td><td>0.0565</td></tr><tr><td>Other upper respiratory disease</td><td>0.0408</td><td>0.0371</td><td>0.0429</td></tr><tr><td>Pleurisy; pneumothorax; pulmonary</td><td>0.0858</td><td>0.09</td><td>0.0905</td></tr><tr><td>collapse</td><td/><td/><td/></tr><tr><td>Pneumonia (except that caused by tu-</td><td>0.14</td><td>0.135</td><td>0.135</td></tr><tr><td>berculosis or sexually transmitted dis-</td><td/><td/><td/></tr><tr><td>ease)</td><td/><td/><td/></tr><tr><td>Respiratory failure; insufficiency; arrest</td><td>0.18</td><td>0.184</td><td>0.177</td></tr><tr><td>(adult)</td><td/><td/><td/></tr><tr><td>Septicemia (except in labor)</td><td>0.142</td><td>0.145</td><td>0.138</td></tr><tr><td>Shock</td><td>0.0783</td><td>0.0745</td><td>0.0811</td></tr><tr><td>Total samples</td><td>29 208</td><td>6359</td><td>6266</td></tr></table>

4: Ablation study of individual components of SEFT. "AUC" denotes the area under the Receiver Operating Characteristic (ROC) curve; "PR AUC" denotes the area under the precision recall curve; "RUNTIME" denotes the runtime of one training epoch in seconds. SEFT (NO ATTENTION) 81.22 ± 0.12 75.95 ± 0.09 74.90 ± 0.11 56.27 ± 2.14 SEFT-ATTN (NO TIME ENC.) 80.46 ± 0.86 74.70 ± 1.12 73.48 ± 1.18 50.17 ± 0.84 SEFT-ATTN 82.00 ± 0.06 76.95 ± 0.09 75.88 ± 0.09 52.32 ± 0.74 NO ATTENTION) 88.65 ± 0.49 36.18 ± 5.07 79.15 ± 3.00 3.72 ± 0.11 SEFT-ATTN (NO TIME ENC.) 89.31 ± 0.08 44.12 ± 0.06 83.72 ± 0.34 17.60 ± 0.43 SEFT-ATTN 89.48 ± 0.16 45.25 ± 0.96 83.79 ± 0.59 16.64 ± 0.20 P-Mortality SEFT (NO ATTENTION) 87.11 ± 0.32 52.07 ± 0.41 84.12 ± 0.32 3.07 ± 0.03 SEFT-ATTN (NO TIME ENC.) 87.03 ± 0.06 51.86 ± 1.04 84.91 ± 0.29 7.04 ± 0.04 SEFT-ATTN 87.62 ± 0.16 54.05 ± 0.27 85.50 ± 0.13 7.54 ± 0.08

Table content: <table><tr><td>DATASET</td><td>MODEL</td><td colspan="2">MICRO AUC MACRO AUC WEIGHTED AUC</td><td>RUNTIME</td></tr><tr><td>H-MNIST</td><td>SEFT</td><td colspan="2">99.76 ± 0.01 99.75 ± 0.01 99.75 ± 0.01</td><td>4.05 ± 0.35</td></tr><tr><td>M3-Phenotyping</td><td/><td/><td/></tr><tr><td/><td/><td>ACCURACY PR AUC</td><td>AUC</td><td>RUNTIME</td></tr><tr><td>M3-Mortality</td><td>SEFT (</td><td/><td/></tr></table>

Peer review (iclr2020):
Summary: The work is focused on classification of irregularly sampled and unaligned multi-modal time series. Prior work has primarily focused on imputation methods, either end-to-end or otherwise. This paper approaches the problem as a set function mapping between the time-series tuples to the class label. The proposed method is uses a set encoding of a multi-modal time series input, followed by mode-specific encoding of the tuples which are then aggregated in multiple ways prior to classification. An attention mechanism is attached in order for the model to automatically weigh the relevance of tuples for the classification. The model is compared to imputation based baselines on clinical ICU time series classification tasks. The performance mostly appears comparable across baselines but the proposed method has much better run-times. The paper is for the most part well written, and related work well characterized. The formulation is interesting and clinically relevant as well so the choice of data-sets makes some sense. I have a few concerns about the architecture formulation and lack of clarification and intuition in what appears to be the main contribution of the paper (Sec 3.2 and 3.3) which I will detail below: a. In the evaluation, I really want to see a decoupling between the *time encoding step* and *attention based aggregation* on the performance to figure out to isolate different sources of performance improvements. That is can there be a SEFT without time encoding? If not, why not? I encourage more ablation like studies that look at different sources of performance gains and demonstrate them in experiments. b. The description of Sec 3.3. is really missing key motivation for the choices made around how the attention formulation is designed. For example why does the dot produce include the set elements? What if it doesn*t? What is Q supposed to capture? c. Is a_{j,i} shared across instances? Then irrespective of the number of observations per instance, the tuple gets similar weights? If not appropriate indexing will help clarify this. d. It would be useful to provide how exactly a label is inferred for a *new* test instance. I have some minor additional feedback (just for presentation and motivation purposes): 1. Authors make a claim in the introduction which should likely be qualified with a citation - *Furthermore, even though a decoupled imputation scheme followed by classification is generally more scalable, it may lose information that is relevant for prediction tasks*. How does decoupled imputation imply loss of relevant information? By losing information about which observations are missing and relying on that for prediction? Does this clinically make sense? Or even generally? 2. In Sec 3.3, you probably mean . That is parenthesis are missing? 3. What are the +- std errors indicating? Is it cross validation error on a held-out test set? 4. Initially is indexing samples and by equation (3), (4) indexes time(?) and in Sec 3.3 indexes observations? How are observations defined here? is it measurement of specific modality at a specific time instance? Can you clear this in the introduction itself? ----------------------------------------------------------------------------------------------------------------------------------------------------------------------- I have read the authors updated draft and response. The experiments section looks much better now. 1. The overall contribution has less clinical utility in my opinion as generally a patient likely deteriorates over time before an adverse outcome and therefore -- to give the model too much flexibility w.r.t. time ordering doesn*t make quite as much sense. This is reflected in the fact that experimental results are not drastically better than other baselines. The authors might be able to show the utility of the method on other time series classification datasets where this is not a limitation of the data itself. However in those settings, it may be a bit hard to beat transformers. Do the authors have a sense of where the benefits of this method really are? 2. Mortality tasks are generally on the simpler side of clinical prediction problems as well. Nonetheless I think the contribution has some utility to the community. I do encourage the authors to try non--clinical datasets for a comparison 3. Please have a discussion that includes limitations and to discuss where the benefits of your methods really lie. A clear and thoughtful discussion is currently missing in your conclusions. With that said, I am updating my score to a 6.