extfieldsofstudy:
Computer Science
Mathematics

text:
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.


extfieldsofstudy:
Computer Science
Engineering

text:
BERT, can HE predict contrastive focus? Predicting and controlling prominence in neural TTS using a language model

Several recent studies have tested the use of transformer language model representations to infer prosodic features for text-to-speech synthesis (TTS). While these studies have explored prosody in general, in this work, we look specifically at the prediction of contrastive focus on personal pronouns. This is a particularly challenging task as it often requires semantic, discursive and/or pragmatic knowledge to predict correctly. We collect a corpus of utterances containing contrastive focus and we evaluate the accuracy of a BERT model, finetuned to predict quantized acoustic prominence features, on these samples. We also investigate how past utterances can provide relevant information for this prediction. Furthermore, we evaluate the controllability of pronoun prominence in a TTS model conditioned on acoustic prominence features.

Introduction
"HE loves her", "he LOVES her", and "he loves HER" all have the same textual content, but three distinct communicative goals. Indeed, such contrastive focus is used by speakers to evoke alternative sets in the discourse [1]. This can be utilized to make explicit intended discourse relations between clauses/paragraphs/sections, to highlight a fact that the listener may find surprising, or to express a specific semantic or pragmatic meaning. The prediction of contrastive focus placement therefore often requires high-level linguistic understanding. Current vanilla neural text-to-speech (TTS) synthesis systems lack this understanding and will always pronounce the above sentences in the same way, irrespective of the context. In this work, we investigate methods to predict the placement of contrastive focus and to control it in a TTS system. Figure 1 illustrates our overall approach which addresses both predicting and controlling prominence.
One way to insert sophisticated linguistic information may be through the use of contextualized word embeddings. While other works have explored the use of transformer language models to predict prosodic and stylistic features [2,3,4,5,6,7,8], it has not been fully explored how much the encoded word representations actually imbue high-level knowledge. In other words, do they provide information about the content and context of the message or do they only provide/reinforce low-level linguistic features such as the likelihood of lexical prominence, parts of speech and position in the sentence? For prominence/pitch accent prediction, a fairly high baseline can be achieved using word majority/accent-ratio alone (i.e., if a lexical item is usually prominent in the training set, it is likely to be prominent in the test set) [9]. Moreover, [10] found that in a binary pitch accent prediction task, using broad word category distinctions (open/content or closed/function) could achieve 68% accuracy; more fine-grained division of the closed class category brought that number up to 77%. In this work, we probe a language model, in the present case BERT [11], by choosing a testing ground that cannot rely on simple heuristics to achieve good results: the prediction of contrastive focus on personal pronouns. Personal pronouns, at least in the corpora typically used to train TTS systems, are majority nonprominent. 1 Contrastive focus in English is mainly realized prosodically through increases in pitch, energy and/or duration. Suni et al. [12] have proposed extracting prosodic information automatically through a combined representation of these signals. They use continuous wavelet transforms, which analyze the speech signal at different timescales, to identify acoustically prominent words. In [4], they tested the prediction capabilities of a language model on these features. Here, we extend this investigation by looking at a challenging set of contrastively focused pronouns and by looking at the use of an extended context (current sentence, +previous sentence, +2 previous sentences). Moreover, because the prediction of contrastive focus is complicated by the one-to-many issue that hinders TTS evaluation in general, we have collected a corpus of audiobooks where we have three separate renditions of each book, each read by a separate speaker. We use these three realizations to study consensus among speakers in the placement of contrastive focus.
Assuming we are able to predict good candidates for contrastive focus, it remains to be verified if we can control the contrastive prominence on pronouns in a neural TTS system. Other works have proposed utterance-level control systems [13,14] but for this task we need to target the word-level. This was achieved previously within an HMM framework [15]. [3] evaluated the expression of prosodic features from ToBI-label conditioned TTS, but only for ground-truth input labels. [16] report successful control over other word categories from systems conditioned on human or automatically annotated data (although this was not evaluated with a perceptual test). More recently, [17] proposed to use control tags for prominence in an end-to-end TTS system. However, their system relied on the availability of a specific (and limited) pre-annotated corpus whereas we demonstrate controllability from prominence labels obtained automatically. Figure 1: TTS overview: The system is split into two modules. The first uses a language model to predict prominence labels. The second controls the prominence in the synthetic speech in accordance with the predicted labels.

Prepared Datasets
The selected audiobooks for our corpus are literary texts sourced from Librivox 2 and from the Blizzard Challenge 2013 dataset. 3 . Criteria for book selection included open-source status, the availability of multiple recordings with different speakers (minimum 3), audio quality and the subject matter (we favoured books dealing with interpersonal relationships as they are more likely to contain contrastively focused pronouns). Five novels (x 3 speakers) were selected for the training set (41,593 utterances, approx. 66 hours of audio/speaker 4 and one novel (x 3 speakers) was selected for testing (6838 utterances, approx. 11 hours of audio per speaker). The corresponding transcriptions were obtained from Project Gutenberg. 5 Transcripts were split into chapters and then sentences using [18] and [19]. The audio files were segmented into utterances using the Aeneas library 6 and phoneme alignment was obtained using the Montreal Forced Aligner [20]. Prosodic feature extraction. The audio files were analyzed with the continuous wavelet transform (CWT) method of [12] implemented in the Wavelet Prosody Toolkit. 7 This method assigns a prominence score to each word in the corpus. This is done by combining f0, energy and duration into a composite signal, performing the CWT, establishing lines of maximum amplitude connecting the various timescales and then calculating a weighted sum of the points in this line (see [12] for details). We then quantize this score into three categories: <p2> (strong prominence), <p1> (intermediate prominence) and <p0> (no prominence). 8 Contrastive personal pronoun subcorpus. With the processed data from the previous section, we searched for utterances containing <p2> labelled personal pronouns (strong prominence) in the test set. With manual verification, we collected positive examples of contrastive focus on pronouns. We randomly selected an equal number of negative samples where the pronoun was tagged as <p0> for all three speakers; these samples were also manually checked. We then enlisted three native English speakers to validate 200 pronouns (x 3 speakers) from the collected samples (100 positive samples and 100 negative samples, just over 20% of the full pronoun corpus). Validators were presented with the audio clips (for all three speakers) and a transcription of the text with the pronoun of interest highlighted in red. They were asked to assign a value of 1 if they deemed the speaker had used prosody to convey contrastive focus and  [21] were used to evaluate inter-annotator agreement between the three raters. These scores range from 0.85-0.90; this shows strong to almost perfect agreement. For evaluation purposes, we sorted the positive samples into two groups: 1) those where the majority of speakers (at least 2 out of 3) used contrastive focus (Pronoun maj.). This group contains 393 pronouns; 310 of which are prosodically contrasted by all three speakers; 2) those where only 1 out of 3 speakers used contrastive focus (Pronoun min.). This group contains 100 pronouns. All contrastive pronouns come from 406 utterances as several utterances contain multiple examples.
In this study, our intention is to find challenging examples for prominence prediction (i.e., words that are not frequently prominent). Subjective pronouns (e.g., I, we), objective pronouns (e.g., me, us) and possessive determiners (e.g., my, our) all fit this requirement, but possessive pronouns (e.g., mine, ours) are more often prominent than not. We decided to include these "easy" words in the corpus because they may be of interest for future work on contrastive focus. However, they do not present a particular challenge to prominence prediction.

Predicting Prominence
BERT. For our prominence prediction task, we used English BERT [11], available on HuggingFace. 9 BERT is a transformer encoder that is trained on a masked language modeling task (i.e., it learns to predict masked words using information from the other words in the sentence). Bert makes use of selfattention layers and positional embeddings to form contextualized representations of words. Prediction task. For the sequence of words {w0, w1, ..., wN } our objective is to predict a sequence of prominence labels {c0, c1, ..., cN }, where cn is either <p0>, <p1> or <p2>. Since knowledge about the previous context is sometimes es-   Figure 2 for an illustration of the model's components. Models and linguistic knowledge. We evaluate three methods for prominence prediction with increasing access to linguistic knowledge: • Our baseline is a simple word majority method: word statistics from the training corpus are computed; we count how often a lexical item belongs to each of the three prominence categories and the majority category is used for all predictions in the test set.
• The second method involves the use of the BERT architecture, but instead of using weights pretrained on a masked language modelling task, we randomly initialize the model and train it to predict the prominence labels for each word in the input sequence. This model can presumably learn the same word statistics available to the word majority method and additionally the self-attention layers and positional embeddings provide the model with information about the surrounding lexical items and the positional context of each word. We expect this model will be able to learn canonical patterns of English prosody (i.e. that prominent, nuclear accents are typically found at the end of an intonational phrase) even if the semantic knowledge about the content of the sentences will be non-existent or at best, very naive.
• The third method involves finetuning a pretrained BERT model on the prominence features (learning rate = 5e−5).
From the beginning of training, this model has access to the syntactic and semantic representations learned from training on huge amounts of textual data; during the finetuning process it must find the optimal way to use this information for the prosody/prominence prediction task. Results. The results for the full dataset are shown in Table 1 and the results for the pronoun subsets are shown in Table 2. Analyzing these results, we notice that performance on the contrastive pronoun sets is significantly lower than the full dataset (Recall with finetuned BERT is 0.239 for the pronoun majority group and 0.552 for all POS). The tokens in the pronoun minority group are very rarely predicted to be prominent (highest recall score= 0.07). Furthermore, we see that BERT does provide some improved prediction accuracy over the two baselines, but the improvement is fairly small. Word majority method correctly predicts possessives (e.g., mine, yours), and the randomly initialized BERT learns structurally/positionally prominent positions; it correctly predicts prominence at the ends of prosodic phrases (immediately preceding punctuation marks) and following the word 'for' (e.g., "For YOU are..."). We performed Mc-Nemar's Test [22] to compare models and we find significant differences (p-values < 0.05) between each of the three models. With regards to the use of previous context, we do not see any improvements in prediction performance (p-values > 0.05).
We can imagine several possible causes for low prediction scores. It may be that we have an insufficient number of samples of contrastively focused pronouns to train the model to recognize focus patterns; there is an average of 7893.6 <p2> labelled personal pronouns/speaker in the training set, and given the complexity of the task, this may not be enough. Alternatively, learnt representations may not be sophisticated enough to encode the higher level linguistic information required for this task. The recent research trend in language modeling is to scale models bigger and bigger and this increase in size results in better quality on tasks such as text generation. In future work, we will explore using larger LMs. Finally, as can be expected with any automatic annotation method, there is some noise in the data: we did find examples for which prominence was questionable, predominately at phrase boundaries (tagged <p2> because of a sharp rise in f0). Hence, human intervention may still be necessary for better fine-grained annotation/control of prosodic data.

Controlling Prominence
Controllable TTS Model. To synthesize speech with controllable prominence, we follow the method proposed in [16] where the TTS is conditioned on prominence labels. Our implementation differs in that we use FastSpeech 2 [23] (as implemented by [24]) instead of DC-TTS and we refrain from using boundary tags as we are primarily interested in prominence control. Fast-Speech 2 is a non-autoregressive, transformer encoder-decoder model that makes explicit duration, f0 and energy predictions between the encoder and decoder. The input to our model is a sequence of phonemes and prominence labels. Each word wn in the utterance is converted into a sequence of phonemes {pn,0, ..., pn,m} and this phoneme sequence is followed by the prominence label cn for wn (phonemes and prominence labels are converted into embeddings in the first layer of the model). The output of the modified FastSpeech 2 model is a Mel-spectrogram, which is converted into a waveform using a Parallel WaveGAN vocoder [25]. 10 To train and test the TTS model, we use the data described in Section 2, but only for a single speaker (Blizzard Challenge 2013; this speaker has read all 6 books). While the training corpus is read by a single speaker, this speaker portrays several different characters with different accents and pitch ranges. This tends to introduce fuzziness into the synthetic speech. To help the model learn these characteristics and improve quality, we added a speaker embedding to the TTS input. To obtain a speaker embedding, we 1) encoded each utterance in the training set with a pretrained speaker identification model (ECAPA-TDNN [26] available at [27]), 2) used k-means clustering on these embeddings to obtain 30 different 'speakers' and 3) used these speaker labels as an additional input to FastSpeech 2; the speaker embeddings are concatenated with the phoneme and prominence label embeddings. Listening Test. To test the controllability of our TTS model in terms of prominence, we conducted an ordinal ranking listening test using the Web Audio Evaluation Tool [28], following this procedure: 1) we randomly sampled 100 utterances containing pronouns from our full test set (not solely from the contrastive subset); 2) from this selection, we took the first 10 utterances containing a subjective pronoun, the first 10 with an objective pronoun and the first 10 with a possessive determiner (for a total of 30 utterances); 3) using our pretrained TTS system, we synthesized three versions of each utterance, changing only the prominence label for the relevant pronoun (c ∈ {<p0>,<p1>,<p2>}). The prominence labels for all other words in the sentence were kept constant with the ground truth values (extracted from the original audio with the CWT method); 4) 30 native English evaluators, recruited on Prolific, 11 were presented with the three versions of the synthesized utterances (in random presentation order) and a transcript of the audio with the pronoun of interest in uppercase letters. Participants were asked to rank the prominence of the pronoun by dragging and dropping the movable audio clips so that they were arranged from most prominent to least; 5) clips ranked as most prominent are assigned a score of 1; clips ranked as second most prominent are assigned a score of 0.5 and clips ranked as least prominent are given a score of 0. Sample audio files are available in the supplementary multimedia materials. Results. The results of the listening test are shown in Figure 3. We see the median values align with the prosody labels used (median: <p0>= 0, <p1>= 0.5, <p2>= 1.0). We do however see a wide distribution in the responses. This, and the examination of the ratings for individual utterances, indicates that this method works, but not consistently (i.e., there are some utterances for which the evaluators could not detect a difference). The only pronoun category for which we see a fairly clear distinction between <p0>, <p1>, and <p2> in perceived prominence is for possessive determiners. This may be because there is more natural variation within the training corpus for this category. Or, it may be due to the labelling errors at phrase bound- Figure 3: Results from the prominence ordinal ranking listening test. Audio ranked as most prominent was assigned a score of 1, second most prominent 0.5 and least prominent 0. Orange lines show the median and green triangles show the mean.
aries discussed in the previous section: the sampled subjective and objective pronouns were found more often at the beginning and end of phrases than the possessive determiners. More work is thus needed to disentangle the global prosodic representations from that of individual words, but this separation is difficult because it may result is less natural utterances.

Conclusion and Perspectives
In this paper, we investigated the difficult tasks of prominence prediction and control for contrastively focused personal pronouns. Prediction: we compared models with varying degrees of access to linguistic knowledge and past context on a wordlevel prominence label prediction task. We found that a finetuned BERT gave the best prediction performance, but that the improvement over baselines was very small. In the future, we will investigate the use of word representations from larger language models with more sophisticated linguistic understanding. Control: with a perceptive test, we evaluated the control of prominence in a TTS model conditioned on prominence labels. The results show the model is able to provide some control but the performance is not consistent over all samples. Natural variation: we used multiple spoken versions of the same written text to see the agreement among speakers on the use of contrastive focus. But we must keep in mind that while the textual context remains the same, the interpretation of the text can vary. For example, we infer that one of the speakers interprets some of the characters in the novel to be passive aggressive and they convey this through the frequent use of contrastive focus on 'I' (e.g., I am going to Gretna Green (intended meaning: and YOU are not). Removing the contrastive focus here is not wrong, but gives a very different impression of character/situation/relationship. This illustrates the difficulty of the prediction task and therefore, depending on the intended usage of the TTS system, it might be fruitful to explore other sources of input to the prominence prediction model (e.g., the source speech in a speech-to-speech translation system) in order to be as faithful to the intended meaning as possible.


extfieldsofstudy:
Mathematics
Computer Science

text:
A branch and bound technique for finding the minimal solutions of the linear optimization problems subjected to Lukasiewicz

In this paper, an optimization model with a linear objective function subject to a system of fuzzy relation equations (FRE) is studied where the feasible region is defined by the Lukasiewicz t-norm. Since the finding of all minimal solutions is an NP-hard problem, designing an efficient solution procedure for solving such problems is not a trivial job. Firstly, the feasible domain is characterized and then the problem is solved with a modified branch-and-bound solution technique based on a new solution set that includes the minimal solutions. After presenting our solution procedure, a concrete example is included for illustration purposes.

I. INTRODUCTION
In this paper, we study the following linear problem in which the constraints are formed as fuzzy relational equations defined by Lukasiewicz t-norm: (1) Where , , is a fuzzy matrix such that ( and ), is a fuzzy vector such that ( ), and " " is the Lukasiewicz t-norm defined by . If is the 'th row of matrix , then problem (1) can be expressed as follows: where the constraints mean: The theory of fuzzy relational equations (FRE) was firstly proposed by Sanchez and applied in problems of the medical diagnosis [13]. Nowadays, it is well known that many issues associated with a body knowledge can be treated as FRE problems [12]. The solution set of FRE is often a non-convex set that is completely determined by one maximum solution and a finite number of minimal solutions [4]. The other bottleneck is concerned with detecting the minimal solutions that is an NP-hard problem [5,7,8,10]. The problem of optimization subject to FRE and FRI is one of the most interesting and on-going research topics among the problems related to FRE and FRI theory [1,2,[4][5][6][7][8]11,14]. Recently, many interesting generalizations of the linear programming subject to a system of fuzzy relations have been introduced and developed based on composite operations used in FRE, fuzzy relations used in the definition of the constraints, some developments on the objective function of the problems and other ideas [3,[7][8][9]11].
In this paper, after carefully study the solution set of system (1), we show that the problem can be converted to a 0-1 integer programming problem and solved by a branch-andbound method.
The rest of the paper is arranged as follows. In Section 2, a special characterization of the feasible domain of problem (1) is derived. Based on the theoretical aspects of the problem, a new solution set is obtained that includes all the minimal solutions. In Section 3, we study the effect of the cost vector c with the special characterization. A modified branch-andbound method is presented and a step-by-step algorithm for solving problem (1) is given in Section 4. Some examples are included in Section 5 to illustrate how the algorithm works. Some concluding remarks are made in Section 6.

III. LINEAR PROGRAMMING AND EFFECT OF COST VECTOR
According to the well-known schemes used for optimization of linear problems such as (1) [5,10], problem (1) is converted to the following two sub-problems: And, where and for . It is easy to prove that is the optimal solution of , and the optimal solution of is for some .

Theorem 2.
Suppose that and and are the optimal solutions for and , respectively. Then is an optimal solution of problem (1), where is defined as follows: for .

Proof.
Let . Then, from Corollary 4 we have . Therefore, for each such that , inequality implies . In addition, for each such that , inequality implies . Hence, .  As mentioned earlier, generating the maximum solution is not a problem. If we know how to find a minimum solution for , then problem (1) can be solved via Theorem 2. In the next section, we provide a branch and bound method for solving .

IV. MODIFIED BRANCH AND BOUND METHOD
A branch-and-bound method implicitly enumerates all possible solutions to a programming problem. For our application, in the beginning, we choose one constraint to branch the original problem into several sub-problems. Each sub-problem is represented by one node. Then branching at each node is done by adding one additional constraint. New sub-problems are created and represented by new nodes. Note that the more constraints added to a sub-problem, the smaller feasible domain it has and, consequently, the larger optimal objective value Z it achieves. Therefore, solving one subproblem could provide the possibility to eliminate many possible solutions from consideration. In other words, once a current candidate solution is obtained, we judge other nodes for further consideration. If the best potential solution of one particular node is no better than the current candidate solution, then there is no need to branch on this node. Otherwise, branching is needed to yield a new bound. In this section, we use one concrete example to illustrate how this method works.
Based on the theory we built in previous sections, here we propose an algorithm for finding an optimal solution of problem (1).
Step 1: Find the maximum solution of system (1). Compute according to relation (2).
Step 4: Arrange cost vector. Define and for and define sub-problems and .
Step 5: Use the branch-and-bound concept based on to solve .
Step 3: We compute , , and . So, the cardinality of and is equal to 60. However, in this example we have also , , and . Therefore, the cardinality of and is equal to 8.
Step 4: The main problem is converted to the following two sub-problems: Step 5: Our branch-and-bound method begins by choosing from . Hence, each feasible solution must satisfy either or . This yields two branches denoted by nodes 1 and 2 in Fig. 1. We then obtain a lower bound on the Z-value associated with each node. For example, when , we know , therefore, any choice with results in . So we write for node 1 of Fig. 1. Similarly, any choice with results in . Since at this moment we have no reason to exclude any one of nodes 1 and 2 from consideration, both nodes have to be further investigated. By using the jump-tracking technique, we branch on the node with a lower bound on Z. In this ease, we choose node 2. Since , any choice associated with node 2 must satisfy either or . Branching on node 2 yields nodes 3 and 4 in Fig. 1. For each new node, we need to evaluate a lower bound for the objective value Z. For example, at node 3, we compute Z . Hence, any choice associated with node 3 must have Z . Similar reasoning shows that Z for node 4. Again, we do not have any reason to exclude any of nodes 3 and 4 from consideration, so we need to branch on one node. The jumptracking technique directs us to branch node 3. Any choice associated with node 3 must satisfy either or , since . This yields nodes 5 and 6 in Fig. 1. By the same reasoning as before, any choice associated with node 5 must have Z and any choice associated with node 6 must have Z . Of course, we are interested in node 5. To branch further on node 5, any choice associated with node 5 must have , since . This yields node 7 in Fig. 1. Note that node 7 corresponds to the sequence . This sequence leads to a value of Z . Therefore, node 7 is a feasible sequence which may be viewed as a candidate solution with Z . Because the Z value for nodes 4 and 6 cannot be lower than 2.75, these nodes can be eliminated from further consideration. However, node 1 cannot be eliminated yet, because it is still possible for node 1 to yield a sequence having Z . Hence, we now branch on node 1. Since , any sequence associated with node 1 must have either or . Correspondingly, we create nodes 8 and 9, respectively. For node 9, we calculate Z . Since any sequence associated with node 9 must have Z , this node is eliminated from consideration. For node 8, we have Z . So it cannot be eliminated yet. We branch on node 8. Since , any sequence associated with node 8 must have either or . This creates nodes 10 and 11 in Fig. 1. For node 10, since Z , we eliminate it. Similarly, for node 11 with Z , it can be eliminated too. Now, with the exception of node 7, every node in Fig. 1 has been eliminated from consideration. Node 7 yields the sequence that generates solution . Thus, we should choose the sequence and solution with Z as an optimal solution for .
Step 6: According to formula (5), we can find the optimal value of the original problem as with optimal objective value .

V. NUMERICAL EXAMPLES
In this section, we present the experimental results for evaluating the performance of our algorithm. For this purpose, we apply our algorithm to the following eight test problems. The test problems have been randomly generated in different sizes.

Test Problem 1:
[-9.3708 , 6.7080 , 6.        Table 1 shows the cardinality of and , the cardinality of and , the number of nodes that the algorithm visits and the number of complete paths visited by the algorithm. Moreover, for each test problem , the best path , the optimal solution and the optimal value is listed as follows: In this paper, we have studied a linear optimization problem subjected to a system of Lukasiewicz fuzzy relation equations and presented a procedure to find an optimal solution. Due to the non-convexity nature of its feasible domain, we tend to believe that there is no polynomial-time algorithm for this problem. The best we can do here is that, after analyzing the properties of its feasible domain, we find a solution set that include all the minimal solutions of the main feasible region, then apply a modified branch-and-bound method to find one solution. The question of how to generate the whole optimal solution set is yet to be investigated. From the analysis of Theorem 2, it is clearly seen that if all minimum solutions of a given system of fuzzy relation equations can be found, then an optimal solution of the optimization problem defined by (5) can be constructed. Therefore, solving this optimization problem is no harder than solving a system of fuzzy relation equations for all minimum solution. Although it is not known whether these two problems are essentially equivalent or not, the basic concepts introduced in Section 2 have been further developed for solving fuzzy relation equations. Extension to other types of objective functions is under investigation.