Efficient learning of nonlinear prediction models with time-series privileged information

In domains where sample sizes are limited, efficient learning algorithms are critical. Learning using privileged information (LuPI) offers increased sample efficiency by allowing prediction models access to auxiliary information at training time which is unavailable when the models are used. In recent work, it was shown that for prediction in linear-Gaussian dynamical systems, a LuPI learner with access to intermediate time series data is never worse and often better in expectation than any unbiased classical learner. We provide new insights into this analysis and generalize it to nonlinear prediction tasks in latent dynamical systems, extending theoretical guarantees to the case where the map connecting latent variables and observations is known up to a linear transform. In addition, we propose algorithms based on random features and representation learning for the case when this map is unknown. A suite of empirical results confirm theoretical findings and show the potential of using privileged time-series information in nonlinear prediction.


Introduction
In data-poor domains, making the best use of all available information is central to efficient and effective machine learning. Despite this, supervised learning is often applied in such a way that informative data is ignored. A good example of this is learning to predict the condition of a patient at a set follow-up time based on information of the first medical examination. Classical supervised learning makes use only of the initial data to predict the disease status at the follow-up, although in many cases data about medications, laboratory tests or vital signs are routinely collected about patients also at intermediate time points. This information is privileged (Vapnik and Vashist, 2009), as it is unavailable at the time of prediction, but can be used for training a model.
Learning using privileged information (LuPI) (Vapnik and Vashist, 2009), generalized distillation (Lopez-Paz et al., 2016) and multi-view learning (Rosenberg and Bartlett, 2007) have been proposed to increase the sample efficiency by leveraging privileged information in learning. Theoretical results guarantee improved learning rates (Pechyony and Vapnik, 2010;Vapnik et al., 2015) or tighter generalization bounds (Wang, 2019) for large sample sizes under appropriate assumptions. However, privileged information is not always beneficial; it must be related to the task at hand (Jonschkowski et al., 2015). Previous works do not identify such settings. Moreover, existing analyses do not state when learning with privileged information is preferable to classical learning for problems with small sample sizes-which is where efficiency is needed the most. Karlsson et al. (2022) studied LuPI in the context of predicting an outcome observed at the end of a time series based on variables collected at the first time step. They showed that making use of data from intermediate time steps in particular settings always leads to lower or equal prediction risk-for any sample size-compared to the best unbiased model which does not make use of this privileged information. However, their method called learning using privileged time series (LuPTS) was limited to settings where the outcome function, and estimators of it, are linear functions of baseline features. Figure 1: Comparison between the linear data generating process and the latent nonlinear generalization in this work. Ψ indicates the observation function of the latent system, with Φ its left inverse.
Moreover, their analysis did not study how variance reduction behaves as a function of increased input dimension. Hayashi et al. (2019) also learned from privileged intermediate time points but their study was limited to empirical results for classification using generalized distillation.
There is an abundance of real-world prediction tasks with fixed follow-up times which can be framed as having access to privileged time-series information. Examples include predicting 90-day patient mortality (Karhade et al., 2019) or patient readmission in healthcare (Mortazavi et al., 2016), the churn of users of an online service over a fixed period (Huang et al., 2012) or yearly crop yields from satellite imagery of farms (You et al., 2017). In these cases, privileged time-series information comprises daily patient vitals, intermediate user interactions and daily satellite imagery respectively. We are motivated by finding sample-efficient learning algorithms that utilize privileged time-series information to improve upon classical learning alternatives in such settings.
Contributions. We extend the LuPTS framework to nonlinear models and prediction tasks in latent dynamical systems (Section 3). In this setting, we prove that learning with privileged information leads to lower risk when the nonlinear map connecting latent variables and observations is known up to a linear transform (Section 3.1). In doing so, we also find that when the representation dimension grows larger than the number of samples, the benefit of privileged information vanishes. We show that a privileged learner using random feature maps can learn optimal models consistently, even when the relationship between latent and observed variables is unknown, and give a practical algorithm based on this idea (Section 3.2). However, random feature methods may suffer from bias in small samples. As a remedy, we propose several representation learning algorithms aimed at trading off bias and variance (Section 3.3). In experiments, we find that privileged time-series learners with either random features or representation learning reduce variance and improve latent state recovery in small-sample settings (Section 4) on both synthetic and real-world regression tasks.

Prediction and privileged information in nonlinear time series
We aim to predict outcomes Y in Y ⊆ R q based on covariates X 1 in X ⊆ R k . X 1 are observed at a baseline time point t = 1, starting a discrete time series on the form X 1 , X 2 , . . . X T , Y . Outcomes are assumed to be compositions of a representation Φ, a linear map θ and Gaussian noise , In addition to observations of X 1 and Y , privileged information (PI) is available in the form of samples of random variables, X 2 , ..., X T , from intermediate time points between X 1 and Y , all taking values in X . Two data generating processes (DGPs) with this structure are illustrated in Figure 1. Unlike baseline variables X 1 , privileged information is observed only at training time, not at test time. Therefore, it can only benefit learning, and not inference. Data sets of m ∈ N + training samples, D := {(x i,1 , x i,2 , ..., x i,T , y i )} m i=1 , are drawn independently and identically distributed from a fixed, unknown distribution p over all random variables in our system. We let X t = [x 1,t , ..., x m,t ] denote the data matrix of features observed at time t = 1, ..., T and Y = [y 1 , ..., y m ] the vector of all outcomes observed in D. A learning algorithm A : D → H maps data sets D ∈ D to hypothesesĥ ∈ H. An efficient algorithm A has small expected risk R p (A ) with respect to a loss L : Y × Y → R over training sets of fixed size m drawn from p, We study the regression setting with L the squared error, L(y, y ) = y − y 2 2 . In analysis, we focus on univariate outcomes, q = 1, but all results extend to multivariate outcomes, q > 1.

Algorithm 1: Generalized LuPTS
Input: Data D = ({X t }, Y); Repr.Φ or kernel κ if using a fixed representationΦ then Z t = [Φ(x 1,1 ), ...,Φ(x m,1 )] for t = 1, ..., T θ P := else if using kernel κ then K t = κ(X t , X t ) for t = 1, ..., T Our main goal is to compare privileged learners A P , which make use of privileged information, to classical learners A C which learn without PI. We seek to identify conditions and algorithms for which privileged learning is more efficient, i.e., it leads to lower risk in expectation, R(A P ) ≤ R(A C ) for the same number of training samples. We describe such a setting next.

Privileged information in latent dynamical systems
We study tasks where data are produced by a latent dynamical process. Observations X t are generated from unobserved latent states Z t through a nonlinear transformation Ψ, see Figure 1b. This means that the target of learning, h(X 1 ) = E[Y | X 1 ], is nonlinear in X 1 . Latent dynamical systems like these have proven successful at modelling a variety of phenomena as for example fluid dynamics in physics (Lee and Carlberg, 2019) and human brain activity in neuroscience (Kao et al., 2015). Assumption 1 (Latent linear-Gaussian system). Let Z 1 , ..., Z T be latent states related as a linear-Gaussian dynamical system in a space Z ⊆ R d with Z 1 of arbitrary distribution. Further, let the observation function be an injective map Ψ : Z → X with left-inverse (representation) Φ : X → Z. With A 1 , ..., A T −1 ∈ R d×d , β ∈ R d assume that Z 1 , ..., Z T , X 1 , ..., X T , Y are generated as where t ∼ N (0, σ 2 t I) and Y ∼ N (0, σ 2 Y ). It is easy to verify that Assumption 1 is consistent with (1), but stronger: there are more systems that map X 1 to Y as θ Φ(X 1 ) than systems where additionally X t = Ψ(Z t ). Nevertheless, it is much more general than the results of Karlsson et al. (2022), limited to the linear setting, i.e. Φ(x 1 ) = x 1 .
Next, we present a generalized version of the LuPTS algorithm of Karlsson et al. (2022) that is provably preferable to classical learning when data is generated according to Assumption 1 and Φ is known up to a linear transform. Further, we discuss how a privileged learner can be made universally consistent for unknown representations Φ when combined with random features.

Efficient learning with time-series privileged information
We analyze and compare learners from the privileged and classical paradigms which produce estimates of the formĥ(x 1 ) =θ Φ (x 1 ), as motivated by (1). We let each use representationsΦ : X →Ẑ ⊆ Rd from a family F shared by both paradigms, so that the hypothesis class H ĥ is shared as well.

Learning with true representations known up to a linear transform
In this section, we assume that privileged and classical learners have access to a common, fixed representation functionΦ which is related to the true representation function through a linear transform, meaning there exists a matrix B such thatΦ(·) = BΦ(·).
Classical learning. As a classical learner to serve as a strong comparison point for our privileged learners, we use the ordinary least-squares (OLS) linear estimator applied to the latent variables at the first time point inferred byΦ. WithẐ 1 = [Φ(x 1,1 ), ...,Φ(x m,1 )] ∈ R m×d , When m ≥d and Φ is known up to linear transform,ĥ C is the minimum-risk unbiased estimator of h which does not use privileged information. To accommodate the underdetermined case where m <d, the matrix inverse is replaced by the Moore-Penrose pseudo inverse (·) † (Penrose, 1956).
Putting this together, we arrive at Algorithm 1 which we call generalized LuPTS. We may apply a simple matrix identity and replace termsẐ tẐ t in Algorithm 1 by the Gram matricesK t = κ(X t , X t ) of a reproducing kernel κ with corresponding (implicit) feature mapΦ, κ(x, x ) = Φ (x),Φ(x ) (Smola and Schölkopf, 1998). This variant allows for learning with unknown Φ.
We now state a result which says that under Assumption 1, for an appropriate fixed representationΦ or kernel κ, generalized LuPTS is never worse in expectation than the classical learner in (3). Theorem 1. Let D be a data set of size m drawn from p, consistent with Assumption 1. Assume that the left inverse Φ : X → Z of the observation function Ψ is known up to a linear transform, explicitly or through a kernel κ(x, x ) = Φ (x),Φ(x ) , i.e., there exists a matrix B with linearly independent columns such thatΦ(x) = BΦ(x) withΦ : X →Ẑ for all x ∈ X . Then, it holds for the privileged learner A P (D) =ĥ P (Algorithm 1) and the classical learner Proof sketch. First, we show that the predictions made by generalized LuPTS are invariant to a linear transform B applied toΦ(·) during training, see Appendix A. Then we consider two cases: (i) When d ≤ m we may re-purpose the proof of Theorem 1 in Karlsson et al. (2022), (ii) When m <d Proposition 1 below directly implies Eĥ P ,X1 [Var D (ĥ C (X 1 ) |ĥ P )] = 0 and R(A P ) = R(A C ).
Theorem 1 implies that the privileged learner is at least as sample efficient as the classical one since Var(·) ≥ 0 and thus R(A P ) ≤ R(A C ) for the same number of training samples m. The result is a direct generalization of the main result in Karlsson et al. (2022). We also observe that generalized LuPTS and the classical learner coincide under certain conditions when m ≤d.
Proposition 1. LetΦ : X →Ẑ ⊆ Rd be any map with corresponding kernel κ. LetK t be the Gram matrix of κ applied to X t and letĥ C ,ĥ P be classical (3) and privileged (Algorithm 1) estimates. Then, K t is noninvertible whenever m >d, assuming linearly independent features.
Proof sketch. When the Gram matrices K t are invertible, the pseudo-inverse coincides with the inverse, and factors K t K † t in the LuPTS estimators cancel, making LuPTS and CL equal. Remarks. Theorem 1 extends the applicability of LuPTS to a) nonlinear prediction through a fixed feature map or b) kernel estimation and c) to the underdetermined case of m <d. While previous work is restricted to observed linear systems our result considers the case of latent linear dynamics which only need to be identified up to a linear transform (see Figure 1). Proposition 1 does not claim that no preferable privileged learner exists; it is a statement only about generalized LuPTS. We may relate the result to the double descent characteristic previously observed for other linear estimators for fixed m and varyingd (Loog et al., 2020). After a phase transition aroundd = m LuPTS's variance reduces for a second time when it becomes equivalent to the classical learner (see Figure 2).

Random feature maps for unknown representations
When the true Φ is entirely unknown, as is often the case in practice, using a poor representation Φ may yield biased results for both classical and privileged learners. A common solution in nonlinear prediction is to use a universal kernel, such as the Gaussian-RBF kernel. These have dense reproducing-kernel Hilbert spaces which allow approximation of any continuous function. However, universal kernels also have positive-definite and thus invertible Gram matrices (Hofmann et al., 2008), which according to Proposition 1 eliminates any gain in sample efficiency of generalized LuPTS.
Instead, we combine our algorithm with an approximation of universal kernels-random feature maps. These methods project inputs x ontod features by a random linear map W ∈ R k×d , and a nonlinear element-wise activation function. By choosingd < m, we can benefit from the function approximation properties of universal kernels (see discussion below) and the variance reduction of generalized LuPTS. Popular random features include random Fourier features (RFF) (Rahimi and Recht, 2007) and random ReLU features (RRF) (Sun et al., 2018) is the rectifier (ReLU) function and γ > 0 is a bandwidth hyperparameter.
For large enough numbers of random features and training samples, any continuous function h can be approximated up to arbitrary precision by a linear map ω applied to the random features, e.g.,  (2017). Applying the same argument to the step-wise estimators of LuPTS we can justify using random features in Algorithm 1 by the following observation: Under appropriate assumptions, we can construct a privileged learner using random feature maps which is a universally consistent estimator of h(x 1 ) = E[Y |x 1 ]. The precise construction deviates somewhat from Algorithm 1, but follows the same structure. We refer to Appendix B for a precise statement. As universal consistency describes the asymptotic behaviour in the limit of infinite samples and random features, this offers only limited insight into the benefits of privileged information in small sample settings, where performance will be a bias-variance trade-off.
Variance reduction & bias amplification. Generalized LuPTS is only guaranteed lower variance compared to the classical estimator under Assumption 1, although our empirical results (Section 4) suggest this applies more widely. WhenΦ is a bad approximation of Φ, generalized LuPTS may amplify bias, increasing with the number of privileged time points, compared to classical learning. We show this theoretically in Appendix C and also empirically in Appendix F. Whether generalized LuPTS is still preferable to classical learning in terms of prediction risk appears to depend on the amount of bias that gets amplified. Our experiments imply the variance reduction mostly dominates when using random features, whereas this is not always the case for linear LuPTS. The phenomenon of bias amplification is familiar from e.g., model-based and model-free reinforcement learning (Kober et al., 2013). As the bias with random features may still be high for small sample settings, we next present privileged representation learning algorithms to trade off bias and variance more efficiently.

Privileged time-series representation learning
Up until now the representationΦ was considered fixed, either because Φ was known up to a linear transform (explicitly or implicitly) or because of the use of random feature methods. Generalized LuPTS (Algorithm 1) produces minimizers {Â t },β of the following objective for fixedΦ, Objective (5) and the systems described by Assumption 1 lend themselves to methods which also learn the representation Φ in addition to the latent dynamics {A t , β}. Next, we present three algorithms which combine the ideas of generalized LuPTS with the expressiveness of deep representation learning. All learners use equivalent encoders to representΦ(·) and linear layers to model the relations between the latent variables {Z t } and the outcome Y . The classical learner predicts the outcome linearly fromẐ 1 . All architectures under consideration are visualized jointly in Figure 3a.
Learner & Φ (a) Classical (left) and privileged (right) representation learners.Φ is shared across all time steps. GRL models the direct mapsθt to the outcome; SRL models the single stepsÂt andβ. CRL combines the two.  SRL. The first privileged representation learner directly optimizes objective (5), just like generalized LuPTS, but now also fitting the representationΦ, parameterized by a neural network. We refer to this model as stepwise representation learner (SRL). As we will see in experiments, a drawback of this approach is that representations may favor predicting transitionsẑ i,t →ẑ i,t+1 with small error, while losing information relevant for the target outcome in the process. At test time, for a new input x 1 , SRL composes the stepwise dynamics to outputĥ P (x 1 ) =β Â T . . .Â 1Φ (x 1 ).

CRL and GRL.
To make sure that the learned representationΦ retains information about Y , we add linear outcome supervision to the representationΦ(X t ) at each time step t. Recall that, by Assumption 1, the expected outcome is linear in the latent state at any time step. We introduce a hyperparameter λ ∈ [0, 1] to trade off the two types of losses and arrive at the combined representation learner (CRL). With α = (Φ, {Â t }, {θ t }) the entire parameter vector, CRL minimizes the objective We make test-time predictions usingĥ P (x 1 ) =θ 1Φ (x 1 ). In experiments, we highlight the case λ = 1 where only outcome supervision is used as greedy representation learner (GRL). For a precise definition of the GRL objective, see Appendix D. GRL is related to multi-view learning, in which prediction of the same quantity is made from multiple "views" (cf. time points) (Zhao et al., 2017).

Experiments
We compare classical learning to variants of generalized LuPTS (Algorithm 1) and the privileged representation learners of Section 3.3 on two synthetic and two real-world time-series data sets. We (i) verify our theoretical findings by analyzing the sample efficiency and bias-variance characteristics of the given algorithms; (ii) demonstrate that generalized LuPTS with random features succeeds in settings where linear LuPTS suffers from large bias; (iii) point out that privileged representation learners offer even greater sample efficiency in practice and (iv) study how well these algorithms recover the true latent variables {Z t } and how this relates to predictive accuracy.
Experimental setup. We report the mean coefficient of determination (R 2 ), proportional to the squared-error risk R, for varying sample sizes, sequence lengths and prediction horizons. Experiments are repeated and averaged over different random seeds. In each repetition, a given model performs hyperparameter tuning on the training data using random search and five-fold cross-validation before being re-trained on all training data. The test set size is 1000 samples for synthetic data and 20% of all data available for real-world data sets. We consider six privileged learners of two groups. The first group comprises generalized LuPTS with the linear kernel (LuPTS) and the two random feature maps shown in Section 3.2: Random Fourier features (Fourier RF) and random ReLU features (ReLU RF). The classical learners for this group are OLS estimators used with the same kernel or feature map. The second group consists of the representation learners SRL, CRL and GRL. For tabular data, their encoder is a multi-layer perceptron with three hidden layers of 25 neurons each. For the image data they use LeNet-5 (LeCun et al., 1989). The classical learner (Classic Rep.) uses the same encoder with a linear output layer. The results presented were found to be robust to small changes in training  parameters such as learning rate. For details on the training process we refer to Appendix E. All experiments required less than 3000 GPU-h to complete using NVIDIA Tesla T4 GPUs. 2 Data sets. We briefly describe evaluation data sets and refer to Appendix E for further details. First, we create two synthetic data sets in which latent states and outcomes are generated from linear-Gaussian systems as in Assumption 1. To produce the observations {X t } we use a deterministic nonlinear function Ψ : Z − → X . In the first synthetic data set, which we call Square-Sign, the nonlinear transformation Ψ : . The second synthetic data set uses the same latent linear system with Z t ∈ R 2 and produces square images (28×28 pixels) as observations. As the images are reminiscent of clocks we refer to this data set as Clocks-LGS. Example sequences of these observations are presented in Figure 3b. The angle, size and fill of the two clock hands encode the value of the corresponding latent variable. The outcome Y is a linear function of Z T with q = 1. For Square-Sign, q = 3.
The Metro Interstate traffic volume data set (Traffic) (Hogue, 2012) contains hourly records of the traffic volume on the interstate 94 between Minneapolis and St. Paul, MN. In addition, the data contains weather features and a holiday indication. We predict the traffic volume for a fixed time horizon given the present observations. Privileged information is observed every four hours.
We also predict the progression of Alzheimer's disease (AD) as measured by the outcome of the Mini Mental State Examination (MMSE) (Galea and Woodward, 2005). The anonymized data were obtained through the Alzheimer's Disease Neuroimaging Initiative (ADNI) (ADNI, 2004). The outcome of interest is the MMSE score 48 months after the first examination. Privileged information are the measurements at 12, 24 and 36 months. In addition we tested our algorithms on the PM 2.5 data set (Liang et al., 2016), where we predict the air quality in five Chinese cities (see Appendix F).

Sample efficiency, bias and variance
The main goal of our work is to improve learning efficiency by incorporating privileged time-series information. Across almost all prediction tasks and sample sizes, nonlinear variants of generalized LuPTS outperform their classical counterpart in terms of sample efficiency as can be seen in Figure 4. On the Traffic prediction task, LuPTS ReLU RF outperforms linear LuPTS as the former appears to exhibit less bias. On the synthetic data of Square-Sign, linear models reach their best accuracy quickly, while they are limited by their lack of expressiveness. Generalized LuPTS amplifies this bias, making linear LuPTS worse than OLS. Random feature methods attain higher accuracy here but are generally less sample efficient. Generalized LuPTS combined with random features manages to decrease this gap significantly. On ADNI, we don't see a benefit of using nonlinear models in general. We make similar observations on the PM 2.5 air quality data set, see Appendix F.  The representation learners proposed in Section 3.3 are evaluated on the same data sets and on the image prediction task Clocks-LGS. Example results are displayed in Figure 5. These demonstrate that directly transferring the LuPTS objective to neural networks in the form of SRL results in subpar performance. SRL does not seem to have a strong enough incentive to learn representations which accurately predict the outcome. Generalized distillation for privileged time series as suggested by Hayashi et al. (2019), does not improve upon classical learning in our tasks. On all experiments displayed in Figure 5, CRL and GRL outperform the classical learner. The predictive accuracy of these models is similar on most tasks, as may be explained by the fact that CRL may reduce to GRL when choosing λ = 1 in objective 6. Noticeably, the general observation of GRL and CRL being more sample efficient than the classic model, neither appears to depend on the neural architecture used for the encoder, nor does the modality of the data play an important role, as the image prediction task Clocks-LGS (Figure 5c) demonstrates. For additional empirical results, we refer to Appendix F.
To analyze the bias and variance characteristics of our algorithms, we can estimate the expected squared prediction bias, DGPs in addition to the variance of the different estimators. Figure 6 depicts bias and variance for all models on the Square-Sign data. On the left panel, all variants of generalized LuPTS exhibit lower variance than classical learning, despite being biased. This holds generally: Across all experiments, we never encountered an example where the use of privileged information has not resulted in lower variance compared to classical learning. On the contrary, the privileged learners suffer higher bias than the comparable classical learners because of the bias compounding over the individual prediction steps as shown in Appendix C. For the random feature variants however, the bias decreases with the number of samples. The representation learners display similar characteristics on the right panel of Figure 6. Learning the transitions between latent variables Z t appears to be associated with low variance and high bias as demonstrated by the results of SRL. GRL however, which does not model these transitions, exhibits the lowest bias and the largest variance of all privileged learners. As CRL is able to trade off between these two objectives, the estimates for its variance and bias lie in between the corresponding values of the other two privileged learners.   Figure 5c).

Latent variable recovery
Under Assumption 1, it is sufficient to identify the representation function Φ up to a linear transform B to have provable gains from privileged information over a classical learner. In synthetic data, we can assess to what extent a representationΦ with this property has been found. As a proxy for the existence of such a transform, we can compute a measure of correlation betweenΦ(X) and Φ(X), such as the Canonical Correlation Analysis (CCA) (Hotelling, 1936). Raghu et al. (2017) introduced a modified version called Singular Vector Canonical Correlation Analysis (SVCCA) to compute correlations when dealing with noisy dimensions in neural network representations.
The mean SVCCA coefficientsρ and predictive accuracy of the representation learners are visualized in Figure 7a. One notices that GRL and CRL produce higher correlation coefficients than the classical learner while also predicting the outcome on the Clocks-LGS task more accurately. Further comparing the representations learned by privileged and classical models, Figure 7b shows a visual example of the improved latent recovery of CRL on the same task. Concluding, the use of privileged time-series information does not only increase the sample efficiency of existing algorithms but can also aid the recovery of latent variables in latent dynamical systems.

Related work
Existing analyses for learning with privileged information guarantee improvements in asymptotic sample efficiency under strong assumptions (Vapnik and Vashist, 2009) but are insufficient to establish a clear preference for LuPI learners for a fixed sample size. For example, Pechyony and Vapnik (2010) showed that for a specialized problem construction, empirical risk minimization (ERM) using privileged information can achieve fast learning rates, O(1/n), while classical (non-privileged) ERM can only achieve slow rates, O(1/ √ n). We are not aware of any generalization theory tight enough to establish a lower bound on the risk of a classical learner larger than an upper bound for a PI learner. Our problem is also related to multi-task (representation) learning (Maurer et al., 2016), see especially (5)-(6). However, our goal is different in that only a single task is of interest after learning.
In estimation of causal effects, learning from surrogate outcomes (Prentice, 1989) has been proposed as a way to increase sample efficiency. Surrogates are variables related to the outcome which may be available even when the outcome is not (Athey et al., 2019). We can view these as privileged information. While the problem shares structure with ours, the goal is to compensate for missing outcomes and analytical results give no guarantees for improved efficiency when both surrogates outcomes are always observed (Chen et al., 2008;Kallus and Mao, 2020). Guo and Perković (2022) showed that, in the context of a linear-Gaussian system on a directed acyclic graph, a recursive least-squares estimator is the asymptotically most efficient estimator of causal effects using only the sample covariance. In the linear case on a path graph, their estimator coincides with ours. However, no analysis is provided for the nonlinear case or for fixed sample sizes.

Discussion
We have presented learning algorithms for predicting nonlinear outcomes by utilizing time-series privileged information during training. We prove that our estimator is preferable to classical learning when data is generated from a latent-variable dynamical system with partially known components. The proof holds for the case when the latent dynamical system is recovered by the representation function up to a linear transform, assuming that a left inverse of the true observation generating function exists. However, this assumption does not appear to be necessary, as our empirical results demonstrate that privileged learning is preferable to classical learning even when these conditions cannot be guaranteed. Consequently, a more general theoretical result where less is known about the latent system seems attainable. For example, one might consider a case in which only a few independent components of the latent variables are recovered by the representation used by the learning algorithm, while other components are treated as noise.
When the latent dynamical system is entirely unknown, we create practical estimators using random feature embeddings which outperform the corresponding classical learner across experiments. We show that a universally consistent learner can be constructed based on this idea, with slightly different form. As a further alternative, we propose representation learning methods of related form using neural networks and demonstrate the empirical benefits also of this estimator over classical learning. In experiments, we analyze how the gap in risk between privileged and classical learning changes for different prediction horizons, as displayed in Figure 16 in Appendix F. The results suggest that the risk advantage of privileged learners grows with the sequence length of the prediction task, despite the fact that the task becomes more difficult at the same time.
Our work focuses on the setting where data is observed as a time series. This setting is chosen for its simple causal structure given by (latent) linear-Gaussian systems and because it can be motivated from many different applications. However, the ideas presented are not specifically tied to time and also apply in the case when all variables are observed simultaneously as long as the causal structure remains sequential. Moreover, we believe that the theory presented here is not limited to sequential settings and generalizes to other causal structures, in particular directed acyclic graphs that connect the baseline covariates to the outcome. In either case, one might only have access to the baseline variable at test time. For example, the timely collection of data for all covariates at test time might be very expensive or even impossible.
As pointed out before, Theorem 1 requires the recovery of the latent variables up to a linear transformation. Nonlinear independent component analysis (ICA) (Hyvarinen and Morioka, 2016) aims to solve precisely this problem and has been applied to time series via time-contrastive learning. This makes for an interesting connection between learning using privileged information and nonlinear ICA, as the experiments of Section 4.2 suggest that privileged time-series information aids the recovery of latent variables. Other remaining challenges include providing risk guarantees for learning with biased representations (including deep neural networks), with regularized estimators, and for more general data generating processes with weaker structural assumptions. We are hopeful that the utility demonstrated in this work will inspire future research to overcome these limitations.  Only original code. (d) Did you discuss whether and how consent was obtained from people whose data you're using/curating? [N/A] (e) Did you discuss whether the data you are using/curating contains personally identifiable information or offensive content? [Yes] ADNI data is anonymized. 5. If you used crowdsourcing or conducted research with human subjects... Our proof requires an additional technical assumption: that the matrix of true latent states Z t = [Φ(x 1,t ), ..., Φ(x m,t )] , for all t, for a random data set D has independent columns with probability 1. This implies that rank(Z t ) = d and m ≥ d. We consider this a minor restriction since it would only be violated if either a) two or more components of Z t were perfectly correlated-in this case, a smaller system with same distributions over observations could always be constructed-or b) if we observe fewer samples than necesary to determine the system (m < d). Note that this does not require that the dimensiond of the estimated representationΦ is smaller than m. We begin by proving that both classical and LuPTS estimators are invariant to a particular form of linear transformation of the representationΦ.
Lemma 1. Assume we have a latent linear Gaussian system as defined in Assumption 1 such that for a data set D of m samples, the matrix of true latent states Z t = [Φ(X 1,t ); ...; Φ(X m,t )] has linearly independent columns with probability 1. Let A Φ P be the LuPTS algorithm using the system's true map Φ(·) with Φ(Ψ(z)) = z ∀z ∈ Z. Let AΦ P be the same algorithm using a different mapΦ(·). We assume that ∃B :Φ(x) = BΦ(x)∀x ∈ X . Analogously, we denote the classical learners AΦ C and A Φ C . If B ∈ Rd ×d has linearly independent columns we havê Proof. Let Z t ∈ R m×d be made up of the rows Z t(i,:) = Φ(X t(i,:) ) when X t ∈ R m×k is the design matrix belonging to data set D. In the same fashion we defineẐ t ∈ R m×d using the mapΦ instead. By assumption, Z t and B ∈ Rd ×d have independent columns such that B † B = I. These assumptions are used for matrix identities involving the Moore-Penrose inverse below. We compute the prediction on a new test point x for the classical learner: The arguments for the privileged learners are analogous: Proof of Theorem 1. We consider the generalized LuPTS estimatorĥΦ P (·) treating different cases for the number of samples m and latent state dimension d in turn. By the added technical assumption, that the true latent state Z t = [Φ(x 1,t ), ..., Φ(x m,t )] has rank d for all t with probability 1, and the assumption thatΦ is linearly close to Φ, by a matrix B ∈ Rd ×d of rank d such thatẐ t = Z t B , we get that rank(Ẑ t ) = d for all t. This also implies d ≤d.
(i) m = d: In this case, the Gram matrix K t =Ẑ tẐ t has full rank and thus is invertible for all t. By Proposition 1, LuPTS coincides with the classical learner. Hence, Eĥ P ,X1 [Var D (ĥ C (X 1 ) |ĥ P )] = 0 and R(A P ) = R(A C ).
(ii) m < d: In this case, there does not exists a linearly close, as defined above, representationΦ to Φ since the rank ofẐ t must be smaller than d. This contradicts that rank(Ẑ t ) = d. Independently, if the conditions of Proposition 1 hold, the same equivalence holds as in the case m = d.
(iii) m > d: In this case, the kernel Gram matrix K t =Ẑ tẐ t has rank d < m and is never invertible.
Three sub-cases remain: a) Whend = d, the matrix B is invertible and square, the covariance matrix Σ t =Ẑ tẐt is invertible for all t and our estimator coincides with linear LuPTS (Karlsson et al., 2022) in the space implied byΦ. To see this, note that Lemma 1 implies thatĥΦ P (·) makes the same predictions as a different generalized LuPTS estimatorĥ Φ P (·) using the true map Φ when the two representation functions are related through B, as defined in the Theorem statement. Consequently, we may analyze the latter estimator instead of the first. It uses the parameter We know by assumption that the covariance matrices Σ t = (Z t Z t ) ∈ R d×d have full rank for all t. This implies that the Moore-Penrose pseudoinverse (·) † may be replaced by the regular matrix inverse (·) −1 in the expression above, yieldinĝ which is equivalent to the LuPTS estimator of Karlsson et al. (2022) used on a linear-Gaussian system in the space Z rather than in X . In this case, Theorem 1 from Karlsson et al. (2022) yields the desired result. b) Ifd > d,Σ is not invertible but, due to Lemma 1, we can instead study a representation which is an appropriate linear transform B ∈ Rd ×d away fromẐ, and apply the Karlsson et al.
(2022) result as described for the cased = d. Note that in this case B is non-square but has linearly independent columns as required. c) Ifd < d, the assumed matrix B cannot exist with the stated conditions (the assumptions of Theorem 1 are not satisfied).

B Universality of random features
A learning algorithm A is said to be universally consistent if, for any continuous function h, the output of A converges in probability to h. That is, for a random dataset D m of m i.i.d. samples drawn from a distribution p, and any > 0, lim Sun et al. (2019) prove that (norm-bounded) linear regression applied to random ReLU features (RRF) is universally consistent. Specifically, for any , δ > 0, there is a finite number of random featuresd and samples m, such that the estimator achieves an error of at most with probability ≥ 1 − δ for univariate continuous functions of x. Universal consistency (asd → ∞, m >d)) follows as a result. We can apply the same idea to a version of generalized LuPTS by considering each parameter estimate of the latent dynamical system givenΦ, with norms restricted by R. We drop the subscript RRF fromΦ moving forward, and continue to use random ReLU features. For the final prediction step of Y , we let Progressing recursively backward from t = T , we let whereΦ γ,dt are random features specific to t. Since the target of each prediction at t + 1 is fixed with respect to the features used as input at t, the result from Sun et al. (2019) can be used to give a learning guarantee for each step. The construction differs from the standard generalized LuPTS formulation asΦ is not shared between time steps, and soÂ t will be non-square in general. We believe that this is merely a limitation of the proof technique and that the results hold for shared random features and square transitions.
Let H ReLU denote the set of linear functions applied tod random ReLU features, with uniform random projection coefficients ω j ∼ U(S d ), j = 1, ...,d.
Corollary 1 (Follows from Proposition 5 in Sun et al. (2019)). Let Assumption 1 hold with noiseless transitions and outcomes, t = 0 for t = 2, ..., T , and Y = 0. Define g * t (x t ) := Ψ(A t Φ(x t )) = x t+1 and assume that for any fixed RRF representationΦ γ,dt+1 , each component j = 1, ...,d t+1 of the transition target satisfiesΦ γ,dt+1 (g * t (·))(j) ∈ H ReLU,dt . LetÂ t be the minimizer of the single-step transitions as defined in (8). Then, for any δ > 0, > 0 there is a number of random featuresd =d( , δ) and samples m = m( , δ), such that with probability  (7), Then, further assume that the largest eigenvalue λ max (Â t ) ≤ C for any t and β ≤ C. Then, with probability at least 1 − δ, Proof. Let · = · L 2 (p) . Then, lettingΦ t =Φ γ,dt , and applying a union bound to each of the T ( , δ)-assumptions, and a series of Cauchy-Schwarz inequalities, Remark 1. Proposition 2 shows that a privileged learner with random ReLU features can be turned into a universally consistent estimator of any (noiseless) continuous function of X 1 by letting each time step have its own random feature representation of appropriate size and adding a norm constraint to each linear transformation. The construction in (8), (7) deviates from Algorithm 1 primarily in that the random feature representations used at each time step are different, but the overall structure is maintained: At training, predictions of Y are made from an embeddingẐ T =Φ T (X T ) of X T , predictions ofẐ T are made from an embeddingẐ T −1 =Φ T −1 (X T −1 ), and so on. At test time, the transition matricesÂ 1 , ...,Â T ,β are multiplied and applied toẐ 1 =Φ 1 (X 1 ). We conjecture that a similar argument can be applied to the construction in Algorithm 1 by letting the dimension of each time step approach ∞.

C Compounding bias
We can describe the compounding bias of the LuPTS estimator due to a biased representationΦ, in comparison with the standard OLS estimator, by propagating the error inΦ through the estimates. Assume that Y = θ Φ(X 1 ) + = (A 1 · · · A T −1 β) Φ(X 1 ) + Then, let Z t = Φ(X t ) and for an estimateΦ, assumed for simplicity to have the same dimension, Fittingθ toΦ using the classical learner (OLS) yields an estimatê where ∆ 1 = −Σ −1 1 Ω 1 (Σ 1 + Ω 1 ) −1 ,θ * c is the OLS estimate of θ for the true Φ and the second line follows from the Woodbury matrix identity. The norm of ∆ 1 is related to the condition number of Σ 1 . The expectation of the first term is θ, and the expectation of the remaining terms is the bias. Now, we can do the same thing for the privileged estimator. Let's start with T = 2.

Thus, the difference in bias between the two estimators is
More generally, we can express this difference recursively as below.
Proposition 3. Letθ p :=Â 1 · · ·Â Tβ be a privileged estimator using a linearly biased representation Φ, and letÂ * t be the same estimator using an unbiased representation Φ * . Then, the bias ofθ p is where E t is the compounded error in transition dynamics, computed recursively as follows

D Privileged time series representation learners
We expand on the description of the greedy representation learner (GRL) described as a special case of CRL in Section 3.3. To avoid the information loss of SRL, we consider its conceptual opposite, using privileged time series information only to predict the outcome. To do this, a linear output layer θ tẐt is used to predict Y at every time step t. Recall that, by Assumption 1, the expected outcome is linear in the latent state at any time step. The method is related to multi-view learning, in which prediction of the same quantity is made from multiple "views" (Zhao et al., 2017). We dub the model greedy privileged representation learner (GRL), which minimizes the objective During inference this algorithm returnsĥ P (·) =θ 1Φ (·). Compared to the objective of CRL in 6, we introduce an additional hyperparameter λ ∈ (0, 1) to place more weight on the loss term that is relevant at inference time. As a consequence, we choose w t = λ for t = 1 and w t = 1 − λ otherwise. We expect GRL to have less bias than SRL, but higher variance since less structure is imposed on the representationΦ.

E Experiment setup & data processing E.1 Detailed experiment setup
In the following we give a detailed description of the experimental setup used to obtain the results presented in Section 4 as well as the additional results that are part of this section. For a given data set, we select a combination of training set sizes and sequence length. For each unique combination of these parameters the models of interest are trained repeatedly with different random sampling. For each repetition the data is split into a train and a test set randomly before hyperparameter tuning and model training are performed. At last each model's predictions on the test set are scored by computing the coefficient of determination R 2 . On synthetic data the test set contains 1000 samples.
In the case of real-world data, where samples are limited, we test on 20% of all available data. The preprocessing used for real-world data and the generation procedure of synthetic data is unique to each data set. We refer to the data set specific subsections for detailed descriptions of how each data set is processed. During the experiments each model uses standardized data for training and inference. To perform the data rescaling we use the StandardScaler implementation that is part of scikit-learn (Pedregosa et al., 2011).
Hyperparameter tuning. The tuning of hyperparameters is carried out for each repetition and is implemented using random search and five-fold cross-validation. Each hyperparameter is sampled from a fixed interval of possible values. An overview of the ranges of different hyperparameters determined through random-search is provided in Table 1. For the experiments with variants of generalized LuPTS we sample ten sets of hyperparameters before retraining on all training data using the best set of parameters. For the representation learners we merely sample five values for λ.
Neural network training. The training of neural networks involves many choices and hyperparameters. We choose PyTorch's implementation (Paszke et al., 2019) of the Adam optimizer (Kingma and Ba, 2017) to train the representation learning models. If not specified otherwise the results shown in this project are obtained using a learning rate of 0.0001, a batch size of 30, leaky ReLU activations and a maximum of 1500 training epochs. In the case of neural network models, the sample sizes reported as part of the experiments denote the combined size of the training and validation set, where the validation set contains 20% of those samples. We use early stopping during the training process by keeping track of the validation loss. If a model does not improve the validation loss over a waiting period of 200 epochs we stop training early and set the network parameters to the values that obtained the lowest validation loss up until that point. In order to make sure that results are not dependent on the specific choice of the parameters just described, we performed additional experiments with different parameter choices. The results were found to be robust to small changes in these parameters.
Generalized distillation. In order to compare our algorithms to the alternative of using generalized distillation for privileged time series as presented by Hayashi et al. (2019), we implemented a model that (i) produces hypotheses of the same class as our other algorithms and (ii) that adopts the learning paradigm of a student model incorporating soft targets produced by a teacher model into its loss function. For tabular data our teacher model is a multi-layer-perceptron (MLP) with T * k input neurons such that all {X t } are concatenated and then used as input. The teacher MLP makes use of five hidden layers, each consisting of 100 neurons. In the case of the image data generated by Clocks-LGS, the teacher uses an implementation of LeNet-5 on all variables X t (while sharing the encoder parameters) before concatenating the 25-dimensional output of this encoder from different time steps. This combined representation is then processed by an MLP with a single hidden layer with 25 neurons. The loss function used for the student model producing the estimateĥ C is architecturally identical to the classic representation learner used in each of the experiments. When training the student model the mean squared error on the data and the error corresponding to the soft targets of the teacher model are combined via a hyperparameter λ: The hyperparameter is determined via hyperparameter tuning in all repetitions as described for other hyperparameters in this section. All other training procedures follow the same logic as described above.
Resources. For the training of the representation learning algorithms we use a cluster of graphics processing units (GPUs) in order to reach the number of experiment repetitions required for our work.
A single experiment like shown in Figure 5c takes several hours on 100 NVIDIA Tesla T4 GPUs. While the random feature methods do not require GPU training, they still require hyperparameter tuning which is why we compute results such as presented in Figure 4 on many CPU cores in parallel.
While the experiments on neural networks cannot reasonably be reproduced on a single desktop machine, this is still possible within a few days for the random feature methods.
Latent variable recovery and SVCCA. In order to assess to what extent a representation has been found that is linearly related to the true latent variables we use SVCCA as described by Raghu et al. (2017), meaning we first use PCA retaining at least 99% of variation before then applying CCA. For the visualization shown in Figure 7b we construct a grid of points (150 × 150) around the origin and assign each point a unique color. Then we compute an image using the observation generating function Ψ : Z − → X of the Clocks-LGS data set for each point. The observations are passed to the encoderΦ of the two representation learners producing estimates of the latent variables. We then map the estimates to the ground truth linearly using SVCCA before plotting the result.

E.2 Alzheimer progression
To test our algorithms on the task of predicting the progression of Alzheimer's disease (AD) we use an anonymized data set obtained through the Alzheimer's Disease Neuroimaging Initiative (ADNI) (ADNI, 2004) under the LONI Research License. The initiative is large multi-site research study on the brains of over 2000 AD patients which collects many features such as genetic, imaging and biospecimen biomarkers. The data consists of measurements taken every 3 months with some observations missing. The outcome of interest in our experiments is the Mini Mental State Experiment (MMSE) score 48 months after the first measurement (Galea and Woodward, 2005). Privileged information are the measurements taken between those time points, at 12, 24 and 36 months into the program.
Data processing. The processing procedure used in this project is borrowed directly from the work of Karlsson et al. (2022). There is a large amount of missing information in the ADNI data set. The missingness varies with the time of when measurements were taken. Further some subjects were not present at some of the follow-up examinations. To deal with the missingness patients without an observation for the final follow up (the outcome Y ) are excluded from our experiments. Further, we also require that patients are present at all intermediate time steps (12,24 and 36 months after the first measurement) which we use as privileged information. We one-hot encode categorical features and exclude features for which more than 70% of the observations are missing. To deal with the remaining missing values, mean imputation is used. Due to the filtering that we apply as a result of the missing data we only obtain 502 suitable sequences that we can use for our experiments.

E.3 Traffic data
The Traffic data set (Hogue, 2012) obtained through the UCI machine learning repository (Dua and Graff, 2017) contains hourly measurements of the traffic volume as well as weather features and holiday information. The raw data contains 48.204 records. An overview of all available features is given in Table 2.
Data processing. We noticed extreme outliers in the data set as well as implausible numerical values for the temperature and rain features. Further, records for some of the hours of the timeframe (2012 -2018) covered by this data set are missing. To deal with the extreme outliers we calculate the mean and standard deviation of each feature and remove records which contain values that are further than six standard deviations from the mean of a particular feature. We also remove a feature entirely if there is no variation left after this filtering. This is the case for the snowfall feature as snow is very rare in Minneapolis. From the date and time of each record we calculate the weekday which we add as a one-hot encoded feature and also represent the hour of the day h ∈ {0, 1, . . . 23} as two separate periodic features given by This ensures that a timestamp just before midnight produces similar features compared to just after midnight. We one-hot encode the holiday information, making no difference between different types of holidays, and make this feature persist over a full calendar day. In the original data set the holiday information is only specified for the first hour of the day. The column Weather_Main contains some weather conditions that are very rare, such as smoke and squall. As a consequence we group the different conditions before encoding them as binary variables. In particular we make drizzle, rain and squall one single feature while also grouping together fog, haze, mist and smoke as they all affect visibility.
Time series selection. After this preprocessing, that leaves only numerical values and one-hot encoded categorical values, we group the data together as time series used for the experiments. In order to do so we specify a desired sequence length T + 1 and a sequence step size in hours. With this information we iterate through the data set assembling time series with i) no values missing ii) the correct length and step size and iii) at least a seven hour gap between each pair of sequences. The third condition is introduced to make sure one does not end up with very similar cases (for short sequences in particular) in training and test set.

E.4 Square-Sign
The Square-Sign data set serves as a test environment for learning from privileged time series information where one can assure the conditions necessary for Assumption 1 to hold. In particular this means creating a linear-Gaussian system which remains unobserved and combining it with an observation generating function ψ : Z − → X .
Latent linear-Gaussian system. The first component that makes up the generation process for Square-Sign (and Clocks-LGS) is the linear-Gaussian system which is latent, just as depicted on the right side of Figure 1b with Z t ∈ R d . The first step in the data creation process is sampling each of the d components of Z 0 from N (0, 5). Then the subsequent latent variables Z t+1 are computed as . For the outcome we use the same form but with different dimensionality: Off-diagonal elements of the transition matrices A t ∈ R d×d are sampled from a Normal distribution N (0, 0.2) while the diagonal elements are set to one. In a second step we compute the spectral radius of the randomly created matrices A t via eigenvalue decomposition, obtaining the components U ΛU . We then set the spectral radius to a predefined value λ max = 1.3 and reassemble the matrix as The coefficients of β are drawn from the same normal distribution as the ones of A t but undergo no further changes. Observation generating function. As the dimensionality d of the latent space Z is not fixed we use an observation generating function that is not restricted to a specific value of d. For each element in Z t ∈ Z = R d we create two elements in X t ∈ X = R 2d by denoting its sign separately from its square. This gives the following nonlinear observation generating function: ,1) ), . . . , Z 2 (t,d) , sgn(Z (t,d) )]

E.5 Clocks-LGS
This data set serves the purpose of testing our algorithms on a different modality with high dimensional data. In particular the idea was to use image data as this is a domain where neural networks have been very successful. For this reason we combine a latent dynamical system with an image generation process which we explain in detail in this section.
Latent linear-Gaussian system. We use exactly the same setup as we do for the Square-Sign latent dynamical system as described in Section E.4. The only difference here is the dimensionality of the latent variables, transition matrices and the outcome. For Clocks-LGS we generally have d = 2 and q ∈ {1, 2}.
Image generation. The second part of Clocks-LGS is creating images from two dimensional latent t ] . The goal was to keep it the process simple while using small black and white images of 28×28 pixels. In addition we wanted each image to have no ambiguity with respect to the latent state it represents. We represent the first component by a clock hand mounted at the image center. One can think of Z (1) t as the angle in radian, meaning the hand points straight up for Z (1) t = 0 or Z (1) t = 2π and straight down for Z (1) t = π. To visualize a full rotation we increase the size of the cirle around the image center in discrete steps for each mutliple of 2π. For negative values the circle is empty (black) while it is filled (white) for positive values. For the second component we make use of the same logic but instead of a clock hand, we only use a circle that orbits the image center. The two hands cannot obscure each other as the orbiting cirle uses a larger radius. Figure 8 shows three examples of pairs of corresponding latent vectors and generated images.

E.6 PM 2.5 air quality
Due to health concerns the air quality in Chinese cities has become an important topic. The PM 2.5 data set contains hourly meteorologic information and the concentration of small particles (PM 2.5 ) for the cities Beijing, Shanghai, Guangzhou, Chengdu and Shenyang (Liang et al., 2016). The individual features available for all cities are listed in Table 3. In addition to the features listed, the data includes the date and time of each record. Just like in the preprocessing of Traffic we compute a periodic time feature using expression 10 to represent the time of day of each record. For each numerical feature we calculate the mean and standard deviation and remove rows with values that are more extreme than six standard deviations from the mean. We also remove rows with missing categorical features, which are then represented as one-hot vectors. Apart from differences in the preprocessing we consider the same prediction task as Karlsson et al. (2022) which is predicting the future particle concentration for a fixed time horizon given current observations.

F Additional experiment results
In the course of this section we present a larger scope of our experimental results. We demonstrate the predictive accuracy of the algorithms introduced in Sections 3.2 and 3.3 in terms of the mean coefficient of determination R 2 over different settings on five data sets. The variation of the results over repetitions is represented by the shaded areas in the visualizations, which denotes one standard deviation above and below the mean value. We also show empirically that the bias of generalized LuPTS increases with the number of privileged time steps when a poor representation function is used. This can be seen in Figure 16a where we test privileged and classical learners over 500 different systems of the Square-Sign type over different sequence lengths. In addition we further illustrate how privileged information can improve the recovery of latent variables of the data generating process by providing more visualizations in the style of Figure 7a on the Clocks-LGS and Square-Sign data set. These can be found in Figures 13 and 16.
In the following we first provide the experiment details for visualizing the phase transition of generalized LuPTS as seen in Figure 2. In the subsections thereafter the material is organized by data set.

F.1 Two regimes of generalized LuPTS
As seen in Figure 2 and demonstrated by Proposition 1, generalized LuPTS becomes equivalent to the corresponding classical learner when the number of featuresd is larger than the number of samples m. In the following we provide the experiment details that led to Figure 2.
In order to evaluate the dependency on the number of features we used linear LuPTS and OLS on a synthetically generated linear-Gaussian system as displayed in Figure 1a. We use the same setup as for the latent dynamics in Square-Sign but without a nonlinear observation generating function. For each number of features k =d = d we sample 50 such systems with different dynamics, each producing a training set with m = 100 and test set of 1000 samples. The systems are all configured with T = 3 and q = 10. We train and score both estimators on all of the data generating systems before computing the mean coefficient of determination over all systems with the same number of features.