Year: 2021

Abstract:
We introduce Latent Temporal Flows (LatTe-Flows), a method for probabilistic multivariate time-series analysis tailored for high dimensional systems whose temporal dynamics are driven by variations in a lower-dimensional discriminative subspace. We perform indirect learning from hidden traits of observed sequences by assuming that the random vector representing the data is generated from an unobserved low-dimensional latent vector. LatTe-Flows jointly learns auto-encoder mappings to a latent space and learns the temporal distribution of lower-dimensional embeddings of input sequences. Since encoder networks retain only the essential information to generate a latent manifold, the temporal distribution transitions can be more efficiently uncovered by time conditioned Normalizing Flows. The learned latent effects can then be directly transferred into the observed space through the decoder network. We demonstrate that the proposed method significantly outperforms the state-of-the-art on multi-step forecasting benchmarks, while enjoying reduced computational complexity on several real-world datasets. We apply LatTe-Flows to a challenging sensor-signal forecasting task, using multivariate time-series measurements collected by wearable devices, an increasingly relevant health application.

heading: 1 INTRODUCTION
text: One of the core objectives in machine learning research is to build models that accurately capture and explain complex structures in real-world systems. Much of unsupervised learning is driven by probabilistic generative modeling, which is notoriously difficult in high dimensions as a consequence of “the curse of dimensionality.” For example, biomarker observations (e.g. heart rate and respiratory rate) collected from wrist-worn devices are naturally collected sequentially from a large number of correlated features.
Patterns and complex interactions in sensor signals have been shown to be associated with changes in physiology (Goodale et al., 2019). Recorded time-series from different sensors can form multivariate time-series data. Joint probability modeling of these multivariate sensors across time allows for analysis of expected or abnormal patterns in physiological vital signs. Expressive time-series generative models can be used to perform sequential data analysis and forecasting. Specifically, these models allow for long-term predictions with uncertainty estimation (Futoma et al., 2016), time-series data augmentation (Wen et al., 2020), out-of-distribution detection (Serrà et al., 2019), missing data imputation (Luo et al., 2018), and improved model interpretability analyses (Ismail et al., 2020; Rooke et al.). Such information can be integrated into health monitoring systems (Dunn et al., 2021).
Given a sequence of data of dimensionality N , with N often large, a multivariate time-series forecasting model aims to effectively predict the future outcomes of each variable n ∈ {1, . . . , N}, given their past values. Typically, this process requires an accurate estimate of a multivariate conditional distribution of temporal transitions, making the task challenging in the high-dimensional setting. This task is challenging in terms of computational complexity and expressiveness due to the number of parameters required to be estimated for the multivariate parameterization. Thus, directly estimating the full predictive distribution of the future values in the high-dimensional observation space is sample-inefficient and may require a large amount of training data. The intractability of estimating such models for large systems, due to the growth in the number of parameters with di-
mension, has limited existing methods to handle at most a few dimensions or requires restrictive assumptions such as tractable distribution classes (Chung et al., 2015) or low-rank approximations (Salinas et al., 2019a). An alternative way to tackle time-series generative modeling is to view the collection of data as N separate time sequences and fit a separate model for each of the dimensions. However, univariate techniques (Zhang et al., 2017; Oreshkin et al., 2019; Montero-Manso et al., 2020), do not benefit from joint leaning of temporal dependencies between features and are thus limited in forecasting accuracy. Since many multivariate time-series in practical applications are highly correlated, it is crucial to learn both intra-series and inter-series patterns. Modeling dependency relationships among the variables between the individual time-series, including indirect relationships through shared latent causes (Wang et al., 2019), shows promise in enhancing the performance of forecasting.
Data dependencies may be hidden in high dimensions due to noise and high variance, and thus can be missed when conducting direct data analysis in the higher dimensional space. In such cases, we aim to learn dynamics in a compact latent space to enable faster forecasting, and potentially reveal latent trends. We assume that the observed, possibly high-dimensional, random vector representing the data of interest, is generated from an unobserved low-dimensional latent vector through a timevarying probabilistic process (Zhang et al., 2019), (Laumer et al., 2020), (Louis et al., 2019), (Chung et al., 2015), (Yoon et al., 2019), (Nguyen & Quanz, 2021).
By modeling a lower-dimensional latent vector by a (stationary) reversible embedding mapping, LatTe-Flows increases parsimony, while reducing computational complexity. These mappings are trained to favor latent representations that uncover latent dynamics over future latent states. Under LatTe-Flows, forecasts are generated by estimating the latent state evolution, along with recovery functions. Nonlinear mappings are trained via deep auto-encoder networks, while latent transition dynamics are learned using a combined structure of Normalizing Flows (Rasul et al., 2020; Dinh et al., 2016; Papamakarios et al., 2017) and a multivariate Recurrent neural network (RNN) (Graves, 2013; Sutskever et al., 2014).
Our extensive experiments on real-world datasets illustrate LatTe-Flows’ state-of-the-art performance and computational tractability on the majority of the datasets considered. To demonstrate the ability of the proposed model to match the ground truth trajectory distribution, we apply LatTe-Flows on the Apple Heart and Movement Study (AH&MS) dataset for a challenging vital sensor-signal forecasting task. We visualize the 2-dimensional and 3-dimensional learned latent representations from our model to showcase the interpretable low-dimensional embeddings. These emebeddings show that LatTe-Flows is able to easily identify participants’ V O2max, a main indicator and summary of cardiorespiratory fitness, while only being trained on relatively lower-level vital signals like resting heart rate and heart rate variability.

heading: 2 RELATED WORK AND OUR CONTRIBUTIONS
text: The literature on time-series forecasting has a long history—in this work, we mainly focus on recent developments in the deep learning context. Simple models for multivariate data, such as general State-Space Models (SSMs) (Hamilton, 2020; Liu et al., 2016), N-BEATS (Oreshkin et al., 2019), Gaussian Processes (Rasmussen, 2003), DeepAR (Salinas et al., 2020; Zhang et al., 2017) learn one model per time-series (univariate methods). As a consequence, they cannot effectively capture complex structure and interdependencies between multiple time-series.
Although multivariate probabilistic time-series forecasting models estimate the full predictive distribution, the number of parameters in these models grows quickly with the number of variables, which results in large computational cost and a high risk of overfitting in high dimensional settings. Methods such as Variational Recurrent Neural Networks (VRNNs) (Chung et al., 2015) or TimeGAN (Yoon et al., 2019), either assume pre-selected tractable distribution classes or another type of structural approximation (Salinas et al., 2019a). In the low-rank Gaussian copula model, for instance (Salinas et al., 2019a), a multitask univariate LSTM (Hochreiter & Schmidhuber, 1997b) is used to output transformed time-series and diagonal and low-rank factors of a Gaussian covariance matrix. These assumptions can limit the distributional expressiveness of low-rank Gaussian copula models.
Recently, more flexible models such as Temporal Conditioned Normalizing Flows (Rasul et al., 2020) have been proposed. This method uses a multivariate RNN to learn temporal dynamics with the state translated to the output joint distribution via a Normalizing Flow (Dinh et al., 2016). However, during forecasting, an invertible flow is applied on the same number of latent dimensions as input dimensions, thus it does not scale to large numbers of time-series (since RNNs have quadratic complexity in N ). The work most closely related to ours is TLAE (Nguyen & Quanz, 2021), where the temporal model is applied across a low dimensional space. This method combines an RNN-based model with auto-encoders to learn a temporal deep learning latent space forecast model. TLAE focuses on the encoder/decoder modeling capabilities and assumes simplistic probabilistic structure on the latent vector (multivariate Gaussian with diagonal covariance matrix), which can be restrictive. We stress that the encoder mapping for TLAE is time-dependent, along with the probabilistic model.
Our key contributions in this work can be summarized as follows:
• LatTe-Flows can model cross-series dependencies across time while scaling to large dimensions. We capture the joint temporal cross-series dynamics, by pairing non-linear dimensionality reduction with latent temporal distribution estimation, while scaling to a large number of time-series in RN as a consequence of reduced dimensionality of the (joint model) learning space.
• Forecasting in LatTe-Flows is performed across a low dimensional space, enabling faster sequence generation. Instead of applying the temporal model across all series, multivariate forecasts in LatTe-Flows are generated across a low dimensional space, and then mapped back to the observation space.
• We allow flexibly and expressive probabilistic models in the latent space. To increase model flexibility, we introduce explicit learning of the latent conditional distribution of temporal transitions without strong assumptions of traditional multivariate models.
• We introduce a flexible end-to-end training process. We are able to harness stochastic gradient descent by combining the objectives for sequence reconstruction and latent probabilistic prediction.

heading: 3 PROBLEM STATEMENT AND PROPOSED APPROACH
text: Consider a collection of high dimensional multivariate time-series yn,t ∈ R, where n ∈ {1, 2, . . . , N} indexes the individual univariate component time-series, and t indexes time. Consequently, the multivariate observation vector at time t is given by yt ∈ RN . Given a sequence of Ttotal vector realizations of yt, a multivariate time-series can be represented as a matrix Y, Y ∈ RN×Ttotal . We focus on the task of multivariate time-series multi-step forecasting. More formally, let us assume that we are given an observed history ( y1, . . . ,yT ) , sampled from the complete time-series history of the training data, where each instance consists of N temporal features (that occur over time, e.g. vital sensor signals). Our goal is to learn future values of the series over a length-τ, τ > 0, forecast horizon and predict a set of plausible future trajectories ( ŷT+1, . . . , ŷT+τ ) by learning the conditional distribution p(yT+1:T+τ |y1:T ) of temporal transitions. We refer to time-series yn,1:T+τ as the target time-series, which for training is split according to a time range ( 1, 2, . . . , T ) referred
to as context window, and to time ( T + 1, T + 2, . . . , T + τ ) as prediction horizon.
The temporal dynamics of complex systems are often driven by fewer and lower-dimensional factors of variation (Laumer et al., 2020), (Louis et al., 2019), (Chung et al., 2015), (Yoon et al., 2019), (Nguyen & Quanz, 2021), (Amiridi et al., 2021). We assume that the observed, possibly high-dimensional, random vector yt ∈ RN representing the data of interest is generated from an unobserved low-dimensional latent vector xt ∈ RD with D N through a time varying probabilistic process p(xt|x1:t−1). To reduce the computational burden while at the same time improving distribution modeling, we propose incorporating representation learning in the generative learning problem to explicitly learn the temporal distribution of compact representations of input sequences. Specifically, LatTe-Flows consists of three key components: an embedding function g : RN → RD, a latent conditional distribution p(xt|x1:t−1) of temporal transitions, and a recovery function q : RD → RN . During training, the model simultaneously learns to produce sequence
representations that will push a latent temporal generative model to predict plausible latent future sequences, recover them back into the observed space, and iterate across time.
Given a model for p(yT+1:T+τ |y1:T ), one can estimate the conditional expectation, which can be expressed as a function of past observations: E[yT+1:T+τ |y1:T ] = f(y1, . . . ,yT ). An indirect estimate of this function can be accomplished using the following strategy: during testing, future trajectories ( ŷT+1, . . . , ŷT+τ ) are generated by embedding the past history via g : ( x1 =
g(y1), . . . ,xT = g(yT ) ) , sampling from the latent future distribution ( x̂T+1, . . . , x̂T+τ ) ∼
p(·|x1:T ), followed by applying a non-linear recovery function q: ( ŷT+1 = q(x̂T+1), . . . , ŷT+τ} =
q(x̂T+τ ) ) . Next, we describe each building block of our approach and the combined training strategy.

heading: 3.1 DIMENSIONALITY REDUCTION OF HIGH-DIMENSIONAL SEQUENCES
text: Learning Sequence Representations
The goal of the proposed method is to obtain suitable representations xt ∈ RD that reveal a reduced search space for future sequence forecasting, where underlying patterns and meaningful information among features is preserved. Latent-space forecasting has been considered both in (Yu et al., 2016) as a result of matrix factorization as well as in (Nguyen & Quanz, 2021). In Yu et al. (2016), a multivariate time-series Y is decomposed into components Q ∈ RN×D and X ∈ RD×Ttotal , with X temporally constrained. Let us denote G as the pseudo-inverse of Q. If Y can be decomposed by G and X, forecasting for the high-dimensional series Y can be performed on a smaller dimensional series. In TLAE (Nguyen & Quanz, 2021), linear mappings are generalized to nonlinear transformations via time-varying auto-encoders. In our framework, each input vector yt ∈ RN is mapped to a condensed representation xt ∈ RD (usually D N ) using a stationary reversible embedding mapping g : RN → RD. Such representations are trained to reveal a lowdimensional structure, which allows an expressive family of temporal constrained distributions to fully uncover this, i.e., mappings that yield a low negative log-likelihood cost over p(xt|x1:t−1).
The first component of our model is estimating the dimensionality reducing mapping. Exploiting the deep neural network’s ability to to extract higher order features (Yang et al., 2017) and approximate any nonlinear function, we replace G by an encoder and Q by a decoder neural network. An encoder network gφ : RN → RD, embeds yt at time t into a low-dimensional latent space, vector xt, using a nonlinear map. Operating within the latent space, we then seek, a dynamical system that prescribes a rule to move forward in time, and a decoder network qφ′ : RD :→ RN , ŷt = qφ′(xt) to reconstruct latent variables in the spatial domain. Although the embedding mapping is assumed to be stationary, latent representations’ progressions over time are captured via an autoregressive deep learning model, where the data distribution is represented by a conditioned Normalizing Flow. During training, the auto-encoder learns by fine-tuning the parameters of a feed-forward Deep Neural Network (DNN) in such a way that the reconstruction error is minimized when back projected with another feed-forward DNN. These networks need to be specified a-priori, in terms of the number of layers and neurons. We note that the specific requirements for g and q are problem dependent, and we detail the particular design we use in the Appendix. Given a batch of time-series B, and with |B| denoting the cardinality of the batch set, the first term of our overall cost function consists of the reconstruction loss: LREC := 1|B|(T+τ) ∑ y1:T+τ∈B ∑T+τ t=1 ||yt − q(g(yt;φ);φ′)|| 2 2.

heading: 3.2 LATENT SEQUENCE MODELING
text: The latent random vector lies at the heart of our overall probabilistic model: It is assumed to “encode” the observed data in a compact manner through xt = gφ(yt), allowing an accurate model of the probabilistic model pθ(xt|x1:t−1), from which new data can be generated. Our goal is to recover the latent space structure through a flexible time-conditioned distribution model in which the most important features are kept. We also learn the mapping ŷt = qφ′(xt) that translates the latent effects to the original data space.
Using the chain rule, the joint distribution of predicted values conditioned on observed values, p(yT+1, . . . ,yT+τ |y1:T ), can be written as a product of conditional distributions. Autoregressive models use a neural network to approximate the conditional distribution p(yt|y1:t−1) by a paramet-
ric distribution pθ(yt|y1:t−1) specified by learnable parameters θ. The prediction at time t is input to the model to predict the value at time t+1: p(yT+1, . . . ,yT+τ |y1:T ) = ∏T+τ t=T+1 p(yt|y1:t−1). We wish to replace this decomposition by a tractable, approximate statistic ht of the past. To represent the history of observations in a compressed state vector ht we use RNNs (Graves, 2013; Sutskever et al., 2014), with the most well-known variants, the LSTM (Hochreiter & Schmidhuber, 1997a) and GRU (Chung et al., 2014).
We assume that a time-related (and recursively updated) vector ht ∈ RH , can summarize the history of the time-series up to time time point t, ht = RNNθ(yt−1,ht−1), where RNNθ is a multilayer LSTM or GRU parameterized by shared weights θ and h0 = 0. The state is compressed as it uses less space than the history of observations. Under this model, we can factorize the joint distribution of the observations as (yT+1, . . . ,yT+τ |y1:T ) = ∏T+τ t=T+1 pθ(yt|ht), where now θ comprises both the weights of the RNN as well as the probabilistic model. This model is autoregressive as it consumes the observations at the time step t − 1 as input to learn the distribution of the next time step t. Then, time conditioning on the latent generative model pθ(xt|ht) can be realized by employing a multivariate RNN to model the series progressions, with the state translated to the output joint distribution via a flow (we focus on Real-NVP, but MAF is also explored in our experiments; background on Normalizing Flows can be found in the Appendix ). This combination retains the power of autoregressive models—such as good performance in extrapolation into the future—with the flexibility of flows as an expressive distribution model.
Following Rasul et al. (2020), we concatenate h to the inputs of the scaling and translation function approximators of the coupling layers, i.e. s(concat(x1:d,h)) and t(concat(x1:d,h)). The model, which is parameterized by both the flow (the weights of the scaling and translation neural networks) and the RNN – θ is trained by minimizing LNegLL := 1|B|(T+τ) ∑ y1:T+τ∈B ∑T+τ t=T+1 log p(yt|ht; θ).

heading: 3.3 THE TRAINING AND TESTING PROCESS
text: A schematic overview of the training procedure is depicted in Figure 1. Each time-series Y is split into a training YTR and test set YTS by using all data prior to a fixed date for training and using rolling windows for the test set. A batch input YB ∈ RD×(T+τ) defines a sub-matrix of YTR with column indices defined by the set B – it contains two components: the first part {yt}Tt=1 is associated with the past input, while the second component {yt}T+τt=T+1 is associated with future observations. A batch of time-series YB is embedded into latent variables via g : xt = g(yt), yielding a matrix of latent codes XB ∈ RD×(T+τ). To discover informative representations in a lower-dimensional space, each vector xt is passed through the decoder layers: ŷt = q(g(yt), yielding a reconstructed matrix ŶB . By minimizing the reconstruction error ‖ŷt − yt‖22, the model is expected to capture feature dependencies across time-series and encode this global information into a few latent variables in xt. Simultaneously, by maximizing the log-likelihood of (latent) future observations given compressed past input, the model targets to capture latent series progressions via time-conditioned Normalizing Flows. The proposed model is meaningful as it encourages the latent variables to capture different complex patterns of the data, which makes the representation more powerful and universal.
The key insight is that the encoding/decoding components and the time conditioned latent temporal generative model are jointly trained by minimizing a combined loss function over a given batch B of time-series consisting of a reconstruction loss-related term and a sequence negative log-likelihood term using stochastic gradient descent-based optimization. The overall loss term is
L(φ, φ′, θ) := 1 |B|(T + τ) ∑ y1:T+τ∈B T∑ t=1 ( ||yt − q(g(yt;φ);φ′)|| 2 2︸ ︷︷ ︸
LREC
−λ log p(q(g(yt;φ);φ′)|ht; θ)︸ ︷︷ ︸ LNegLL
) ,
(1) By minimizing a combined error of two tractable losses LREC and LNegLL, the model is given the capability to predict the future from latent representations of the observed history that preserve only the essential information, which is transferred to the decoder by minimizing the reconstruction loss. At the same time, it serves to reduce the dimensions of the temporal generative learning space. The
overall optimization problem is minφ,φ′,θ L(φ, φ′, θ) = 1|B|(T+τ) ∑ y1:T+τ∈B ∑T t=1 LREC(φ, φ′) + λLNegLL(θ), (2)
where λ ≥ 0 is a hyperparameter that balances the two losses. Given a trained model, our goal is to produce a set of future trajectories given the past (see Figure 2). Note that, instead of producing samples directly in a possibly very high dimensional feature space, the generator first produces Monte Carlo samples from the predictive distribution p(xT+1, . . .xT+τ |x1, . . . ,xT ) in the lower dimensional embedding space, by sequentially sampling from p(xt|ht), updating ht, and passing the latent sequences through qφ′ to map xt back to the original domain.

heading: 4 EXPERIMENTAL RESULTS

heading: 4.1 RESULTS ON PUBLICLY AVAILABLE DATASETS
text: We use publicly available datasets to evaluate the accuracy of forecasts of future observations yt ∈ RN compared to to a wide range of existing baselines (described in the Appendix). The Solar dataset, with N = 137 and Ttotal = 4832, contains hourly solar power production records. The Electricity dataset with N = 370 and Ttotal = 32303, is a hourly time-series of the electricity consumption of 370 customers. The Traffic dataset, with N = 963 and Ttotal = 4049, contains hourly occupancy rates of car lanes. Wiki, with N = 2000 and Ttotal = 792, is a daily time-series daily page views of Wikipedia articles. Taxi, with N = 1214 and Ttotal = 1488, is a spatio-temporal traffic time-series of New York taxi rides taken every 30 minutes.
Batch sets are formed by randomly sampling context and adjoining prediction horizon interval sized windows from the complete time-series history of the training data. For the test set, we perform rolling prediction evaluation. Following (Salinas et al., 2019a), Solar, Electricity, and Traffic, accuracy is measured on 7 rolling time windows, for Traffic we use 5 time windows, and for Taxi 57 windows are used in order to cover the full test set. We perform rolling prediction evaluation:
Data set Solar Electricity Traffic Taxi Wikipedia
VRNN (Chung et al., 2015) 0.133±0.009 0.051±0.001 0.181±0.002 0.139±0.015 3.400±0.322 KVAE (Fraccaro et al., 2017) 0.344±0.002 0.051±0.025 − − 0.099±0.012
TCNF LSTM - RealNVP (Rasul et al., 2020) 0.365±0.002 0.027±0.001 0.026±0.001 0.182±0.013 0.097±0.024 TCNF LSTM - MAF (Rasul et al., 2020) 0.378±0.032 0.025±0.002 0.028±0.002 0.172±0.001 0.097±0.033
Vec-LSTM - ind-scaling (Salinas et al., 2019a) 0.391±0.017 0.025±0.001 0.087±0.041 0.506±0.005 0.133±0.002 Vec-LSTM -lowrank-Copula (Salinas et al., 2019a) 0.319±0.011 0.064±0.008 0.103±0.006 0.326±0.007 0.241±0.033
GP - Scaling (Salinas et al., 2019a) 0.368±0.012 0.022±0.000 0.079±0.000 0.183±0.395 1.483±1.034 GP - Copula (Salinas et al., 2019a) 0.337±0.024 0.024±0.002 0.078±0.002 0.208±0.183 0.086±0.004
24 time-points per window, last 7 windows for testing for traffic and electricity, and 14 per window with last 4 windows for wiki.
We evaluate forecasting performance starting on equally spaced time points after the last observed training point. Since the marginal continuous ranked probability score (CRPS) cannot assess whether dependencies across time-series are accurately captured, we report the the Continuous Ranked Probability Score (CRPS) (CRPS-Sum). This metric is obtained by first summing across the N (ground-truth data and forecasted) time-series, which yields a CDF estimate F̂sum(t) for each time point. The results are then averaged over the prediction horizon. We take 100 samples to estimate the empirical CDF in practice. Table 2 lists the mean CRPS-Sum values averaged over 10 independent runs with standard deviations and shows that the model sets the new state-of-the-art on most of the benchmark data sets.

heading: 4.2 RESULTS ON AH&MS DATASET
text: Biometrics collected from wearable devices, including heart rate measurements and activity levels (e.g., step counts, standing hours) throughout the day, are rich sources of information that can yield crucial insights into the health trajectory of an individual. Vital signals from wearables show feasibility for accurate prediction of clinical laboratory measurements (Dunn et al., 2021). These measurements also show promise to detect acute infections—cosinor models fit to diurnal heart rate variability (HRV) patterns can detect pre-symptomatic COVID-19 infection (Hirten et al., 2021), and anomaly detection algorithms using resting heart rate (RHR) and step counts can identify presymptomatic COVID-19 (Alavi et al., 2021). Biometrics also show sensitivity to detect common colds (H1N1 and rhinovirus) (Grzesiak et al., 2021). Beyond acute infections, wearable measurements can provide insight into cardiovascular health, as higher RHR is significantly associated with coronary artery disease, stroke, and sudden death (Zhang et al., 2016). All-cause and cardiovascu-
lar mortality is also indicated by lower HRV (Singh et al., 2018). HRV is also found to indicate dysregulation of the autonomic nervous system, and thus is implicated in the development of hypertension (Schroeder et al., 2003). Joint modeling of vital signs from wearable devices continues to yield powerful signals for early detection of disease and health decline, as well as indicate general well-being and fitness. The relatively low cost of wearables and longitudinal nature of measurements compared to in-clinic evaluations motivates continued research in this area.
To this target, we apply LatTe-Flows for multivariate signal forecasting: given the past 30-day measurements of a group of vital sensor signals as input, we wish to predict the future 20-day values of the multivariate input time-series. We evaluate the empircal performance of our method on the Apple Heart and Movement Study (AH&MS) dataset. The AH&MS study was sponsored by Apple and conducted in collaboration with the American Heart Association and Brigham and Women’s Hospital. The study was approved by Advarra IRB and data were collected in accordance with the IRB approved consent form.
AH&MS dataset description: The dataset contains signals collected in real-world environments, which constitute a measure of cardiovascular and autonomic nervous system activity as well as movement and activity metrics. For the participants enrolled in this study, this dataset consists of passively collected sensor signal data from wearable devices in addition to self-reported measurements. This rich collection of signals comprises a partial view of individual biometric information as it evolves over time (see Figure 3).
Our data consists of passively collected observations: Resting Heart Rate (RHR), Heart Rate Variability (HRV), Step Count (SC), and user logged information such as Diastolic Blood Pressure (DBP), Systolic Blood Pressure (SBP), and Basal Body Temperature (BBT). HRV is calculated in the time-domain as the standard deviation between heartbeat measurements (defined by measured R-R intervals). Each time-series measurement is aggregated at the daily level and averaged if there are multiple measurements per day (for RHR, HRV, DBP, SBP, BBT) or summed (for SC).
AH&MS dataset analysis and results: Given a training set comprised of multiple context windows of 30 daily observations of BBT, DBP, SBP, HRV, RHR, and SC measurements of 1630 individuals,
the goal is to jointly learn each signal output for the next 20 days, by modeling the temporal dependencies of latent representations of the signals while also learning a meaningful lower-dimensional space. We project another signal captured in our dataset, V O2max, (omitted from training) onto the learned latent space to show that our model captures meaningful representations(Figure 5).
Using the Normalized Mean Square Error (NMSE) on each individual time-series, and CRPS-Sum as an evaluation metric, we qualitatively assess test-time signal predictions produced by our method for different latent-space dimensionality, D = {3, 4, 5}. We compare our approach against existing classical multivariate methods described in the previous section. Because of different scales in each signal, we first normalize each signal by the sum of its absolute values before computing this metric. The results are reported in Table 2, where in average, LatTe-Flows with D = 5 achieves the best signal predictions.
In Figure 4, we demonstrate the quality of predictions produced by the proposed model—we show the predicted median, 50% and 90% distribution intervals of for HRV, RHR, SC, BBT, DBP, and SBP in the future 20 day window over a randomly chosen individual. In this case, the reduced dimensionality is 4 (see more forecasting results in the Appendix). Interpretable representation learning on time-series is a fundamental problem for uncovering the latent structure in complex systems, such as sensor streams of vital signals. To reveal the underlying factors controlling accurate forecasts, we train the model with D = 2 and D = 3, and visualize the representations across 50 time points
for eight different individuals. These individuals belong to two disjoint cardiorespiratory fitness levels: the first group consists of four subjects with low values of estimated VO2max (Apple, 2021), VO2max < 25, and the second group consists of four subjects with high values of VO2max,VO2max > 40. Maximum oxygen consumption (VO2max) is the measurement of the maximum amount of oxygen a person can utilize during intense exercise and is considered the gold standard for determining an individual’s cardiorespiratory fitness and performance capacity. The latent representations for both D = 2 and D = 3 indicate that the mapping learned by the optimization leads to a space with easily visualized clusters that well-separates individuals according to VO2max (Figure 5). Our latent space estimated using RHR, HRV, BBT, DBP, SBP, and SC corroborates other recent methods that use step count, sedentary time, and moderate-vigorous physical activity to associate with VO2max (Nayor et al., 2021).

heading: 5 DISCUSSION
text: In this paper, we developed a general tool for multivariate forecasting that scales to large collections of dependent time-series. We have shown that the combination of representation learning, flexible density models and auto-regressive structures can produce accurate forecasts by modeling dynamics on a low-dimensional subspace. Our model is competitive with and often outperforms state-of-theart time-series forecasting methods on benchmark applications on publicly available datasets. We further analyzed a challenging practical application of predicting sensor signal trajectories over time. While here our focus is on time-series forecasting, future work will extend our framework to many tasks that require learning over high-dimensional time-series: imputation, interpolation, anomaly and out of distribution detection, and general sequence generation.

heading: 6 ETHICS STATEMENT
text: We study the joint trajectory of biometrics over time, a domain that has seen success across health applications. With this success and subsequent increase in research comes an increased need for model, algorithm, and data validation. We acknowledge that any algorithm deployed in a healthcare setting should be subject to validation efforts. We analyze data from the Apple Heart and Movement Study, a study containing information from human participants. The Apple Heart and Movement Study is a collaboration between Apple and Brigham and Women’s hospital to understand the factors underlying physical health. Participants in the study must sign an Advarra IRB approved informed consent form before being able to participate in the study. Participation in the study is entirely voluntary and without compensation.

heading: 7 REPRODUCIBILITY STATEMENT
text: In the Appendix we provide detailed information and description on the architecture and hyperparameters to reproduce the experiments for each publicly available dataset considered. Although a significant part of this paper considers the Apple Heart and Movement Study, which is not publicly accessible, we provide necessary insights and perform various experiments to showcase the applicability of the proposed model on practical healthcare-related tasks (see details in the Appendix).

Review summaries:
The authors propose a novel approach that scales well for multivariate sequence forecasting tasks. They offer to use conditional normalizing flow to capture the sequence dynamics in the latent space and then use auto-regressive architecture to decode the signal in the original space. The authors clearly explain why they prefer latent forecasting and experimentally verify the benefits of the proposed method on both widely used public datasets and the AH&MS dataset.