Functional Principal Component Analysis for Extrapolating Multi-stream Longitudinal Data

TL;DR

This paper introduces a non-parametric functional principal component analysis method for predicting multi-stream longitudinal data by leveraging historical data and real-time updates, outperforming existing approaches.

Contribution

It presents a novel Gaussian process-based framework that combines functional PCA with Bayesian updating for accurate real-time predictions of multi-stream data.

Findings

Outperforms state-of-the-art methods in synthetic and real data

Effectively captures heterogeneity among units

Achieves high predictive accuracy with real-time data

Abstract

The advance of modern sensor technologies enables collection of multi-stream longitudinal data where multiple signals from different units are collected in real-time. In this article, we present a non-parametric approach to predict the evolution of multi-stream longitudinal data for an in-service unit through borrowing strength from other historical units. Our approach first decomposes each stream into a linear combination of eigenfunctions and their corresponding functional principal component (FPC) scores. A Gaussian process prior for the FPC scores is then established based on a functional semi-metric that measures similarities between streams of historical units and the in-service unit. Finally, an empirical Bayesian updating strategy is derived to update the established prior using real-time stream data obtained from the in-service unit. Experiments on synthetic and real world data show that the proposed framework outperforms state-of-the-art approaches and can effectively account for heterogeneity as well as achieve high predictive accuracy.

Figures

Tables

Table 1. Table 1: Mean and standard deviation (STD.) of comparative models performed on the NASA dataset. The values for the sensor 15 are scaled by × 10 − 2 absent superscript 10 2 \times 10^{-2} . F-GP and F-B indicates FPCA-FP and FPCA-B, respectively. The best result in each case is boldfaced.

	Sensor 4			Sensor 15 ($\times {10}^{-2}$)
Model	25%	50%	75%	25%	50%	75%
F-GP	3.26	3.21	3.19	1.62	1.62	1.57
(std.)	(0.36)	(0.42)	(0.48)	(0.18)	(0.19)	(0.33)
F-B	3.49	3.37	3.31	1.76	1.75	1.63
(std.)	(0.45)	(0.56)	(0.54)	(0.28)	(0.31)	(0.36)
ME	3.51	3.38	3.34	1.79	1.77	1.65
(std.)	(0.45)	(0.59)	(0.55)	(0.28)	(0.32)	(0.37)

Equations

$$G(t_1,t_2)=\sum _{k=1}^{\infty }{\lambda }_k{\varphi }_k(t_1){\varphi }_k(t_2),$$

$$X(t)-\mu (t)=\sum _{k=1}^{\infty }{\xi }_k{\varphi }_k(t)+ϵ(t),$$

$$X_i(t)=\mu (t)+f_i(t)+ϵ(t),i\in I^{\mathcal{H}},$$

$$X_i(t)=\mu (t)+\sum _{k=1}^K{\xi }_{ik}{\varphi }_k(t)+ϵ(t),i\in I^{\mathcal{H}}.$$

$$\begin{split}[&\bm{\xi^{\prime}}_{k},\xi_{rk}]^{\prime}\\ &=\mathcal{G}_{k}\big{(}\bm{X}^{(-s)}_{1}(t),...,\bm{X}^{(-s)}_{N}(t),\bm{X}^{(-s)}_{r}(t)\big{)}+\bm{\epsilon}_{\xi_{k}},\end{split}$$

$$\begin{bmatrix}\bm{\xi}_{k}\\ \xi_{rk}\end{bmatrix}\sim\mathcal{N}\bigg{(}\mathbf{0}_{N+1},\begin{bmatrix}\mathbf{C}_{k}+\sigma_{k}^{2}\mathbf{I}_{N}&\mathbf{c}_{k}\\ \mathbf{c}^{\prime}_{k}&c^{(r)}_{k}\end{bmatrix}\bigg{)},$$

$$h(i,j;\bm{\theta}_{k})=\alpha_{k}\exp\bigg{(}-\frac{1}{2}\sum_{l\in L^{-s}}\frac{\big{\|}X^{(l)}_{i}(t)-X^{(l)}_{j}(t)\big{\|}^{2}_{l}}{\big{(}\beta^{(l)}_{k}\big{)}^{2}}\bigg{)},$$

$$\displaystyle\big{\|}$$

$$\displaystyle=\sqrt{\sum_{k=1}^{K(l)}\bigg{(}\int_{\mathcal{T}_{t^{*}}}\big{[}X^{(l)}_{i}(t)-X^{(l)}_{j}(t)\big{]}\psi^{(l)}_{k}(t)dt\bigg{)}^{2}},$$

$$P({\mathbf{\xi }}_k$$

$$=\int P({\mathbf{\xi }}_k|{\mathbf{z}}_k,\mathcal{X},{\mathbf{\Theta }}_k)P({\mathbf{z}}_k|\mathcal{X},{\mathbf{\Theta }}_k)d{\mathbf{z}}_k$$

$$=\int \prod _{i\in I^{\mathcal{H}}}P({\xi }_{ik}|z_{ik},\mathcal{X},{\mathbf{\Theta }}_k)P({\mathbf{z}}_k|\mathcal{X},{\mathbf{\Theta }}_k)d{\mathbf{z}}_k$$

$$\log P({\mathbf{\xi }}_k|\mathcal{X},{\mathbf{\Theta }}_k)$$

$$\displaystyle\qquad\>-\log\big{|}\mathbf{C}_{k}+\sigma_{k}^{2}\mathbf{I}_{N}\big{|}-\frac{n}{2}\log 2\pi,$$

$${\mathbf{\Theta }}_k^{\ast }=[{\mathbf{\theta }}_k^{{\ast }^{\prime }},{\sigma }_k^{\ast }{]}^{\prime }=\underset{{\mathbf{\Theta }}_k}{\mathrm{argmax}}\log P({\mathbf{\xi }}_k|\mathcal{X},{\mathbf{\Theta }}_k).$$

$$P(\xi_{rk}|\bm{\xi}_{k},\mathcal{X},\bm{\Theta}_{k}^{*})=\mathcal{N}\big{(}\hat{\xi}_{rk},\hat{\sigma}_{rk}^{2}\big{)}$$

$${\hat{\xi }}_{rk}={\hat{\mathbf{c}}}_k^{\prime }({\hat{\mathbf{C}}}_k+({\sigma }_k^{\ast }{)}^2{\mathbf{I}}_N{)}^{-1}{\mathbf{\xi }}_k,$$

$${\hat{\sigma }}_{rk}^2={\hat{c}}_k^{(r)}-{\hat{\mathbf{c}}}_k^{\prime }({\hat{\mathbf{C}}}_k+({\sigma }_k^{\ast }{)}^2{\mathbf{I}}_N{)}^{-1}{\hat{\mathbf{c}}}_k.$$

$$\begin{array}{rl}& {\hat{X}}_r(t)=\mu (t)+\sum _{k=1}^K{\hat{\xi }}_{rk}{\varphi }_k(t),\\ & {\hat{\sigma }}_r^2(t)={\sigma }_{\mu }^2(t)+\sum _{k=1}^K{\hat{\sigma }}_{rk}^2{\varphi }_k^2(t)+{\sigma }_{ϵ}^2.\end{array}$$

$$[{\xi }_{r1}^{\ast },...,{\xi }_{rK}^{\ast }{]}^{\prime }\sim \mathcal{N}({\mathbf{\xi }}^{\ast },{\mathbf{\Sigma }}^{\ast })$$

$$\displaystyle\bm{\xi}^{*}=\bm{\Sigma}^{*}\bigg{(}\bm{\Sigma}_{0}^{-1}\bm{\mu}_{0}+\frac{1}{\sigma_{\epsilon}^{2}}\bm{\Phi}(\mathbf{t})^{\prime}(\bm{X}_{r}(\mathbf{t})-\bm{\mu}(\mathbf{t}))\bigg{)},$$

$$\displaystyle\bm{\Sigma}^{*}=\bigg{(}\frac{1}{\sigma_{\epsilon}^{2}}\bm{\Phi}(\mathbf{t})^{\prime}\bm{\Phi}(\mathbf{t})+\bm{\Sigma}_{0}^{-1}\bigg{)}^{-1}$$

$$\begin{split}&\bm{\mu}_{0}=\big{[}\hat{\xi}_{r1},...,\hat{\xi}_{rK}\big{]}^{\prime},\\ &\bm{\Sigma}_{0}=\text{diag}\big{(}\hat{\sigma}_{r1}^{2},...,\hat{\sigma}_{rK}^{2}\big{)},\\ &\bm{X}_{r}(\mathbf{t})=[X_{r}(t_{r1}),...,X_{r}(t_{rp_{r}})]^{\prime},\\ &\bm{\mu}(\mathbf{t})=[\mu(t_{r1}),...,\mu(t_{rp_{r}})]^{\prime},\end{split}$$

$$\mathbf{\Phi }(\mathbf{t})=\left[\begin{array}{ccc}{\varphi }_1(t_{r1})& \dots & {\varphi }_K(t_{r1})\\ ⋮& \ddots & ⋮\\ {\varphi }_1(t_{r{p}_r})& \dots & {\varphi }_K(t_{r{p}_r})\end{array}\right].$$

$${\hat{X}}_r^{\ast }(\tilde{t})=\mu (\tilde{t})+\sum _{k=1}^K{\xi }_{rk}^{\ast }{\varphi }_k(\tilde{t}).$$

$$\displaystyle\big{(}$$

$$\displaystyle\qquad\sigma^{2}_{\mu}(t)+\sum_{k_{1}=1}^{K}\sum_{k_{2}=1}^{K}\big{(}\bm{\Sigma}^{*}\big{)}_{k_{1}k_{2}}\phi_{k_{1}}(\tilde{t})\phi_{k_{2}}(\tilde{t})+\sigma_{\epsilon}^{2},$$

$$\text{MAE}=\frac{1}{U}\sum_{u=1}^{U}\big{|}\hat{X}^{(s)}_{r}(t_{u})-X^{T}_{r}(t_{u})\big{|},\ t_{u}\in[t^{\ast},b].$$

Full text

Functional principal component analysis, Gaussian process

1 Introduction

Among various environments where longitudinal data is gathered, the environment covered in this study is a multi-stream and real-time environment. Recent progress in sensor and data storage technologies has facilitated data collection from multiple sensors in real-time as well as the accumulation of historical signals from multiple similar units during their operational lifetime. This data structure where multiple signals across different units are collected is referred to as multi-stream longitudinal data. Examples include: vital health signals from patients collected through wearable devices (Caldara et al., 2014; Magno et al., 2016), battery degradation signals from cars on the road (Meeker & Hong, 2014; Salamati et al., 2018) and energy usage patterns from different smart home appliances (Hsu et al., 2017).

In this article, we propose an efficient approach to extrapolate multi-stream data for an in-service unit through borrowing strength from other historical units. An illustrative example is provided in Figure 1. In this figure, there are $N$ historical units and an in-service unit whose index is denoted by $r$. Each unit has $M$ identical sensors from which each respective signal forms a stream. Multi-stream data from the in-service unit is partially observed up to the current time instance $t^{\ast}$. Our goal is to extrapolate stream data from the in-service unit $r$ over a future period $t\geq t^{*}\in\mathcal{T}$ where $\mathcal{T}$ is the time domain of interest.

In mathematical notation, let $I=\{1,...,N,r\}$ and $I^{\mathcal{H}}=\{1,...,N\}$ be the respective index sets for all available units including the in-service unit $r$ and the units in our historical dataset. For each unit $i\in I$, we have $M$ streams of data where $l\in L=\{1,...,M\}$. For the $i$th unit, the history of observed data for a specific stream $l$ is denoted as $X^{(l)}_{i}(t_{iu})$, where $\{t_{iu}:u=1,...,p^{(l)}_{i}\}\subset\mathcal{T}$ represents the observation time points and $p^{(l)}_{i}$ represents the number of observations for signal $l$ of unit $i$. The underlying principle of our model is borrowing strength from a sample of curves $\{X_{i}^{(l)}(t):i\in I^{\mathcal{H}},l\in L\}$ to predict individual trajectories $X_{r}^{(l)}(t)$ over a future time period $t\geq t^{*}\in\mathcal{T}$. Without loss of generality, throughout the article we focus on predicting stream $s\in L$, which we refer to as the target stream. Note that the target stream in Figure 1 is stream $M$. To achieve this goal we exploit functional principal component analysis (FPCA), which is a non-parametric tool for functional data analysis (Ramsay & Silverman, 2005). Indeed, FPCA has recently drawn increased attention due to its flexibility, uncertainty quantification capabilities and the ability to handle sparse and irregular data (Peng & Paul, 2009; Di et al., 2009; Wang et al., 2016; Xiao et al., 2018). However, advances in FPCA fall short of handling multi-stream data and real-time predictions.

Our overall framework is summarized at the bottom of Figure 1. Specifically, historical signals $X^{(s)}_{i}(t)$, $i\in I^{\mathcal{H}}$, from the target stream $s$ are decomposed into a linear combination of orthonormal eigenfunctions that form their functional space. The coefficients of the linear combination are called functional principal component (FPC) scores. With the assumption that the target signal of the in-service unit $X^{(s)}_{r}(t)$ lies in the same functional space, a proper estimation of FPC scores associated with $X^{(s)}_{r}(t)$ is required. Here, we propose to establish a prior on these FPC scores using information from streams $l\in L^{-s}=L\backslash\{s\}$. Specifically, a Gaussian process (GP) prior for FPC scores of $X^{(s)}_{r}(t)$ is built using a functional semi-metric that measures similarities of streams $L^{-s}$ between historical units and the in-service unit $r$. The underlying principle is that unit $X^{(s)}_{r}(t)$ will exhibit more commonalities with historical units that exhibit similar trends in streams $L^{-s}$. For example, if stream $s$ denotes degradation trajectories and $L^{-s}$ denotes external factors such as temperature. Then $X^{(s)}_{r}(t)$ will share more commonalities with a subset of historical signals $X^{(s)}_{i}(t)$, $i\in I^{\mathcal{H}}$, degrading under similar external factors (similar temperatures). This approach allows us to address heterogeneity in the data. Lastly, an empirical Bayesian updating strategy is derived to update the established prior using real-time stream data obtained from the in-service unit.

2 Literature Review

There has been extensive literature on the extrapolation of longitudinal signals under a single stream setting. However, literature has mainly focused on parametric models due to their computational efficiency and ease of implementation (Gebraeel et al., 2005; Gebraeel & Pan, 2008; Si et al., 2012, 2013; Kontar et al., 2017). Such models have been applied in healthcare, manufacturing and mobility applications specifically to understand the remaining useful life of operational units. Unfortunately, in real-world applications, parametric modeling is vulnerable to model misspecifications, and if the specified form is far from the truth, predictive results will be misleading. For instance, parametric representations are specifically challenging when data is sparse or when the underlying physical and chemical theories guiding the process are unknown.

To address this issue, recent attempts at non-parametric approaches have been based on FPCA (Zhou et al., 2011, 2012; Fang et al., 2015) or multivariate Gaussian processes (Álvarez & Lawrence, 2011; Saul et al., 2016; Kontar et al., 2018, 2019). These studies show that such non-parametric approaches outperform parametric models in case where functional forms are complex and exhibit heterogeneity. Nevertheless, the foregoing works have dealt with only single stream cases.

On the other hand, the few literature that addressed multi-stream settings have focused on data fusion approaches. Data fusion in this case refers to aggregating all streams into a single stream using fusion mechanisms. In health related applications, this fused stream is coined as a health-index which is often derived through a weighted combination of the $M$ data streams (Liu et al., 2013; Song & Liu, 2018). Such methods require regularly sampled observations and enforce strong parametric assumptions. An alternative data fusion approach includes multivariate FPCA (Fang et al., 2017a, b). However, since data fusion methods are operated by aggregating multi-streams into a single or a smaller group of streams, they are not capable of predicting individual stream trajectories and thus have limited applications.

Compared to current literature, our contribution can be summarized as follows. We propose an FPCA-based model that provides individualized predictions in a multi-stream environment. Our model is able to automatically account for heterogeneity in the data and screen the sharing of information between the in-service unit and units in our historical dataset. We then derive a computationally efficient Bayesian updating strategy to update predictions when data is collected in real-time. We demonstrate the advantageous features of our approach compared to state-of-the-art methods using both synthetic and real-world data.

The rest of this paper is structured as follows. In section 3, we briefly revisit the FPCA. In section 4, we discuss our proposed model. Numerical experiments using synthetic data and real-world data are provided in section 5. Finally, section 6 discusses the computational complexity of our model. Technical proofs, a detailed code and additional numerical results are available in the supplementary materials.

3 Brief Review of FPCA

From an FPCA perspective, longitudinal signals observed in a given time domain can be decomposed into a linear combination of orthonormal basis functions with corresponding FPC scores as coefficients. Therefore, FPCA can be regarded as a dimensionality reduction method in which a signal corresponds to a vector in a functional space defined by the basis functions. The basis functions are referred to as eigenfunctions. Let us assume that the longitudinal signals, over a given time domain $\mathcal{T}$, are generated from a square-integrable stochastic process $X(t)$ with its mean $\mathbb{E}[X(t)]=\mu(t)$ and covariance defined by a positive semi-definite kernel $G(t_{1},t_{2})=\text{Cov}(X(t_{1}),X(t_{2}))$ for $t_{1},t_{2}\in\mathcal{T}$. Using Mercer’s theorem on $G(t_{1},t_{2})$, we have

[TABLE]

where $\phi_{k}(t)$ presents the $k$th eigenfunction of the linear Hilbert-Schmidt operator $G:L^{2}(\mathcal{T})\rightarrow L^{2}(\mathcal{T}),G(f)=\int_{\mathcal{T}}G(t_{1},t_{2})f(t_{1})dt_{1}$ ordered by the corresponding eigenvalues $\lambda_{k}$, $\lambda_{1}\geq\lambda_{2}\geq\dots\geq 0$. The eigenfunctions $\phi_{k}(t)$ form a set of orthonormal basis in the Hilbert space $L^{2}(\mathcal{T})$. Following the Karhunen-Loéve decomposition, the centered stochastic process $X(t)-\mu(t)$ can then be expressed as

[TABLE]

where \xi_{k}=\int_{\mathcal{T}}\big{(}X(t)-\mu(t)\big{)}\phi_{k}(t) presents FPC scores associated with $\phi_{k}(t)$. The scores are uncorrelated normal random variables with zero-mean and variance $\lambda_{k}$; that is, $\mathbb{E}[\xi_{k}]=0,\forall k\in\mathbb{N}$ and $\mathbb{E}[\xi_{k_{1}}\xi_{k_{2}}]=\delta_{k_{1}k_{2}}\lambda_{k_{1}},\forall k_{1},k_{2}\in\mathbb{N}$ where $\delta_{k_{1}k_{2}}$ denotes the Kronecker delta. Also, $\epsilon(t)$ is additive Gaussian noise.

This idea of projecting signals onto a functional space spanned by eigenfunctions was first introduced by Rao (1958) for growth curves in particular. Basic principles (Castro et al., 1986; Besse & Ramsay, 1986) and theoretical characteristics (Silverman, 1996; Boente & Fraiman, 2000; Kneip & Utikal, 2001) were then developed. These ideas were expanded to longitudinal data settings in the seminal work of Yao et al. (2005). After that, the FPCA was applied and extended to a wide variety of applications, where multiple works tackled fast and efficient estimation of the underlying covariance surface (Huang et al., 2008; Di et al., 2009; Peng & Paul, 2009; Goldsmith et al., 2013; Xiao et al., 2016).

4 Extrapolation of Multivariate Longitudinal Data

4.1 FPCA for Signal Approximation

Now we discuss our proposed non-parametric approach for extrapolation of multi-stream longitudinal data. Hereon, unless there is ambiguity, we suppress subscripting the target stream with $(s)$. Using historical signals $X_{i}(t)$, $i\in I^{\mathcal{H}}$, we decompose the target stream $s$ as

[TABLE]

where $f_{i}(t)$ represents random effects characterizing stochastic deviations across different historical signals in stream $s$ and $\epsilon(t)$ denotes additive noise. We assume $f_{i}(t)$ and $\epsilon(t)$ are independent. Through an FPCA decomposition, we have that $f_{i}(t)=\sum_{k=1}^{\infty}\xi_{ik}\phi_{k}(t)$. This decomposition is an infinite sum, however, only a small number of eigenvalues are commonly significantly non-zero. For these values the corresponding scores $\xi_{ik}$ will also be approximately zero. Therefore, we approximate this decomposition as $f_{i}(t)=\sum_{k=1}^{K}\xi_{ik}\phi_{k}(t)$, where $K$ is the number of significantly nonzero eigenvalues.

[TABLE]

Here we follow the standard estimation procedures in Di et al. (2009) and Goldsmith et al. (2018) to estimate the model parameters where $\mu(t)$ is obtained by local linear smoothers (Fan & Gijbels, 1996), while $K$ is selected to minimize the modified Akaike criterion. Now given that the in-service unit $r$ lies in the same functional space spanned by $\phi_{k}(t)$, our task is to find the individual distribution of $\xi_{rk}$ using the partially observed multi-stream data from unit $r$. Specifically, we aim to find $X_{r}(t)=\mu(t)+\sum_{k=1}^{K}\xi_{rk}\phi_{k}(t)+\epsilon(t)$.

4.2 Estimation for Prior Distribution of FPC scores via GP

Next, we estimate the prior distribution of $\xi_{rk}$ based on the key premise that $X_{r}(t)$ will behave more similarly to $X_{i}(t)$ for some units $i\in I^{\mathcal{H}}$ whose signals $X^{(l)}_{i}(t)$ for $l\in L^{-s}$ are similar to the corresponding signals of the in-service unit $X^{(l)}_{r}(t)$. To this end, for $k\in\{1,...,K\}$, we model a functional relationship between $[\xi_{1k},...,\xi_{Nk},\xi_{rk}]^{\prime}$ and $X_{i}^{(l)}(t),...,X_{N}^{(l)}(t),X_{r}^{(l)}(t)$ for $l\in L^{-s}$ as

[TABLE]

where $\bm{X}^{(-s)}_{i}(t)=[X_{i}^{(l)}(t)]_{l\in L^{-s}}$ for $i\in I$, $\bm{\xi}_{k}=[\xi_{ik}]_{i\in I^{\mathcal{H}}}$, and $\bm{\epsilon}_{\xi_{k}}\sim\mathcal{N}(\mathbf{0},\sigma^{2}_{k})$.

The idea here, is to model $\mathcal{G}_{k}$ as a GP with a covariance function defined by a similarity measure between the observed signals, i.e., a functional similarity measure. Specifically, for any $k\in\{1,...,K\}$, the vector of FPC scores $[\bm{\xi}^{\prime}_{k},\xi_{rk}]^{\prime}$ will follow a multivariate Gaussian distribution

[TABLE]

where $\mathbf{C}_{k}\in\mathbb{R}^{N\times N}$ is constructed such that its $(i,j)$th element is $h(i,j;\bm{\theta}_{k})$ for $i,j\in I^{\mathcal{H}}$, $\mathbf{c}_{k}=[h(i,r,\bm{\theta}_{k})]_{i\in I^{\mathcal{H}}}$, $c^{(r)}_{k}=h(r,r,\bm{\theta}_{k})$, and $h$ denotes a covariance function defined as

[TABLE]

in which $\|\cdot\|_{l}$ is a semi-metric providing as similarity measure across functions, and $\alpha_{k}$ and $\beta_{k}^{(l)}$ are hyperparameters for streams $l\in L^{-s}$. For notational simplicity, we introduce \bm{\theta}_{k}=\big{[}\alpha_{k},\beta_{k}^{(1)},...,\beta_{k}^{(s-1)},\beta_{k}^{(s+1)},...,\beta_{k}^{(M)}\big{]}^{\prime}.

To show the validity of the GP (4), we provide the following lemma.

Lemma 1**.**

The matrix $\begin{bmatrix}\mathbf{C}_{k}+\sigma_{k}^{2}\mathbf{I}_{N}&\mathbf{c}_{k}\\ \mathbf{c}^{\prime}_{k}&c^{(r)}_{k}\end{bmatrix}$ corresponding to the covariance function $h(i,j;\bm{\theta}_{k})$ a valid covariance matrix.

Proof.

See Section A.1 in supplementary document. ∎

One possible semi-metric that represents the similarity between two signals can be derived based on FPCA. Let $\mathcal{T}_{t^{*}}\subset\mathcal{T}$ denote the time domain for observations up to $t^{*}$. Note that we define $\mathcal{T}_{t^{*}}$ since the signals of the in-service unit $r$ are available only for $t\in\mathcal{T}_{t^{*}}$. For $i,j\in I$, $l\in L^{-s}$, and $t\in\mathcal{T}_{t^{*}}$, the semi-metric based on FPCA for two signals $X^{(l)}_{i}(t)$ and $X^{(l)}_{j}(t)$ can be represented as

[TABLE]

where $\psi_{k}^{(l)}(t)$ is $k$th eigenfunction derived by the FPCA on $X^{(l)}_{i}(t)$ for $i\in I$ and $t\in\mathcal{T}_{t^{*}}$, and $K(l)$ is the number of eigenfunctions. We would like to point out that \int_{\mathcal{T}_{t^{*}}}\big{[}X^{(l)}_{i}(t)-X^{(l)}_{j}(t)\big{]}\psi^{(l)}_{k}(t)dt is the difference between the FPC scores of $X^{(l)}_{i}(t)$ and $X^{(l)}_{j}(t)$ associated with $\psi^{(l)}_{k}(t)$, which implies that this metric measures the Euclidean distance between two vectors composed of the corresponding FPC scores.

In order to optimize the hyperparameter $\bm{\Theta}_{k}=[\bm{\theta}^{\prime}_{k},\sigma_{k}]^{\prime}$ for the multivariate Gaussian distribution (4), we maximize the marginal log-likelihood function of $\bm{\xi}_{k}=[\xi_{1k},...,\xi_{Nk}]^{\prime}$ given $\bm{X}^{(-s)}_{1}(t),...,\bm{X}^{(-s)}_{N}(t)$. Let $\mathcal{X}$ denote the observations of signals $X^{(l)}_{i}(t)$ for $l\in L^{-s}$ units $i\in I^{\mathcal{H}}$, that is \mathcal{X}=\big{\{}X^{(l)}_{i}(t)|{t\in T_{i}^{(l)},l\in L^{-s},i\in I^{\mathcal{H}}}\big{\}} where $T_{i}^{(l)}=\{t_{iu}|i\in I^{\mathcal{H}},l\in L^{(-s)},u=1,...,p^{(l)}_{i}\}$. Also, let $z_{ik}$ denote the true underlying latent values corresponding to the FPC scores $\xi_{ik}$ and let $\bm{z}_{k}=[z_{1k},...,z_{Nk}]^{\prime}$. Then the marginal likelihood is given as

[TABLE]

where $P(\xi_{ik}|z_{ik},\mathcal{X},\bm{\Theta}_{k})=\mathcal{N}(0,\sigma_{k}^{2})$ and P(\bm{z}_{k}|\mathcal{X},\bm{\Theta}_{k})=\mathcal{N}\big{(}\bm{0}_{N},\mathbf{C}_{k}\big{)}. The second equality follows from the fact that the error is an additive Gaussian noise. Thus, $\prod_{i\in I^{\mathcal{H}}}P(\xi_{ik}|z_{ik},\mathcal{X},\bm{\Theta}_{k})=\mathcal{N}(\bm{0}_{N},\sigma_{k}^{2}\mathbf{I}_{N})$ and the log-likelihood of $P(\bm{\xi}_{k}|\mathcal{X},\bm{\Theta}_{k})$ is

[TABLE]

where $\langle\bm{A},\bm{B}\rangle_{tr}=trace(\bm{A}\bm{B})$ and $\bm{\Xi}_{k}=\bm{\xi}_{k}\bm{\xi}^{\prime}_{k}$. As a consequence, the optimized hyperparameters denoted by $\bm{\Theta}_{k}^{*}=[\bm{\theta}^{*^{\prime}}_{k},\sigma_{k}^{*}]^{\prime}$ are found by maximizing the marginal log-likelihood. More formally, we have

[TABLE]

Following multivariate normal theory, the posterior predictive distribution of $\xi_{rk}$, given (4) and $\bm{\Theta}_{k}^{*}$, is derived as

[TABLE]

with

[TABLE]

Here we note that for each $k\in\{1,...,K\}$, we can derive $\hat{\xi}_{rk}$ and $\hat{\sigma}_{rk}^{2}$ using an independent GP as the FPC scores from different orthonormal basis functions are uncorrelated. This facilitates scalability of computation as for different $k$ we can derive $\{\hat{\xi}_{rk}\}$ and $\{\hat{\sigma}_{rk}^{2}\}$ in parallel. As shown in the computational complexity derivations in section 6, this aspect is important specifically in a real-time environment where predictions need to be continuously updated.

Now combining equations (2) and (6), we obtain the predictive mean $\hat{X}_{r}(t)$ and variance $\hat{\sigma}_{r}^{2}(t)$ of ${X}_{r}(t)$ as follows

[TABLE]

Here note that $\mu(t)$, $\sigma^{2}_{\mu}(t)$ and $\sigma^{2}_{\epsilon}$ are model parameters corresponding to the estimated FPCA model in (2), where $\sigma^{2}_{\epsilon}$ denotes the estimated variance of $\epsilon(t)$. Here we recall that $\hat{X}_{r}(t)=\hat{X}^{(s)}_{r}(t)$ as the index $(s)$ is dropped for the target stream.

4.3 Empirical Bayesian Updating with On-line Data

In the previous section, we derive a prior for $\hat{X}_{r}(t)$ and $\hat{\sigma}_{r}^{2}(t)$ using data observed from streams $l\in L^{-s}$. Here, we develop an empirical Bayesian approach to update $\hat{X}_{r}(t)$ and $\hat{\sigma}_{r}^{2}(t)$ given the target stream ($s$) observations from the in-service unit $r$. Specifically, given the prior distributions P(\xi_{rk})=\mathcal{N}\big{(}\hat{\xi}_{rk},\hat{\sigma}^{2}_{rk}\big{)} for each $k$ and given the observations $\bm{X}_{r}(\bm{t})$ at $\bm{t}=[t_{r1},...,t_{rp_{r}}]^{\prime}$, the posterior $P(\xi_{rk}\>|\>\bm{X}_{r}(\bm{t}))$ is given in Proposition 2.

Proposition 2**.**

Given that $X_{r}(t)=\mu(t)+\sum_{k=1}^{K}\xi_{rk}\phi_{k}(t)+\epsilon(t)$, where the prior distribution of $\xi_{rk}$ is \mathcal{N}\big{(}\hat{\xi}_{rk},\hat{\sigma}_{rk}^{2}\big{)}, and the FPC scores are pairwise independent. Then, the posterior distribution of the FPC scores, such that $\xi^{*}_{rk}=P(\xi_{rk}\>|\>\bm{X}_{r}(\bm{t}))$, is given as

[TABLE]

where

[TABLE]

with

[TABLE]

[TABLE]

Proof.

See Section A.2 in supplementary document. ∎

Based on the updated FPC scores for in-service unir $r$, the posterior predicted mean $\hat{X}^{*}_{r}(\tilde{t})$, of $X_{r}(t)$ for any future time point $\tilde{t}\geq t^{\ast}$ where $\tilde{t}\in\mathcal{T}$ is given as

[TABLE]

Similarly, the posterior variance \big{(}\hat{\sigma}^{*}_{r}(\tilde{t})\big{)}^{2} can be computed as

[TABLE]

where \big{(}\bm{\Sigma}^{*}\big{)}_{k_{1}k_{2}} indicates the $(k_{1},k_{2})$th element of the covariance matrix $\bm{\Sigma}^{*}$.

Despite our focus on the target stream $s$ we note that our framework can predict every individual stream for the in-service unit $r$. This ability to provide individualized predictions is a key feature of the proposed methodology compared the data fusion literature that predicts a single aggregated signal. Further, one differentiating factor is that we allow irregularly sampled data from each stream where time points of each signal $\{t^{(l)}_{iu}:u=1,...,p^{(l)}_{i}\}\subset\mathcal{T}$ do not need to be identical or regularly spaced across streams. Indeed, such situations are quite common in practice because most multi-stream data is gathered from different types of sensors. Therefore, the proposed approach is applicable to a wide array of practical situations.

5 Numerical Case Study

5.1 General Settings

In this section, we discuss the general settings used to assess the proposed model, denoted as FPCA-GP. We evaluate the model by performing experiments with both synthetic and real-world data. We report the prediction accuracy at varying time points $t^{\ast}$ for the partially observed unit $r$. Specifically, for the time domain $\mathcal{T}=[a,b]$, we assume that the on-line signals from the in-service unit $r$ are partially observed in the range of $[a,t^{\ast}=a+\gamma(b-a)]$, referred as $\gamma$-observation. We set $\gamma=25\%,50\%$, and $75\%$ for every case study. In the extrapolation interval $[t^{\ast},b]$, we use the mean absolute error (MAE) between the true signal value $X^{T}_{r}(t)$ and its predicted value $\hat{X}^{(s)}_{r}(t)$ at $U$ evenly spaced test points (denoted as $t_{u}$ for $u=1,...,U$) as the criterion to evaluate our prediction accuracy.

[TABLE]

We report the distribution of the errors across $G$ repetitions using a group of boxplots representing the MAE for the testing unit $r$ at diffrent $\gamma$-observation percentiles. Further, we benchmark our method with two other reference methods for comparison: (i) The FPCA approach for single stream settings denoted as FPCA-B. In this method we only consider the target stream $s$ (Zhou et al., 2011; Kontar et al., 2018). Note that we incorporate our Bayesain updating procedure to update predictions as new data is observed. (ii) The Bayesian mixed effect model with a general polynomial function whose degree is determined through an Akaike information criteria (AIC) (Rizopoulos, 2011; Son et al., 2013; Kontar et al., 2017). We denote this methods as ME. The ME model intrinsically applies a Bayesian updating scheme as more data is obtained from the in-service unit. Detailed codes for both reference methods are included in our supplementary materials.

5.2 Numerical Study with Synthetic Data

First, we show the numerical results of the proposed model performed on synthetic data. For this experiment, we assume that two streams $l\in\{1,2\}$ of data are observed from two different sensors embedded in each unit. The target stream of interest is $l=1$. To generate signals possessing heterogeneity, we suppose there are two separating environments, denoted by environment I and II. We generate signals for each unit using different underlying functions depending on which environment the unit is in. This is illustrated in Figure 2. As shown in the figure, the underlying trend of the target stream ($l=1$) will vary under different profiles of stream $l=2$. To relate this setting with real-world application, consider $l=1$ as the degradation level and $l=2$ as the temperature profile. Then from Figure 2, we have that units operating under different temperature profiles will exhibit different trends.

We generate a training set of $N=50$ units and one testing unit $r$ whose signals are partially observed. Also we repeat the experiment $G=100$ times. Historical units are operated in either environment I or II whereas the in-service unit is operated in environment II. The population of historical units is created according to three levels of heterogeneity: (i) 0% heterogeneity where all units in the historical database are operated under environment II (similar to that of the testing unit) (ii) 50% heterogeneity where 25 units are distributed to each environment (iii) 90% heterogeneity where only 5 units are assigned to environment II. Conducting the experiments across a homogeneous setting and a heterogeneous setting, where the in-service unit belongs to the minority group with only 10% ratio, will allow us to investigate the robustness of our approach.

For units in environment I, the signals from respective streams $l=1,2$ are generated according to $X_{i}^{(1)}(t)=0.3t^{2}-2\sin{(w_{1,i}\pi t)}+w^{\text{I}}_{2,i}$ and $X_{i}^{(2)}(t)=2w_{1,i}\sin t$, where $w_{1,i}\sim\mathcal{N}(0.4,0.03^{2})$ and $w^{\text{I}}_{2,i}\sim\mathcal{U}(0,5)$, where $\mathcal{U}$ denotes the uniform distribution. For units in environment II, we generate the signals as $X_{i}^{(1)}(t)=0.3t^{2}-2\sin{(w_{1,i}\pi t^{0.85})}+3(\arctan(t-5)+\frac{\pi}{2})+w^{\text{II}}_{2,i}$ and $X_{i}^{(2)}(t)=2w_{1,i}\sin 0.3t$ where $w^{\text{II}}_{2,i}\sim\mathcal{U}(1.5,6.5)$. Measurement error $\epsilon\sim\mathcal{N}(0,0.05^{2})$ is assumed similar across both streams.

Figure 2 illustrates training signals $X^{(1)}_{i}(t)$ in the case of 90% heterogeneity. It is crucial to note that at early stages (ex: $\gamma=25\%$), it is hard to distinguish between the two different trends in stream $l=1$. We model that on purpose to check if our model can leverage information from stream $l=2$ to uncover the underlying heterogeneity at early stages. This in fact is a common feature in many health related applications, as many diseases remain dormant at early stages and it is only through measuring other factors we can predict there evolution early on.

The results are illustrated in both Figure 3 and 4. Based on the figures we can obtain some important insights. First, the FPCA-GP clearly outperforms the FPCA-B. This is specifically obvious at early stages ($\gamma=25\%$) and when the data exhibits heterogeneity (90% and 50% heterogeneity). This confirms the ability of our model to borrow strength from information across different streams to discern the heterogeneity and enhance predictive accuracy at early stages. This result is very motivating specifically since at $\gamma=25\%$ data from in-service unit $r$ is sparse and all signals in stream $s$ have similar behaviour which makes it hard to uncover future heterogeneity. It further implies that the FPC scores of the testing unit are appropriately estimated by the proposed approach, as shown in the first column of Figure 3. From the figure, we observe that the estimated prior mean from the FPCA-GP appropriately follows the signals in environment II, whereas the prior mean from the FPCA-B follows the signals in environment I, which is the majority. Second, as expected, prediction errors significantly decrease as the percentiles increase. Thus, our Bayesian updating framework is able to efficiently utilize new collected data and provide more accurate predictions as $t^{\ast}$ increases. Third, the results show that ME behaved the worst and its predictions accuracy merely decreases at later stages. This result illustrates the vulnerabilities of parametric modeling and demonstrates the ability of our non-parametric modeling to avoid model misspecifications. Fourth, the results confirm that even in the case where other streams have no effect on the target stream (0% heterogeneity) the FPCA-GP is competitive compared to FPCA-B. This highlights the robustness of the FPCA-GP.

5.3 Numerical Study with Real-world Data

In this section, we discuss the numerical study using real-world data provided by the National Aeronautics and Space Administration (NASA). The dataset contains degradation signals collected from multiple sensors on an aircraft turbofan engine. This dataset was generated from a simulation model, developed in Matlab Simulink, called commercial modular aero-propulsion system simulation (C-MAPSS). This system simulates degradation signals from multiple-sensors, installed in several components of an aero turbofan engine, under a variety of environmental conditions. The list of the components includes Fan, LPC, HPC, and LPT, and are illustrated in Figure 5. Refer to Saxena & Simon (2008) for more details about turbofan engine data. The dataset is available at Saxena & Goebel (2008). The dataset is composed of 21 sensor streams from 100 training and 100 testing units. Following the analysis of Liu et al. (2013), we select the 11 most crucial streams. Some signals from these streams are shown in Figure 6. We provide the detailed list and description of sensors in the supplementary materials. In our analysis we truncate the time range $(0,100]$ and predict the testing signal over the time range $(100,160]$.

Table 1 demonstrates that the MAE results of stream 4 and 15. Note that these two streams have shown to have the largest impact on failure (Fang et al., 2017b) and therefore, due to space limitation, we only focus on their predictive results. Results for other signals are provided in the supplementary materials. Note that we include the standard deviation of MAE across the testing units.

The results clearly show that our approach is far more superior than benchmarks for the real-world data. For all provided cases, the mean of MAE for the FPCA-GP is less than that of the FPCA. Once again this highlights the importance of leveraging information from all streams of the data. Another important insight from this study is that our model was able to outperform the ME even though the curves from Figure 6 seem to exhibit a clear parametric trend. This further highlights the robustness of our method and its ability to safeguard against parametric misspecifications.

6 Discussion

In this study, we developed a non-parametric statistical model that can extrapolate individual signals in a multi-stream data setting. Using both synthetic and real-world data, we demonstrate our models ability to borrow strength across all streams of data, predict individual streams, account for heterogeneity and provide accurate real-time predictions where an empirical Bayesian approach updates our predictor as new data is observed in real-time. Since we work in the regime of streaming data, the frequency with which we receive data is very high. Due to this, our model needs to be efficient in terms of the time taken to make each update. With the assumption that all signals from $M$ streams have $Q$ observations, the complexity of multivariate FPCA for multi-stream data is $\mathcal{O}(M^{2}NQ^{2}+M^{3}Q^{3})$ (Fang et al., 2017b). In our model, the computationally expensive steps are the FPCA for the target stream (Section 4.1) and the implementation of GP for estimating the FPC scores (Section 4.2). Following Xiao et al. (2016), the complexity of the former is $\mathcal{O}(QN^{2}+N^{3})$. While complexity of a GP with an $N\times N$ covariance matrix is $\mathcal{O}(N^{3})$ (Rasmussen & Williams, 2005). Given that we implement $K$ independent GP models the complexity of estimating the FPC scores is $\mathcal{O}(KN^{3})$. Combing the above observations, we conclude that the complexity of our procedure is $\mathcal{O}(QN^{2}+N^{3}+KN^{3})$. Typically, we have that $M,N,K\ll Q$, also, in real-time $Q$ is increasing. Thus, our model is clearly more efficient than multivariate the FPCA and applicable in practice in a real-time streaming environment.

7 Software and Data

Technical proofs, the used dataset, a detailed code and additional numerical results are available in the supplementary materials.