# A Supervised Tensor Dimension Reduction-Based Prognostic Model for
Applications with Incomplete Imaging Data

Chengyu Zhou
Edward P. Fitts Department of Industrial and Systems Engineering,
North Carolina State University
and
Xiaolei Fang
Edward P. Fitts Department of Industrial and Systems Engineering,
North Carolina State University

###### Abstract

Imaging data-based prognostic models focus on using an asset’s degradation
images to predict its time-to-failure (TTF). Most image-based prognostic
models have two common limitations. First, they require degradation images to
be complete (i.e., images are observed continuously and regularly over time).
Second, they usually employ an unsupervised dimension reduction method to
extract low-dimensional features, and then use the features for TTF
prediction. Since unsupervised dimension reduction is conducted on the
degradation images without the involvement of TTFs, there is no guarantee that
the extracted features are effective for failure time prediction. To address
these challenges, this article develops a supervised tensor dimension
reduction-based prognostic model. The model first proposes a supervised
dimension reduction method for tensor data. It uses historical TTFs to guide
the detection of a tensor subspace to extract low-dimensional features from
high-dimensional incomplete degradation imaging data. Next, the extracted
features are used to construct a prognostic model based on (log)-location-
scale regression. An optimization algorithm is proposed for parameter
estimation, and closed-form solutions are derived under certain distributions.
Simulated data and a real-world data set are used to validate the performance
of the proposed model.

Keywords: Supervised dimension reduction, missing data, tensor, failure time
prediction

##  1 Introduction

Degradation is an irreversible process of damage accumulation that results in
the failure of engineering systems/assets/components (Bogdanoff and Kozin,
1985). Although it is usually challenging to observe a physical degradation
process, there often are some manifestations associated with degradation
processes that can be monitored by sensing technology, which yields data known
as degradation data/signals. Degradation data contains the health condition of
engineering assets; thus, if modeled properly, they can be used to predict the
assets’ time-to-failure (TTF) via a process known as prognostic. Many
prognostic models have been developed in the literature, most of which focus
on using time series-based degradation data (Ye et al., 2014; Ye and Chen,
2014; Hong and Meeker, 2010, 2013; Shu et al., 2015; Liu et al., 2013;
Gebraeel et al., 2005). Recently, prognostic models with imaging-based
degradation data have been investigated and attracted more and more attention
(Fang et al., 2019; Aydemir and Paynabar, 2019; Yang et al., 2021; Dong et
al., 2021; Tang et al., 2021; Jiang et al., 2022). This is because comparing
with time-series data, imaging data usually contains much richer information
of the object being monitored, and imaging sensing technologies are noncontact
and thus they can usually be easily deployed. One example of imaging-based
degradation data is the infrared image stream that measures the change of
temperature distribution of a thrust bearing during its degradation process
over time (Fang et al., 2019; Aydemir and Paynabar, 2019; Dong et al., 2021;
Jiang et al., 2022). Another example is the images used to measure the
performance degradation of infrared systems such as rotary-wing drones (Dong
et al., 2021).

The existing imaging-based prognostic methods include deep learning-based
models and statistical learning methods. Examples of the deep learning-based
models designed for TTF prediction using imaging data include the ones
developed by Aydemir and Paynabar (2019); Yang et al. (2021); Dong et al.
(2021), and Jiang et al. (2022). Although these models have worked relatively
well, they usually provide point estimations of failure times, and it is
challenging for them to quantify the uncertainty of predicted TTFs. This
limits their applicability since the subsequent decision-making analysis such
as maintenance/inventory/logistic optimization requires prognostic models to
provide a distribution of the predicted TTF. One example of statistical
learning methods for image-based prognostic model is the penalized
(log)-location-scale (LLS) tensor regression proposed by Fang et al. (2019).
The model first employs multilinear principal component analysis (MPCA) (Lu et
al., 2008) to reduce the dimension of high-dimensional imaging-based
degradation data, which yields a low-dimensional feature tensor for each
asset. Next, it constructs a prognostic model by regressing an asset’s TTF
against it low-dimensional feature tensor using LLS regression. In the same
article, Fang et al. (2019) have also proposed several benchmarking prognostic
models that use imaging-based degradation data for TTF prediction. These
models also first employ a dimension reduction method such as functional
principal component analysis (FPCA) (Ramsay and Silverman, 2005), principal
component analysis (PCA) (Abdi and Williams, 2010), or B-Spline (Prautzsch et
al., 2002), to reduce the dimension of degradation data and then use low-
dimensional features to build an LLS regression model for prognostic. Although
the aforementioned statistical learning-based prognostic models can provide a
distribution for the predicted TTF, and their effectiveness have been well
investigated, they share two common limitations.

The first limitation is that they assume imaging-based degradation data
(including historical data for model training and real-time data for model
test) are complete, which means images from all the assets should be collected
continuously and regularly with the same sampling time interval (see Figure 1
(a) for an example). In reality, however, engineering assets often operate in
harsh environments that significantly impact the quality of collected data due
to errors in data acquisition, communication, read/write operations, etc. As a
result, degradation images often contain significant levels of missing
observations, which is known as incomplete/missing imaging data (see Figure 1
(b) for an example). Such data incompleteness poses a significant challenge
for the parameter estimation of existing statistical learning-based prognostic
models.

(a) Complete Image Streams

(b) Incomplete Image Streams

Figure 1:  Degradation stream images with and without missing data.

The second common limitation for the existing statistical learning-based
prognostic models for applications with imaging data is that they employ
unsupervised dimension reduction methods for feature extraction, so there is
no guarantee that the extracted features are effective for the subsequent TTF
prediction. Specifically, they first use unsupervised dimension reduction
methods such as FPCA, PCA, and B Spline to extract features, which are then
used to construct prognostic models. Since feature extraction and prognostic
model construction are two sequential steps, and no TTF information gets
involved in the feature extraction process, it is possible that the extracted
features may not be most suitable for predicting TTFs.

To address the aforementioned challenges, this article proposes a supervised
dimension reduction-based prognostic model that uses an asset’s incomplete
degradation images to predict its TTF. Similar to the existing statistical
learning-based prognostic models, the proposed model also consists of two
steps: feature extraction and prognostic model construction. However, unlike
the existing models, feature extraction in this article is achieved by
developing a new supervised tensor dimension reduction method, which uses
historical TTFs to supervise the feature extraction process such that the
extracted features are more effective for the subsequent TTF prediction. In
addition, unlike the existing unsupervised dimension reduction methods that
only work for complete imaging data, the proposed supervised dimension
reduction method works for both complete and incomplete degradation image
streams.

The proposed supervised dimension reduction method works as follows: First, it
detects a low-dimensional tensor subspace in which the high-dimensional
degradation image streams are embedded. This is achieved by constructing an
optimization criterion that comprises a feature extraction term and a
regression term. The first term extracts low-dimensional features from
complete/incomplete degradation image streams of training assets, and the
second term builds the connection between these assets’ TTFs and the extracted
features using LLS regression. LLS regression has been widely used in
reliability engineering and survival analysis. It includes a variety of TTF
distributions, such as (log)normal, (log)logistic, smallest extreme value, and
Weibull, etc., which cover most of the TTF distributions in reality (Doray,
1994). Since historical TTFs are used to supervise the feature extraction
process, it is expected that the extracted features are more effective for the
subsequent prognostic. Solving the optimization criterion of the proposed
supervised tensor dimension reduction method yields a set of tensor basis
matrices that span the low-dimensional tensor subspace for dimension
reduction. We then expand both the historical degradation images in the
training data set and real-time degradation images from an asset operating in
the field (i.e., test data) using the set of tensor basis matrices to extract
the low-dimensional tensor features of the training assets and the test asset.
The TTFs of the training assets are then regressed against their tensor
features using LLS regression, and the parameters are estimated using maximum
likelihood estimation. After that, the tensor features of the test asset are
fed into the LLS regression model, and its TTF distribution is predicted.

To solve the optimization criterion of the proposed supervised dimension
reduction method, we will first transfer the criterion into a block
multiconvex problem. Next, we will propose a block updating algorithm, which
cyclically optimizes one block parameters while keeping other blocks fixed
until convergence. In addition, we will demonstrate that when TTFs follow
normal or lognormal distributions, each sub problem of the block updating
algorithm has a closed-form solution, no matter the degradation image streams
are complete or incomplete.

The rest of this paper is orginized as follows. Section 2 presents the
supervised tensor dimension reduction-based prognostic method. Section 3
introduces the block updating algorithm and closed-form solutions when TTFs
follow normal/lognormal distributions. Sections 5 and 6 validate the
effectiveness of the proposed prognostic model using a simulated dataset and
data from a rotating machinery, respectively. Section 7 concludes.

##  2 The Methodology

In this section, we will introduce the proposed supervised tensor dimension
reduction-based prognostic model for applications with incomplete imaging
data. In Subsection 2.1, we will present some basic tensor notations and
definitions. Subsection 2.2 introduces the supervised tensor dimension
reduction method. In Subsection 2.3, we will discuss the construction of
prognostic model and how to predict the TTF of an asset operating in the field
using its real-time degradation imaging data.

###  2.1 Preliminaries

In this section, we introduce some basic notations and definitions of tensor
operations that are used throughout the article. The order of a tensor is the
number of dimensions, also known as ways or modes. Vectors (1 order tensors)
are denoted by lowercase boldface letters, e.g., $\bm{s}$. Matrices (2 order
tensors) are denoted by boldface uppercase letters, e.g., $\bm{S}$. Higher-
order tensors (order is 3 or larger) are denoted by calligraphic letters,
e.g., $\mathcal{S}$. Indices are denoted by lowercase letters whose range is
from 1 to the uppercase letter of the index, e.g., $n=1,2,\ldots,N$. An $N$th-
order tensor is denoted as $\mathcal{S}\in\mathbb{R}^{I_{1}\times
I_{2}\times\cdot\cdot\cdot\times I_{N}}$, where $I_{n}$ represents the $n$th
mode of $\mathcal{S}$. The $(i_{1},i_{2},\ldots,i_{N})$th entry of
$\mathcal{S}\in\mathbb{R}^{I_{1}\times I_{2}\times\cdot\cdot\cdot\times
I_{N}}$ is denoted by $s_{i_{1},i_{2},\ldots,i_{n}}$. A fiber of $\mathcal{S}$
is a vector defined by fixing every index but one. A matrix column is a mode-1
fiber and a matrix row is a mode-2 fiber. The vectorization of $\mathcal{S}$,
denoted by vec($\mathcal{S}$), stacks all the entries of $\mathcal{S}$ into a
column vector. The mode-$n$ matricization of a tensor
$\mathcal{S}\in\mathbb{R}^{I_{1}\times I_{2}\times\cdot\cdot\cdot\times
I_{N}}$ is denoted by $\bm{S}_{(n)}$, which arranges the mode-$n$ fibers to be
the columns of the resulting matrix. The $n$th mode product of a tensor
$\mathcal{S}\in\mathbb{R}^{I_{1}\times I_{2}\times\cdot\cdot\cdot\times
I_{N}}$ and a matrix $\bm{U}_{n}\in\mathbb{R}^{J_{n}\times I_{n}}$, denoted by
$\mathcal{S}\times_{n}\bm{U}_{n}$, is a tensor whose entry is
$(\mathcal{S}\times_{n}\bm{U}_{n})_{i_{1},\ldots,i_{n-1},j_{n},i_{n+1},\ldots,i_{N}}=\Sigma_{I_{n}=1}^{I_{N}}s_{i_{1},\ldots,i_{N}}u_{j,i_{n}}$.
The Kronecker product of two matrices $\bm{A}\in\mathbb{R}^{m\times n}$ and
$\bm{B}\in\mathbb{R}^{p\times q}$ is an $mp\times nq$ block matrix:

$\bm{A}\otimes\bm{B}=\begin{bmatrix}\bm{a}_{11}\bm{B}&\ldots&\bm{a}_{1n}\bm{B}\\\
\vdots&\ddots&\vdots\\\
\bm{a}_{m1}\bm{B}&\ldots&\bm{a}_{mn}\bm{B}\end{bmatrix}.$

If $\bm{A}$ and $\bm{B}$ have the same number of columns $n=q$, then the
Khatri-Rao product is defined as the $mp\times n$ column-wise Kronecker
product:
$\bm{A}\odot\bm{B}=[\bm{a}_{1}\otimes\bm{b}_{1}\quad\bm{a}_{2}\otimes\bm{b}_{2}\quad\cdot\cdot\cdot\quad\bm{a}_{n}\otimes\bm{b}_{n}]$.
If $\bm{a}$ and $\bm{b}$ are vectors, then
$\bm{A}\otimes\bm{B}=\bm{A}\odot\bm{B}$. More details about tensor notations
and operators can be found in Kolda and Bader (2009).

###  2.2 The Supervised Tensor Dimension Reduction Method

We assume that there exists a historical data set for model training. The data
set consists of the degradation image streams of $M$ failed assets along with
their TTFs, which are denoted as $\mathcal{X}_{m}\in\mathbb{R}^{I_{1}\times
I_{2}\times I_{3}}$ and $y_{m}\in\mathbb{R}$, respectively, where
$m=1,2,\ldots,M$. For the convenience of introducing the dimension reduction
method, we convert the 3D degradation image streams from all the $M$ assets to
a 4D tensor $\mathcal{X}\in\mathbb{R}^{I_{1}\times I_{2}\times I_{3}\times
M}$, where the sample size $M$ is the 4th mode. Similarly, we let
$\bm{y}=(y_{1},\ldots,y_{M})^{\top}\in\mathbb{R}^{M\times 1}$ be the vector
containing all the TTFs of the $M$ assets.

Out of the $I_{1}\times I_{2}\times I_{3}\times M$ entries of $\mathcal{X}$,
we use a subset $\Omega\subseteq\\{(i_{1},i_{2},i_{3},m),1\leq i_{1}\leq
I_{1},1\leq i_{2}\leq I_{2},1\leq i_{3}\leq I_{3},1\leq m\leq M\\}$ to denote
the indices of the missing ones. To model the missing data, we define a
projection operator $\mathcal{P}_{\Omega}(\cdot)$ as follows:

$\mathcal{P}_{\Omega}(\mathcal{X})_{(i_{1},i_{2},i_{3},m)}=\begin{cases}\mathcal{X}_{(i_{1},i_{2},i_{3},m)},&\text{if}~{}(i_{1},i_{2},i_{3},m)\not\in\Omega,\\\
0,&\text{if}~{}(i_{1},i_{2},i_{3},m)\in\Omega,\end{cases}$ (1)

where $\mathcal{X}_{(i_{1},i_{2},i_{3},m)}$ is the $(i_{1},i_{2},i_{3},m)$-th
entry of the 4D tensor $\mathcal{X}\in\mathbb{R}^{I_{1}\times I_{2}\times
I_{3}\times M}$. To recover the missing entries in tensor $\mathcal{X}$, we
may employ the following Tucker decomposition-based tensor completion method
(Filipović and Jukić, 2015; Liu et al., 2012; Xu et al., 2013):

$\min_{\mathcal{S},\bm{U}_{1},\bm{U}_{2},\bm{U}_{3}}~{}~{}||\mathcal{P}_{\Omega}(\mathcal{X}-\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{2}\bm{U}_{2}^{\top}\times_{3}\bm{U}_{3}^{\top})||_{F}^{2}.$
(2)

where $||\cdot||_{F}^{2}$ is the Frobenius norm,
$\bm{U}_{1}\in\mathbb{R}^{P_{1}\times I_{1}}$,
$\bm{U}_{2}\in\mathbb{R}^{P_{2}\times I_{2}}$,
$\bm{U}_{3}\in\mathbb{R}^{P_{3}\times I_{3}}$ are three factor matrices,
$\mathcal{S}\in\mathbb{R}^{P_{1}\times P_{2}\times P_{3}\times M}$ is the low-
dimensional core tensor, and $\times_{n}$ is the n-mode product of a tensor
with a matrix. The tensor completion criterion (2) can be seen as an
unsupervised dimension reduction method for tensor data with missing entries.
This is because the degradation image tensor
$\mathcal{X}\in\mathbb{R}^{I_{1}\times I_{2}\times I_{3}\times M}$ is a
4-order tensor, which resides in the tensor (multilinear) space
$\mathbb{R}^{I_{1}}\otimes\mathbb{R}^{I_{2}}\otimes\mathbb{R}^{I_{3}}\otimes\mathbb{R}^{M}$,
where
$\mathbb{R}^{I_{1}},\mathbb{R}^{I_{2}},\mathbb{R}^{I_{3}},\mathbb{R}^{M}$ are
the 4 vector (linear) spaces; $\mathcal{S}\in\mathbb{R}^{P_{1}\times
P_{2}\times P_{3}\times M}$ can be seen as a feature tensor that resides in
the tensor space
$\mathbb{R}^{P_{1}}\otimes\mathbb{R}^{P_{2}}\otimes\mathbb{R}^{P_{3}}\otimes\mathbb{R}^{M}$.
Usually, we have $P_{1}\ll I_{1}$, $P_{2}\ll I_{2}$, and $P_{3}\ll I_{3}$ for
degradation imaging data due to the high spatio-temporal correlation among the
pixels. This implies that the dimension of the image stream from the $m$th
asset is reduced from $\mathbb{R}^{I_{1}\times I_{2}\times I_{3}}$ to
$\mathbb{R}^{P_{1}\times P_{2}\times P_{3}}$, where $m=1,\ldots,M$.

Although criterion (2) can be seen as a dimension reduction method, there is
no guarantee that the extracted low-dimensional feature tensor $\mathcal{S}$
is effective for the subsequent TTF prediction. To address this challenge, we
propose the following supervised dimension reduction method by combining a
tensor completion term and an LLS regression term:

$\min_{\bm{U}_{1},\bm{U}_{2},\bm{U}_{3},\bm{\beta}_{1},\beta_{0},\mathcal{S}}\alpha||\mathcal{P}_{\Omega}(\mathcal{X}-\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{2}\bm{U}_{2}^{\top}\times_{3}\bm{U}_{3}^{\top})||_{F}^{2}+(1-\alpha)\ell(\frac{\bm{y}-\bm{1}_{M}\beta_{0}-\bm{S}_{(4)}\bm{\beta}_{1}}{\sigma}),$
(3)

where $\bm{y}\in\mathbb{R}^{M\times 1}$ is the vector containing all the TTFs
of the $M$ assets in the training data set. The matrix
$\bm{S}_{(4)}\in\mathbb{R}^{M\times(P_{1}\times P_{2}\times P_{3})}$ is the
mode-4 matricization of the low-dimensional feature tensor $\mathcal{S}$, the
$m$th row of which represents the vectorization of the $m$th asset’s feature
tensor, $m=1,\ldots,M$. $\beta_{0}$ is the intercept, and
$\bm{\beta_{1}}\in\mathbb{R}^{(P_{1}\times P_{2}\times P_{3})\times 1}$ is the
regression coefficient vector. $\bm{1}_{m}\in\mathbb{R}^{M\times 1}$ is an
$M\times 1$ vector whose elements are all ones. $\ell(\cdot)$ is the negative
log-likelihood function of a location-scale distribution. For example, if TTFs
follow normal distributions, then
$\ell(\frac{\bm{y}-\bm{1}_{M}\beta_{0}-\bm{S}_{(4)}\bm{\beta}_{1}}{\sigma})=\frac{M}{2}\log
2\pi+M\log\sigma+\frac{1}{2}\Sigma_{m=1}^{M}\omega_{m}^{2}$, where
$\omega_{m}=\frac{y_{m}-\beta_{0}-\bm{s}_{(4)}^{m}}{\sigma}$, where
$\bm{s}_{(4)}^{m}$ is the $m$th row of $\bm{S}_{(4)}$, and $y_{m}$ is the TTF
of asset $m$; if TTFs follow logistic distributions, then
$\ell(\frac{\bm{y}-\bm{1}_{M}\beta_{0}-\bm{S}_{(4)}\bm{\beta}_{1}}{\sigma})=M\log\sigma-\Sigma_{m=1}^{M}\omega_{m}+2\Sigma_{m=1}^{M}\log(1+\exp(\omega_{m}))$;
if TTFs follow small extreme value (SEV) distributions, then
$\ell(\frac{\bm{y}-\bm{1}_{M}\beta_{0}-\bm{S}_{(4)}\bm{\beta}_{1}}{\sigma})=n\log\sigma-\Sigma_{m=1}^{M}\omega_{m}+\Sigma_{m=1}^{M}\exp(\omega_{m})$.
For assets whose TTFs follow log-location-scale distributions, we may transfer
them to the corresponding location-scale distributions by taking their
logarithm such that criterion (3) can still be used. For example, log-normal,
log-logistics, and Weibull distributions can be transferred to normal,
logistics, and SEV distributions, respectively. $\alpha\in[0,1]$ is a weight
and $||\cdot||_{F}^{2}$ is the Frobenius norm.

In criterion (3), the first term
$||\mathcal{P}_{\Omega}(\mathcal{X}-\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{2}\bm{U}_{2}^{\top}\times_{3}\bm{U}_{3}^{\top})||_{F}^{2}$
is tensor completion from (2), which reduces the dimension of high-dimensional
incomplete degradation image streams and extracts low-dimensional tensor
features. The second term
$\ell(\frac{\bm{y}-\bm{1}_{M}\beta_{0}-\bm{S}_{(4)}\bm{\beta}_{1}}{\sigma})$
is LLS regression, which regresses each asset’s TTF against its tensor
features extracted by the first term. By jointly optimizing the two terms, it
is expected that the extracted features are effective for TTF prediction.
However, it is challenging to solve criterion (3) since it is neither convex
nor block multi-convex. An optimization problem is block multi-convex when its
feasible set and objective function are generally non-convex but convex in
each block of variables (Xu and Yin, 2013). Thus, to simplify the development
of optimization algorithms for model parameter estimation, we first transform
criterion (3) to be a block multi-convex one. Specifically, we apply the
following re-parameterization:
$\tilde{\sigma}=1/\sigma,\tilde{\beta_{0}}=\beta_{0}/\sigma,\tilde{\bm{\beta}_{1}}=\bm{\beta}_{1}/\sigma$.
As a result, criterion (3) can be re-expressed as follows:

$\min_{\bm{U}_{1},\bm{U}_{2},\bm{U}_{3},\tilde{\bm{\beta_{1}}},\tilde{\beta_{0}},\bm{S}_{(4)}}\alpha||\mathcal{P}_{\Omega}(\mathcal{X}-\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{2}\bm{U}_{2}^{\top}\times_{3}\bm{U}_{3}^{\top})||_{F}^{2}+(1-\alpha)\ell(\tilde{\sigma}\bm{y}-\bm{1}_{M}\tilde{\beta_{0}}-\bm{S}_{(4)}\tilde{\bm{\beta}_{1}}),$
(4)

where,
$\ell(\tilde{\sigma}\bm{y}-\bm{1}_{M}\tilde{\beta_{0}}-\bm{S}_{(4)}\tilde{\bm{\beta}_{1}})=\frac{M}{2}\log
2\pi-M\log\tilde{\sigma}+\frac{1}{2}\Sigma_{m=1}^{M}\tilde{\omega}_{m}^{2}$
for TTFs following normal distributions, and
$\tilde{\omega}_{m}=\tilde{\sigma}y_{m}-\tilde{\beta}_{0}-{\bm{s}}_{(4)}^{m}$,
where $\bm{s}_{(4)}^{m}$ is the $m$th row of $\bm{S}_{(4)}$ and $y_{m}$ is the
TTF of asset $m$;
$\ell(\tilde{\sigma}\bm{y}-\bm{1}_{M}\tilde{\beta_{0}}-\bm{S}_{(4)}\tilde{\bm{\beta}_{1}})=-M\log\tilde{\sigma}-\Sigma_{m=1}^{M}\tilde{\omega}_{m}+2\Sigma_{m=1}^{M}\log(1+\exp(\tilde{\omega}_{m}))$
for TTFs following logistics distributions, and
$\ell(\tilde{\sigma}\bm{y}-\bm{1}_{M}\tilde{\beta_{0}}-\bm{S}_{(4)}\tilde{\bm{\beta}_{1}})=-n\log\tilde{\sigma}-\Sigma_{m=1}^{M}\tilde{\omega}_{m}+\Sigma_{m=1}^{M}\exp(\tilde{\omega}_{m})$
for TTFs following SEV distributions.

The optimization algorithm to solve criterion (4) will be discussed later in
Section 3. The value of the weight $\alpha$ and the dimension of the low-
dimensional tensor subspace $\\{P_{1},P_{2},P_{3}\\}$ can be determined using
a model selection method such as cross-validation. Solving the optimization
criterion (4) using historical training data yields a set of basis matrices
$\hat{\bm{U}}_{1}\in\mathbb{R}^{P_{1}\times
I_{1}},\hat{\bm{U}}_{2}\in\mathbb{R}^{P_{2}\times
I_{2}},\hat{\bm{U}}_{3}\in\mathbb{R}^{P_{3}\times I_{3}}$, which contains
$P_{1}$ basis vectors of the $1$-mode linear space $\mathbb{R}^{I_{1}}$,
$P_{2}$ basis vectors of the $2$-mode linear space $\mathbb{R}^{I_{2}}$, and
$P_{3}$ basis vectors of the $3$-mode linear space $\mathbb{R}^{I_{3}}$,
respectively. The three linear subspaces form the low-dimensional tensor
subspace
$\mathbb{R}^{P_{1}}\otimes\mathbb{R}^{P_{2}}\otimes\mathbb{R}^{P_{3}}$
detected by the proposed supervised dimension reduction method.

###  2.3 Prognostic Model Construction and Real-Time TTF Prediction

In this subsection, we discuss how to build a prognostic model based on the
supervised dimension reduction method proposed in Section 2.2, and how to
predict the TTF distribution of an asset operating in the field using its
real-time degradation image data.

Similar to Section 2.2, we denote the training data set as
$\\{\mathcal{X}_{m}\in\mathbb{R}^{I_{1}\times I_{2}\times
D_{m}},y_{m}\\}_{m=1}^{M}$, where $M$ is the number of failed assets in the
training data set. Notice that $D_{m}$ might not be the same as
$D_{m^{\prime}}$ for two assets $m$ and $m^{\prime}$,
$m=1,\ldots,M,m^{\prime}=1,\ldots,M,m\neq m^{\prime}$. This is because
different asset’s failure times (i.e., TTFs) are different, and usually no
image data can be collected beyond an asset’s failure time since the asset is
stopped for maintenance or replace once it is failed. In addition to the
training data, we denote the degradation image stream of a test asset by time
$t$ as $\mathcal{X}_{t}\in\mathbb{R}^{I_{1}\times I_{2}\times I_{t}}$. The
objectives of this subsection include 1) constructing a prognostic model and
estimating its parameters using $\\{\mathcal{X}_{m},y_{m}\\}_{m=1}^{M}$ in the
training data set, and 2) using the estimated prognostic model to predict the
TTF (denoted as $\hat{y}_{t}$) of the test asset based on its degradation
image stream $\mathcal{X}_{t}$.

We first use the proposed supervised dimension reduction method to extract
low-dimensional features of both the training and test assets. Specifically,
as discussed in Section 2.2, we first construct a 4D tensor
$\mathcal{X}\in\mathbb{R}^{I_{1}\times I_{2}\times I_{3}\times M}$ using the
degradation image streams of the training assets, where
$I_{3}=\max(\\{D_{m}\\}_{m=1}^{M})$. Note that $\mathcal{X}$ is an incomplete
tensor no matter the image streams from the training assets are complete or
incomplete. This is because the TTFs of training assets are different, and
thus not all the training assets have $I_{3}$ images. To detect the low-
dimensional tensor subspace in which the high-dimensional degradation images
are embedded, we solve optimization criterion (4) by using training data
$\\{\mathcal{X},\bm{y}\\}$, where $\bm{y}=(y_{1},\ldots,y_{M})^{\top}$. This
yields basis matrices $\hat{\bm{U}}_{1}\in\mathbb{R}^{P_{1}\times
I_{1}},\hat{\bm{U}}_{2}\in\mathbb{R}^{P_{2}\times
I_{2}},\hat{\bm{U}}_{3}\in\mathbb{R}^{P_{3}\times I_{3}}$, which form the low-
dimensional tensor subspace
$\mathbb{R}^{P_{1}}\otimes\mathbb{R}^{P_{2}}\otimes\mathbb{R}^{P_{3}}$. To
extract the low-dimensional features of the training and test assets, we
expand the image streams in the low-dimensional tensor subspace
$\mathbb{R}^{P_{1}}\otimes\mathbb{R}^{P_{2}}\otimes\mathbb{R}^{P_{3}}$ using
the basis matrices $\hat{\bm{U}}_{1},\hat{\bm{U}}_{2},\hat{\bm{U}}_{3}$. This
is achieved by solving the following optimization criteria:

$\hat{\mathcal{S}}_{m}=\arg\min_{\bm{S}_{m}}||\mathcal{P}_{\Omega}(\mathcal{X}_{m}-\mathcal{S}_{m}\times_{1}\bm{U}_{1}^{\top}\times_{2}\bm{U}_{2}^{\top}\times_{3}\bm{U}_{3}^{\top})||_{F}^{2}.$
(5)
$\hat{\mathcal{S}}_{t}=\arg\min_{\bm{S}_{t}}||\mathcal{P}_{\Omega}(\mathcal{X}_{t}-\mathcal{S}_{t}\times_{1}\bm{U}_{1}^{\top}\times_{2}\bm{U}_{2}^{\top}\times_{3}\bm{U}_{3}^{\top})||_{F}^{2},$
(6)

where $\\{\hat{\mathcal{S}}_{m}\\}_{m=1}^{M}$ are the low-dimensional feature
tensors of the $M$ assets in the training data set, and
$\hat{\mathcal{S}}_{t}$ is the low-dimensional feature tensor of the test
asset.

Next, we construct a prognostic model using the low-dimensional feature
tensors of the $M$ assets in the training data set (i.e.,
$\\{\hat{\mathcal{S}}_{m}\\}_{m=1}^{M}$) along with their TTFs
$\\{y_{m}\\}_{m=1}^{M}$. Specifically, we build the following LLS regression
model:

$y_{m}=\gamma_{0}+vec(\hat{\mathcal{S}}_{m})^{\top}\bm{\gamma}_{1}+\sigma\epsilon_{m},$
(7)

where $vec(\hat{\mathcal{S}}_{m})$ is the vectorization of
$\hat{\mathcal{S}}_{m}$. $\gamma_{0}\in\mathbb{R}$ and
$\bm{\gamma}_{1}\in\mathbb{R}^{(P_{1}\times P_{2}\times P_{3})\times 1}$ are
the regression coefficients, $\sigma$ is the scale parameter, and
$\epsilon_{m}$ is the random noise term with a standard location-scale
probability density function $f(\epsilon)$. For example,
$f(\epsilon)=1/\sqrt{2\pi}\exp(-\epsilon^{2}/2)$ for a normal distribution and
$f(\epsilon)=\exp(\epsilon-\exp(\epsilon))$ for an SEV distribution. The
parameters in criterion (7) can be estimated by solving the following
optimization problem:

$\min_{\bm{y},\gamma_{0},\bm{\gamma}_{1},\sigma}\
\ell(\frac{\bm{y}-\bm{1}_{M}\gamma_{0}-\hat{\bm{S}}_{(4)}\bm{\gamma}_{1}}{\sigma}),$
(8)

where $\ell(\cdot)$ is the negative log-likelihood function of a location-
scale distribution,$\bm{y}=(y_{1},y_{2},\cdot\cdot\cdot,y_{m})^{\top}$ and
$\hat{\bm{S}}_{(4)}=(vec(\hat{\mathcal{S}}_{1})^{\top},vec(\hat{\mathcal{S}}_{2})^{\top},\cdot\cdot\cdot,vec(\hat{\mathcal{S}}_{M})^{\top})^{\top}$,
and $\ell(\cdot)$ is the negative log-likelihood function. We conduct the
following reparameterization to transform the optimization to be a convex one:
$\tilde{\sigma}=1/\sigma,\tilde{\gamma}_{0}=\gamma_{0}/\sigma,\tilde{\bm{\gamma}}_{1}=\bm{\gamma}_{1}/\sigma$:

$\\{\hat{\tilde{\gamma_{0}}},\hat{\tilde{\bm{\gamma}}}_{1},\hat{\tilde{\sigma}}\\}=\arg\min_{\tilde{\gamma}_{0},\tilde{\bm{\gamma}}_{1},\tilde{\sigma}}\
\ell(\tilde{\sigma}\bm{y}-\bm{1}_{M}\tilde{\gamma}_{0}-\hat{\bm{S}}_{(4)}\tilde{\bm{\gamma}}_{1}).\
$ (9)

Solving (9) provides the estimated parameters
$\\{\hat{\tilde{\gamma_{0}}},\hat{\tilde{\bm{\gamma}}}_{1},\hat{\tilde{\sigma}}\\}$,
which can be transformed back to the estimation of the parameters in the LLS
regression model:
$\hat{\gamma}_{0}=\hat{\tilde{\gamma}}_{0}/\hat{\tilde{\sigma}}$,
$\hat{\bm{\gamma}}_{1}=\hat{\tilde{\bm{\gamma}}}_{1}/\hat{\tilde{\sigma}}$ and
$\hat{\sigma}=1/\hat{\tilde{\sigma}}$. As a result, the fitted LLS regression
model is $\hat{y}_{m}\sim
LLS(\hat{\gamma}_{0}+vec(\hat{\mathcal{S}}_{m})^{\top}\hat{\bm{\gamma}}_{1},\hat{\sigma})$,
where
$\hat{\gamma}_{0}+vec(\hat{\mathcal{S}}_{m})^{\top}\hat{\bm{\gamma}}_{1}$ and
$\hat{\sigma}$ are respectively the estimated location and scale parameters.

Finally, we fed the extracted low-dimensional feature tensor of the test asset
into the estimated LLS regression model to predict the asset’s TTF
distribution: $\hat{y}_{t}\sim
LLS(\hat{\gamma}_{0}+vec(\hat{\mathcal{S}}_{t})^{\top}\hat{\bm{\gamma}}_{1},\hat{\sigma})$.

##  3 The Optimization Algorithm for General LLS Distributions

In this section, we develop a Block Updating Algorithm to solve the supervised
tensor dimension reduction method proposed in Section 2.2. The algorithm
splits the unknown parameters in criterion (4) into several blocks, and it
cyclically optimizes one block parameters in criterion (4) while keeping other
blocks fixed until convergence. The sub optimization problem for each block is
convex, so the convergence of the block updating algorithm is guaranteed. In
addition, we will show that if TTFs follow normal (or log-normal)
distributions, the sub optimization problem for each block has a closed-form
solution, no matter the degradation image streams are complete or incomplete.

###  3.1 The Block Updating Algorithm for General LLS Distributions

The Block Updating Algorithm first splits the unknown parameters in criterion
(4) into 5 blocks, i.e., $\bm{U}_{1},\bm{U}_{2},\bm{U}_{3},\mathcal{S}$, and
$\\{\tilde{\beta_{0}},\tilde{\bm{\beta}_{1}}\\}$. It then cyclically optimizes
one block of parameters each time while keeping other blocks fixed.

Specifically, at the $k$th iteration, $\bm{U}_{1}$ is updated by solving the
following optimization problem while keeping other blocks (i.e.,
$\bm{U}_{2}^{k-1},\bm{U}_{3}^{k-1},\tilde{\beta}_{0}^{k-1},\tilde{\bm{\beta}}_{1}^{k-1},\mathcal{S}^{k-1}$)
fixed:

$\begin{split}\bm{U}_{1}^{k}=\arg\min_{\bm{U}_{1}}~{}~{}\alpha||\mathcal{P}_{\Omega}(\mathcal{X}-\mathcal{S}^{k-1}\times_{1}\bm{U}_{1}^{\top}\times_{2}{\bm{U}_{2}^{k-1}}^{\top}&\times_{3}{\bm{U}_{3}^{k-1}}^{\top})||_{F}^{2}\\\
&+(1-\alpha)\ell(\tilde{\beta}_{0}^{k-1},\tilde{\bm{\beta}}_{1}^{k-1},\bm{S}_{(4)}^{k-1})\end{split}$
(10)

Similarly, the remaining blocks are updated as follows:

$\begin{split}\bm{U}_{2}^{k}=\arg\min_{\bm{U}_{2}}~{}~{}\alpha||\mathcal{P}_{\Omega}(\mathcal{X}-{\mathcal{S}}^{k-1}\times_{1}{\bm{U}_{1}^{k}}^{\top}\times_{2}{\bm{U}_{2}^{\top}}&\times_{3}{\bm{U}_{3}^{k-1}}^{\top})||_{F}^{2}\\\
&+(1-\alpha)\ell(\tilde{\beta}_{0}^{k-1},\tilde{\bm{\beta}}_{1}^{k-1},\bm{S}_{(4)}^{k-1})\end{split}$
(11)
$\begin{split}\bm{U}_{3}^{k}=\arg\min_{\bm{U}_{3}}~{}~{}\alpha||\mathcal{P}_{\Omega}(\mathcal{X}-{\mathcal{S}}^{k-1}\times_{1}{\bm{U}_{1}^{k}}^{\top}\times_{2}{\bm{U}_{2}^{k}}^{\top}&\times_{3}{\bm{U}_{3}^{\top}})||_{F}^{2}\\\
&+(1-\alpha)\ell(\tilde{\beta}_{0}^{k-1},\tilde{\bm{\beta}}_{1}^{k-1},\bm{S}_{(4)}^{k-1})\end{split}$
(12)
$\begin{split}\\{\tilde{\beta}_{0}^{k},\tilde{\bm{\beta}}_{1}^{k}\\}=\arg\min_{\tilde{\beta}_{0},\tilde{\bm{\beta}}_{1}}~{}~{}\alpha||\mathcal{P}_{\Omega}(\mathcal{X}-{\mathcal{S}}^{k-1}\times_{1}{\bm{U}_{1}^{k}}^{\top}\times_{2}{\bm{U}_{2}^{k}}^{\top}&\times_{3}{\bm{U}_{3}^{k}}^{\top})||_{F}^{2}\\\
&+(1-\alpha)\ell(\tilde{\beta}_{0},\tilde{\bm{\beta}}_{1},\bm{S}_{(4)}^{k-1})\end{split}$
(13)
$\begin{split}{\mathcal{S}}^{k}=\arg\min_{\mathcal{S}}~{}~{}\alpha||\mathcal{P}_{\Omega}(\mathcal{X}-\mathcal{S}\times_{1}{\bm{U}_{1}^{k}}^{\top}\times_{2}{\bm{U}_{2}^{k}}^{\top}&\times_{3}{\bm{U}_{3}^{k}}^{\top})||_{F}^{2}\\\
&+(1-\alpha)\ell(\tilde{\beta}_{0}^{k},\tilde{\bm{\beta}}_{1}^{k},\bm{S}_{(4)})\end{split}$
(14)

We summarize the Block Updating Algorithm in Algorithm 1 below. The
convergence criterion can be set as
$\Psi(\bm{U}_{1}^{k},\bm{U}_{2}^{k},\bm{U}_{3}^{k},\tilde{\beta}_{0}^{k},\tilde{\bm{\beta}}_{1}^{k},\mathcal{S}^{k})-\Psi(\bm{U}_{1}^{k+1},\bm{U}_{2}^{k+1},\bm{U}_{3}^{k+1},\tilde{\beta}_{0}^{k+1},\tilde{\bm{\beta}}_{1}^{k+1},\mathcal{S}^{k+1})<\epsilon$,
where $\Psi$ is the value of the objective function in criterion (4), and
$\epsilon$ is a small number. It is easy to show that all the sub problems
(10), (11), (12), (13), and (14) are convex. Thus, the Block Updating
Algorithm converges to a stationary point of criterion (4).

1

2Input: Tensor $\mathcal{X}$ constructed from the (incomplete) degradation
image streams of $M$ assets and the TTF vector $\bm{y}$; the dimension of the
low-dimensional tensor subspace $\\{P_{1},P_{2},P_{3}\\}$

3Initialization: Initialize
$(\bm{U}_{1}^{0},\bm{U}_{2}^{0},\bm{U}_{3}^{0},\tilde{\beta}_{0}^{0},\tilde{\bm{\beta}}_{1}^{0},\mathcal{S}^{0})$
randomly or heuristically

4While convergence criterion not met do

5$\bm{U}_{1}^{k}\longleftarrow$ (10)

6$\bm{U}_{2}^{k}\longleftarrow$ (11)

7$\bm{U}_{3}^{k}\longleftarrow$ (12)

8$(\tilde{\beta}_{0}^{k},\tilde{\bm{\beta}}_{1}^{k})\longleftarrow$ (13)

9$\mathcal{S}^{k}\longleftarrow$ (14)

10$k=k+1$

11End While

12Output: Basis matrices of the low-dimensional tensor subspace
$\\{\bm{U}_{1}^{k},\bm{U}_{2}^{k},\bm{U}_{3}^{k}\\}$

Algorithm 1 Block Updating Algorithm for solving criterion (4)

The initialization of Algorithm 1 can be conducted randomly or heuristically.
In this article, we propose a heuristic initialization method. Specifically,
MPCA (Lu et al., 2008) is applied to tensor $\mathcal{X}$, which provides
$\\{\bm{U}_{1}^{0},\bm{U}_{2}^{0},\bm{U}_{3}^{0}\\}$. Next, we compute
${\mathcal{S}}^{0}$ by solving
${\mathcal{S}}^{0}=\arg\min_{\mathcal{S}}~{}~{}||\mathcal{P}_{\Omega}(\mathcal{X}-\mathcal{S}\times_{1}{\bm{U}_{1}^{0}}^{\top}\times_{2}{\bm{U}_{2}^{0}}^{\top}\times_{3}{\bm{U}_{3}^{0}}^{\top})||_{F}^{2}$.
Finally, $\tilde{\beta}_{0}^{0},\tilde{\bm{\beta}}_{1}^{0}$ are computed by
solving
$\min_{\tilde{\beta}_{0}^{0},\tilde{\bm{\beta}}_{1}^{0},\tilde{\sigma}}\
\ell(\tilde{\sigma}\bm{y}-\bm{1}_{M}\tilde{\beta}_{0}^{0}-\bm{S}_{(4)}^{0}\tilde{\bm{\beta}}_{1}^{0})$,
where $\bm{S}_{(4)}^{0}$ is the mode-4 matricization of $\mathcal{S}^{0}$.
Algorithm 1 requires the dimension of the tensor subspace
$\\{P_{1},P_{2},P_{3}\\}$ to be known in advance. This can be seen as a model
selection problem such that the values of $\\{P_{1},P_{2},P_{3}\\}$ can be
determined using cross-validation. To be specific, we may try a certain number
of candidate values for $\\{P_{1},P_{2},P_{3}\\}$ to run Algorithm 1 to
extract low-dimensional features, which are then used to build the prognostic
model discussed in Section 2.3 for TTF prediction. The values which achieve
the smallest prediction error will be chosen. It is known that there usually
exists high spatio-temporal correlations among the degradation image streams
(Fang et al., 2019), so the dimension of the tensor subspace is usually low,
which helps reduce the computation intensity of model selection.

##  4 Closed-Form Solutions for Normal (or Lognormal) Distributions

In this section, we discuss the closed-form solutions for sub problems (10),
(11), (12), (13), and (14) for applications whose assets’ TTFs follow normal
(or lognormal) distributions. Specifically, if the TTFs follow normal (or
lognormal) distributions, we may replace the log-likelihood function-based
loss function in criterion (4) with the mean squared errors. As a result,
criteria (10), (11), (12),(13), and (14) can be re-expressed as follows:

$\begin{split}\bm{U}_{1}^{k}=\arg\min_{\bm{U}_{1}}~{}~{}\alpha||\mathcal{P}_{\Omega}(\mathcal{X}-\mathcal{S}^{k-1}\times_{1}\bm{U}_{1}^{\top}&\times_{2}{\bm{U}_{2}^{k-1}}^{\top}\times_{3}{\bm{U}_{3}^{k-1}}^{\top})||_{F}^{2}\\\
&+(1-\alpha)\|\bm{y}-\bm{1}_{M}\tilde{\beta}_{0}^{k-1}-\bm{S}_{(4)}^{k-1}\tilde{\bm{\beta}}_{1}^{k-1}\|_{2}^{2}\end{split}$
(15)
$\begin{split}\bm{U}_{2}^{k}=\arg\min_{\bm{U}_{2}}~{}~{}\alpha||\mathcal{P}_{\Omega}(\mathcal{X}-{\mathcal{S}}^{k-1}\times_{1}{\bm{U}_{1}^{k}}^{\top}&\times_{2}{\bm{U}_{2}^{\top}}\times_{3}{\bm{U}_{3}^{k-1}}^{\top})||_{F}^{2}\\\
&+(1-\alpha)\|\bm{y}-\bm{1}_{M}\tilde{\beta}_{0}^{k-1}-\bm{S}_{(4)}^{k-1}\tilde{\bm{\beta}}_{1}^{k-1}\|_{2}^{2}\end{split}$
(16)
$\begin{split}\bm{U}_{3}^{k}=\arg\min_{\bm{U}_{3}}~{}~{}\alpha||\mathcal{P}_{\Omega}(\mathcal{X}-{\mathcal{S}}^{k-1}\times_{1}{\bm{U}_{1}^{k}}^{\top}&\times_{2}{\bm{U}_{2}^{k}}^{\top}\times_{3}{\bm{U}_{3}^{\top}})||_{F}^{2}\\\
&+(1-\alpha)\|\bm{y}-\bm{1}_{M}\tilde{\beta}_{0}^{k-1}-\bm{S}_{(4)}^{k-1}\tilde{\bm{\beta}}_{1}^{k-1}\|_{2}^{2}\end{split}$
(17)
$\begin{split}\\{\tilde{\beta}_{0}^{k},\tilde{\bm{\beta}}_{1}^{k}\\}=\arg\min_{\tilde{\beta}_{0},\tilde{\bm{\beta}}_{1}}~{}~{}\alpha||\mathcal{P}_{\Omega}(\mathcal{X}-{\mathcal{S}}^{k-1}\times_{1}{\bm{U}_{1}^{k}}^{\top}&\times_{2}{\bm{U}_{2}^{k}}^{\top}\times_{3}{\bm{U}_{3}^{k}}^{\top})||_{F}^{2}\\\
&+(1-\alpha)\|\bm{y}-\bm{1}_{M}\tilde{\beta}_{0}-\bm{S}_{(4)}^{k-1}\tilde{\bm{\beta}}_{1}\|_{2}^{2}\end{split}$
(18)
$\begin{split}{\mathcal{S}}^{k}=\arg\min_{\mathcal{S}}~{}~{}\alpha||\mathcal{P}_{\Omega}(\mathcal{X}-\mathcal{S}\times_{1}{\bm{U}_{1}^{k}}^{\top}&\times_{2}{\bm{U}_{2}^{k}}^{\top}\times_{3}{\bm{U}_{3}^{k}}^{\top})||_{F}^{2}\\\
&+(1-\alpha)\|\bm{y}-\bm{1}_{M}\tilde{\beta}_{0}^{k}-\bm{S}_{(4)}\tilde{\bm{\beta}}_{1}^{k}\|_{2}^{2}\end{split}$
(19)

We will discuss the analytic solutions to problems (15), (16), (17),(18), and
(19) when degradation image streams are complete and incomplete in Sections
4.1 and 4.2, respectively. For simplicity, we will remove the superscripts $k$
and $k-1$.

###  4.1 Analytical Solutions for Complete Data

When degradation image streams are complete (i.e., the 4D image tensor
$\mathcal{X}$ in criterion (4) has no missing entries), we have the following
Propositions 1-4, which provide the analytical solutions to problems (15),
(16), (17), and (19), respectively. The solution to problem (18) is the
ordinary least squares, which does not need to be discussed.

###### Proposition 1

If the 4D tensor $\mathcal{X}$ has no missing values, the optimization problem
(15) has the following analytical solution

$\bm{U}_{1}=(\bm{X}_{(1)}\cdot\bm{S}_{{U}_{1}(1)}^{\top}\cdot(\bm{S}_{{U}_{1}(1)}\cdot\bm{S}_{{U}_{1}(1)}^{\top})^{-1})^{\top}$

where $\bm{X}_{(1)}$ is the mode-1 matricization of $\mathcal{X}$,
$\mathcal{S}_{{U}_{1}}=\mathcal{S}\times_{2}{\bm{U}}_{2}^{\top}\times_{3}{\bm{U}}_{3}^{\top}$,
$\bm{S}_{{U}_{1}(1)}$ is the mode-1 matricization of $\mathcal{S}_{{U}_{1}}$,
and the operator “$\cdot$” represents multiplication.

###### Proposition 2

If the 4D tensor $\mathcal{X}$ has no missing values, the optimization problem
(16) has the following analytical solution

$\bm{U}_{2}=(\bm{X}_{(2)}\cdot\bm{S}_{{U}_{2}(2)}^{\top}\cdot(\bm{S}_{{U}_{2}(2)}\cdot\bm{S}_{{U}_{2}(2)}^{\top})^{-1})^{\top}$

where $\bm{X}_{(2)}$ is the mode-2 matricization of $\mathcal{X}$,
$\mathcal{S}_{{U}_{2}}=\mathcal{S}\times_{1}{\bm{U}}_{1}^{\top}\times_{3}{\bm{U}}_{3}^{\top}$,
and $\bm{S}_{{U}_{2}(2)}$ is the mode-2 matricization of
$\mathcal{S}_{{U}_{2}}$.

###### Proposition 3

If the 4D tensor $\mathcal{X}$ has no missing values, the optimization problem
(17) has the following analytical solution

$\bm{U}_{3}=(\bm{X}_{(3)}\cdot\bm{S}_{{U}_{3}(3)}^{\top}\cdot(\bm{S}_{{U}_{3}(3)}\cdot\bm{S}_{{U}_{3}(3)}^{\top})^{-1})^{\top}$

where $\bm{X}_{(3)}$ is the mode-3 matricization of $\mathcal{X}$,
$\mathcal{S}_{{U}_{3}}=\mathcal{S}\times_{1}{\bm{U}}_{1}^{\top}\times_{2}{\bm{U}}_{2}^{\top}$,
and $\bm{S}_{{U}_{3}(3)}$ is the mode-3 matricization of
$\mathcal{S}_{{U}_{3}}$.

###### Proposition 4

If the 4D tensor $\mathcal{X}$ has no missing values, the optimization problem
(19) has the following analytical solution

$\begin{split}\bm{S}_{(4)}&=\big{[}\alpha\cdot\bm{X}_{(4)}\cdot(\bm{U}_{3}\otimes\bm{U}_{2}\otimes\bm{U}_{1})^{\top}+(1-\alpha)\cdot(\bm{y}-\bm{1}_{M}\cdot\tilde{\beta}_{0})\cdot\tilde{\bm{\beta}}_{1}^{\top}\big{]}\cdot\\\
&\big{[}\alpha\cdot(\bm{U}_{3}\otimes\bm{U}_{2}\otimes\bm{U}_{1})\cdot(\bm{U}_{3}\otimes\bm{U}_{2}\otimes\bm{U}_{1})^{\top}+(1-\alpha)\cdot\tilde{\bm{\beta}}_{1}\cdot\tilde{\bm{\beta}}_{1}^{\top}\big{]}^{-1}\end{split}$

where $\bm{X}_{(4)}$ is the mode-4 matricization of $\mathcal{X}$,
$\bm{S}_{(4)}$ is the mode-4 matricization of $\mathcal{S}$.

The proof of Propositions 1-4 can be found in the Appendix.

###  4.2 Analytical Solutions for Incomplete Data

In this subsection, we discuss the closed-form solutions for optimization
problems (15), (16), (17), (18), and (19) when the degradation tensor
$\mathcal{X}$ in criterion (4) is incomplete. We consider two data missing
patterns: entry-wise missing and image-wise missing.

Entry-wise missing means that any entry of
$\mathcal{X}\in\mathbb{R}^{I_{1}\times I_{2}\times I_{3}\times M}$ can be
missing. Thus, the indices of missing entries are denoted as a subset
$\Omega\subseteq\\{(i_{1},i_{2},i_{3},m),1\leq i_{1}\leq I_{1},1\leq i_{2}\leq
I_{2},1\leq i_{3}\leq I_{3},1\leq m\leq M\\}$. For image-based applications,
the missing elements in tensor $\mathcal{X}$ are images but not entries, as
illustrated in Figure 2(a)(b), which yields image-wise missing. Thus, the
subset containing the indices of missing entries can be expressed as
$\Omega\subseteq\\{(:,:,i_{3},m),1\leq i_{3}\leq I_{3},1\leq m\leq M\\}$,
where “$:$” denotes all the indices in a specific dimension.

Since the solution to problem (18) is the ordinary least squares, which does
not need to be discussed, we will only discuss the solutions to problems (15),
(16), (17), and (19).

####  4.2.1 Entry-wise Missing

When the tensor $\mathcal{X}$ has an entry-wise missing structure, there is no
direct closed-form solutions for optimization criteria (15), (16), (17), and
(19). However, we may decompose them into multiple sub optimization problems,
each of which yields an analytical solution. To be specific, we first give the
following lemma.

###### Lemma 1

Let $\bm{A}\in\mathbb{R}^{M\times N}$, $\bm{C}\in\mathbb{R}^{P\times N}$, and
$\bm{B}\in\mathbb{R}^{M\times P}$, the solution to criterion
$\arg\min_{B}\|\bm{A}-\bm{B}\bm{C}\|_{F}^{2}$ can be found by solving each row
of $\bm{B}$ separately–that is–solving $\\{\bm{b}_{m}\\}_{m=1}^{M}$ as
follows:

$\arg\min_{b_{m}}\|\bm{a}_{m}-\bm{b}_{m}\bm{C}\|_{F}^{2},\quad m=1,\ldots,M$

where $\bm{a}_{m}\in\mathbb{R}^{1\times N}$ is the $m$th row of $\bm{A}$,
$\bm{b}_{m}\in\mathbb{R}^{1\times P}$ is the $m$th row of $\bm{B}$.

The proof of Lemma 1 can be found in the Appendix. Lemma 1 enables us to solve
each column of matrices $\bm{U}_{1}$, $\bm{U}_{2}$, and $\bm{U}_{3}$ as well
as each row of matrix $\bm{S}_{(4)}$ (the mode-4 matricization of tensor
$\mathcal{S}$) separately.

Denote the $i_{1}$th column of matrix $\bm{U}_{1}\in\mathbb{R}^{P_{1}\times
I_{1}}$ as $\bm{u}_{1}^{i_{1}}\in\mathbb{R}^{P_{1}\times
1},i_{1}=1,\ldots,I_{1}$, we replace optimization problem (15) with $I_{1}$
sub problems by separately optimizing
$\bm{u}_{1}^{1},\bm{u}_{1}^{2},\ldots,\bm{u}_{1}^{I_{1}}$. Proposition 5
suggests that there is an analytical solution when optimizing
$\bm{u}_{1}^{i_{1}}$.

###### Proposition 5

When optimizing the $i_{1}$th column of $\bm{U}_{1}$ in Problem (15), we have
the following analytical solution

$\bm{u}_{1}^{i_{1}}=(\bm{x}_{(1)}^{i_{1},{\pi_{i_{1}}}}\cdot{\bm{S}_{{U}_{1}(1)}^{\pi_{i_{1}}}}^{\top}\cdot(\bm{S}_{{U}_{1}(1)}^{\pi_{i_{1}}}\cdot{\bm{S}_{{U}_{1}(1)}^{\pi_{i_{1}}}}^{\top})^{-1})^{\top}$

where $\bm{x}_{(1)}^{i_{1}}$ denotes the $i_{1}$th row of $\bm{X}_{(1)}$,
$\pi_{i_{1}}$ is a set consisting of the indices of available elements of
$\bm{x}_{(1)}^{i_{1}}$, $\bm{x}_{(1)}^{i_{1},\pi_{i_{1}}}$ denotes a vector
consisting of the available elements in the $i_{1}$th row of $\bm{X}_{(1)}$,
$\mathcal{S}_{{U}_{1}}=\mathcal{S}\times_{2}\bm{U}_{2}^{\top}\times_{3}\bm{U}_{3}^{\top}$,
$\bm{S}_{{U}_{1}(1)}$ is the mode-1 matricization of $\mathcal{S}_{{U}_{1}}$,
and $\bm{S}_{{U}_{1}(1)}^{\pi_{i_{1}}}$ denotes a matrix comprises the
$\pi_{i_{1}}$ columns of $\bm{S}_{{U}_{1}(1)}$.

Denote the $i_{2}$th column of matrix $\bm{U}_{2}\in\mathbb{R}^{P_{2}\times
I_{2}}$ as $\bm{u}_{2}^{i_{2}}\in\mathbb{R}^{P_{2}\times
1},i_{2}=1,\ldots,I_{2}$, we replace optimization problem (16) with $I_{2}$
sub problems by separately optimizing
$\bm{u}_{2}^{1},\bm{u}_{2}^{2},\ldots,\bm{u}_{2}^{I_{2}}$. Proposition 6
suggests that there is an analytical solution when optimizing
$\bm{u}_{2}^{i_{2}}$.

###### Proposition 6

When optimizing the $i_{2}$th column of $\bm{U}_{2}$ in Problem (16), we have
the following analytical solution

$\bm{u}_{2}^{i_{2}}=(\bm{x}_{(2)}^{i_{2},\pi_{i_{2}}}\cdot{\bm{S}_{{U}_{2}(2)}^{\pi_{i_{2}}}}^{\top}\cdot(\bm{S}_{{U}_{2}(2)}^{\pi_{i_{2}}}\cdot{\bm{S}_{{U}_{2}(2)}^{\pi_{i_{2}}}}^{\top})^{-1})^{\top}$

where $\bm{x}_{(2)}^{i_{2}}$ denotes the $i_{2}$th row of $\bm{X}_{(2)}$,
$\pi_{i_{2}}$ is a set consisting of the indices of available elements of
$\bm{x}_{(2)}^{i_{2}}$, $\bm{x}_{(2)}^{i_{2},\pi_{i_{2}}}$ is a vector
consisting of the available elements in the $i_{2}$th column of
$\bm{X}_{(2)}$,
$\mathcal{S}_{{U}_{2}}=\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{3}\bm{U}_{3}^{\top}$,
$\bm{S}_{{U}_{2}(2)}$ is the mode-2 matricization of $\mathcal{S}_{{U}_{2}}$,
and $\bm{S}_{{U}_{2}(2)}^{\pi_{i_{2}}}$ denotes a matrix comprises the
$\pi_{i_{2}}$ columns of $\bm{S}_{{U}_{2}(2)}$.

Denote the $i_{3}$th column of matrix $\bm{U}_{3}\in\mathbb{R}^{P_{3}\times
I_{3}}$ as $\bm{u}_{3}^{i_{3}}\in\mathbb{R}^{P_{3}\times
1},i_{3}=1,\ldots,I_{3}$, we replace optimization problem (17) with $I_{3}$
sub problems by separately optimizing
$\bm{u}_{3}^{1},\bm{u}_{3}^{2},\ldots,\bm{u}_{3}^{I_{3}}$. Proposition 7
suggests that there is an analytical solution when optimizing
$\bm{u}_{3}^{i_{3}}$.

###### Proposition 7

When optimizing the $i_{3}$th column of $\bm{U}_{3}$ in Problem (17), we have
the following analytical solution

$\bm{u}_{3}^{i_{3}}=(\bm{x}_{(3)}^{i_{3},\pi_{i_{3}}}\cdot{\bm{S}_{{U}_{3}(3)}^{\pi_{i_{3}}}}^{\top}\cdot(\bm{S}_{{U}_{3}(3)}^{\pi_{i_{3}}}\cdot{\bm{S}_{{U}_{3}(3)}^{\pi_{i_{3}}}}^{\top})^{-1})^{\top}$

where $\bm{x}_{(3)}^{i_{3}}$ denotes the $i_{3}$th row of $\bm{X}_{(3)}$,
$\pi_{i_{3}}$ is a set consisting of the indices of available elements of
$\bm{x}_{(3)}^{i_{3}}$, $\bm{x}_{(3)}^{i_{3},\pi_{i_{3}}}$ is a vector
consisting of the available elements in the $i_{3}$th row of $\bm{X}_{(3)}$,
$\mathcal{S}_{{U}_{3}}=\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{2}\bm{U}_{2}^{\top}$,
$\bm{S}_{{U}_{3}(3)}$ is the mode-3 matricization of $\mathcal{S}_{{U}_{3}}$,
and $\bm{S}_{{U}_{3}(3)}^{\pi_{i_{3}}}$ denotes a matrix comprises the
$\pi_{i_{3}}$ columns of $\bm{S}_{{U}_{3}(3)}$.

Similarly, we denote the $m$th row of matrix
$\bm{S}_{(4)}\in\mathbb{R}^{M\times(P_{1}\times P_{2}\times P_{3})}$ as
$\bm{s}_{(4)}^{m}\in\mathbb{R}^{1\times(P_{1}\times P_{2}\times P_{3})}$,
$m=1,\ldots,M$, and replace optimization problem (19) with $M$ sub
problems–that is–we separately optimize
$\bm{s}_{(4)}^{1},\bm{s}_{(4)}^{2},\ldots,\bm{s}_{(4)}^{M}$. Proposition 8
suggests that there is an analytical solution when optimizing
$\bm{s}_{(4)}^{m}$.

###### Proposition 8

When optimizing the $m$th row of matrix $\bm{S}_{(4)}$ in Problem (19), we
have the following analytical solution

$\begin{split}\bm{s}_{(4)}^{m}&=\big{[}\alpha\cdot\bm{x}_{(4)}^{m,\pi_{m}}\cdot(\bm{U}_{3}\otimes\bm{U}_{2}\otimes\bm{U}_{1})^{\pi_{m}\top}+(1-\alpha)\cdot(y_{m}-\tilde{\beta}_{0})\cdot\tilde{\bm{\beta}}_{1}^{\top}\big{]}\cdot\\\
&\big{[}\alpha\cdot(\bm{U}_{3}\otimes\bm{U}_{2}\otimes\bm{U}_{1})^{\pi_{m}}\cdot(\bm{U}_{3}\otimes\bm{U}_{2}\otimes\bm{U}_{1})^{\pi_{m}\top}+(1-\alpha)\cdot\tilde{\bm{\beta}}_{1}\cdot\tilde{\bm{\beta}}_{1}^{\top}\big{]}^{-1}\end{split}$

where $\bm{x}_{(4)}^{m}$ represents the $m$th row of $\bm{X}_{(4)}$, $\pi_{m}$
denotes the set consisting of the indices of available elements in
$\bm{x}_{(4)}^{m}$, $\bm{x}_{(4)}^{m,\pi_{m}}$ is a vector consisting of the
available elements in the $m$th row of $\bm{X}_{(4)}$, and
$(\bm{U}_{3}\otimes\bm{U}_{2}\otimes\bm{U}_{1})^{\pi_{m}}$ denotes a matrix
comprising the $\pi_{m}$ columns of the matrix
$\bm{U}_{3}\otimes\bm{U}_{2}\otimes\bm{U}_{1}$.

The proof of Propositions 5-8 can be found in the Appendix.

####  4.2.2 Image-wise Missing

When the tensor $\mathcal{X}$ has an image-wise missing structure,
Propositions 9 and 10 suggest that optimization criteria (15) and (16) have
analytical solutions.

###### Proposition 9

If the indices of tensor $\mathcal{X}$’s missing entries can be denoted as
$\Omega\subseteq\\{(:,:,i_{3},m)$, $1\leq i_{3}\leq I_{3},1\leq m\leq M\\}$,
where “$:$” denotes all the indices in a dimension, then $\mathcal{X}$’s
mode-1 matricization $\bm{X}_{(1)}$ has missing columns. Let $\pi$ be the set
consisting of the indices of available columns in $\bm{X}_{(1)}$, then the
optimization problem (15) has the following analytical solution

$\bm{U}_{1}=(\bm{X}_{(1)}^{\pi}\cdot{\bm{S}_{{U}_{1}(1)}^{\pi}}^{\top}\cdot(\bm{S}_{{U}_{1}(1)}^{\pi}\cdot{\bm{S}_{{U}_{1}(1)}^{\pi}}^{\top})^{-1})^{\top},$

where $\bm{X}_{(1)}^{\pi}$ is a matrix consisting of the $\pi$ columns of
$\bm{X}_{(1)}$,
$\mathcal{S}_{{U}_{1}}=\mathcal{S}\times_{2}\bm{U}_{2}^{\top}\times_{3}\bm{U}_{3}^{\top}$,
$\bm{S}_{{U}_{1}(1)}$ is the mode-1 matricization of $\mathcal{S}_{{U}_{1}}$,
and $\bm{S}_{{U}_{1}(1)}^{\pi}$ denotes a matrix constituting the $\pi$
columns of $\bm{S}_{{U}_{1}(1)}$.

###### Proposition 10

If the indices of tensor $\mathcal{X}$’s missing entries can be denoted as
$\Omega\subseteq\\{(:,:,i_{3},m)$, $1\leq i_{3}\leq I_{3},1\leq m\leq M\\}$,
where “$:$” denotes all the indices in a dimension, then $\mathcal{X}$’s
mode-2 matricization $\bm{X}_{(2)}$ has missing columns. Let $\pi$ be the set
consisting of the indices of available columns in $\bm{X}_{(2)}$, then the
optimization problem (16) has the following analytical solution

$\bm{U}_{2}=(\bm{X}_{(2)}^{\pi}\cdot{\bm{S}_{{U}_{2}(2)}^{\pi}}^{\top}\cdot(\bm{S}_{{U}_{2}(2)}^{\pi}\cdot{\bm{S}_{{U}_{2}(2)}^{\pi}}^{\top})^{-1})^{\top}$

where $\bm{X}_{(2)}^{\pi}$ is a matrix consisting of the $\pi$ columns of
$\bm{X}_{(2)}$, $\bm{S}_{{U}_{2}(2)}$ is the mode-2 matricization of
$\mathcal{S}_{{U}_{2}}$, and $\bm{S}_{{U}_{2}(2)}^{\pi}$ denotes a matrix
constituting the $\pi$ columns of $\bm{S}_{{U}_{2}(2)}$.

The proof of Propositions 9-10 can be found in the Appendix.

##  5 Numerical Studies

In this section, we validate the effectiveness of our proposed supervised
tensor dimension reduction-based prognostic model using simulated data.

###  5.1 Data Generation

We generate degradation image streams for $500$ assets. The image stream from
asset $m$, which is denoted by $\mathcal{X}_{m}(x,y,t),~{}m=1,2,\ldots,500$,
is generated from the following heat transfer equation:

$\frac{\partial\mathcal{X}_{m}(x,y,t)}{\partial
t}=\alpha_{m}\Big{(}\frac{\partial^{2}\mathcal{X}_{m}}{\partial
x^{2}}+\frac{\partial^{2}\mathcal{X}_{m}}{\partial y^{2}}\Big{)},$ (20)

where $(x,y),0\leq x,y\leq 0.2$, represents the location of each image pixel.
$\alpha_{m}$ is the thermal diffusivity coefficient, which is randomly
generated from a uniform distribution $\mathcal{U}(0.5\times 10^{-4},1\times
10^{-4})$. $t$ is the time index. The initial and boundary conditions are set
such that $\mathcal{X}|_{t=1}=0$ and
$\mathcal{X}_{m}|_{x=0}=\mathcal{X}_{m}|_{x=0.2}=\mathcal{X}_{m}|_{y=0}=\mathcal{X}_{m}|_{y=0.2}=30$.
At each time $t$, the image is recorded at locations
$x=\frac{j}{n+1},y=\frac{k}{n+1},j,k=1,\ldots,n$, resulting in an $n\times n$
matrix. Here, we set $n=21$ and $t=1,2,\ldots,150$, which yields 150 images of
size 21 $\times$ 21 for each asset. This implies that the degradation image
stream of each asset can be represented by a 21 $\times$ 21 $\times$ 150
tensor. In addition, an independent and identically distributed random noise
$\epsilon\sim N(0,0.1)$ is added to each pixel. Figure 2(a) demonstrates an
example of some images with and without noise from one of the assets simulated
in this study.

(a) Without noise (b) With noise Figure 2: Simulated degradation images based
on heat transfer process.

To determine the TTF of an asset, we first transform the asset’s $21\times
21\times 150$ tensor to a $1\times 150$ time series by taking the average
pixel intensity of each image. The time series signal indicates how the
average heat of the asset involves over time. Next, we let the TTF of the
asset be the time point where the amplitude of the time series signal crosses
a pre-defined soft failure threshold, which is set as $23$ in this study.
Since the images of different assets are generated with different thermal
diffusivity coefficients, the time points where their time series signals go
beyond the threshold may be different. Thus, the TTF of different asset may
also be different. To mimic reality, we truncate the image stream of each
asset by keeping only the images observed before its TTF. In other words, any
images observed after an asset’s TTF are removed from the image tensor of the
asset. Such a truncation is normal in reality since an asset usually gets
maintained or replaced once its degradation signal crosses the soft failure
threshold. Consequently, the third dimension of the tensor of different assets
might be different. In addition, to reduce the computation load, we keep one
of every 10 images in the truncated image stream of each asset.

###  5.2 The Benchmark and Performance Comparison

We randomly split the generated data into a training data set consisting of
$400$ assets and a test data set consisting of the remaining $100$ assets. To
test the robustness of the proposed method, we consider four levels of data
incompleteness: (1) $0\%$ missing, (2) $10\%$ missing, (3) $50\%$ missing, (4)
$90\%$ missing. For the first scenario, (1) $0\%$ missing, we use all the
generated data for model training and testing. Please notice that even though
all the available images are used, the image tensor $\mathcal{X}$ is still
incomplete due to failure time truncation–that is–different assets may have
different TTFs and thus different number of images (see the discussion in the
second paragraph of Section 2.3). For the remaining scenarios, we randomly
remove some images from each asset’s image stream. For example, when $10\%$
missing, we randomly remove $10\%$ of the images (rounding to the nearest
integer) from the image stream of each asset.

We compare the performance of our proposed method with an unsupervised tensor
dimension reduction-based benchmark. Considering image streams are incomplete,
the baseline model first applies a tensor completion method known as TMac
developed by Xu et al. (2013) to impute the missing values of the image
tensor. Next, an unsupervised tensor dimension reduction method, MPCA (Lu et
al., 2008), is employed to reduce the dimension of the imputed image tensor to
reduce dimension and extract low-dimension features, which are then used to
build an LLS-based prognostic model as we discussed in Section 2.3. MPCA is a
widely used dimension reduction method for tensor data. It projects a high-
dimensional tensor into a subspace but maximizes the total tensor scatter
which is assumed to measure the variations in the original tensor objects. Lu
et al. (2008) proposed a fraction-of-variation-explained (FVE) method to
determine the dimension of the low-dimension tensor subspace/features, which
represents the percentage of variation of the original high-dimensional tensor
is preserved by the low-dimensional tensor features. Since the optimal FVE
suggested by Lu et al. (2008) is $97\%$, we will first set FVE as $97\%$ in
this study, and the corresponding baseline model is designated as “MPCA
(97%)”. In addition to the FVE method, we also use cross-validation (CV) to
select an appropriate dimension for the tensor subspace. Specifically, we use
the training data to conduct a 10-fold CV for various combinations of
$(P_{1},P_{2},P_{3})$, where $P_{1}=1,\ldots,4$, $P_{2}=1,\ldots,4$, and
$P_{3}=1,\ldots,4$. The baseline model is referred to as “MPCA_CV”.

We also use 10-fold CV to determine the value of the weight parameter $\alpha$
and the appropriate dimension of the tensor subspace of our proposed method.
Same as “MPCA_CV”, we try various combinations of $(P_{1},P_{2},P_{3})$ by
letting $P_{1},P_{2},P_{3}$ vary from 1 to 4. We use the heuristic method
discussed in Section 3.1 to initialize the block updating algorithm. In this
study, we use lognormal regression to build the prognostic model. The proposed
method is denoted as “Proposed_CV”. The prediction errors of our proposed
method and two benchmarks are calculated by using the equation below and
reported in Figures 4, 4, 6, and 6.

$\text{Prediction Error}=\frac{|\text{Estimated TTF}-\text{True
TTF}|}{\text{True TTF}}.$ (21)

###  5.3 Results and analysis

Figure 4 reports the prediction errors of the two benchmarks and our proposed
method when data is complete, which means no image is removed on purpose.
Figure 4 shows the prediction errors when 10% elements in the 3rd mode (time)
of degradation image streams are missing, while Figures 6 and 6 demonstrate
the errors when 50% and 90% images are missing, respectively.

Figures 4, 4, 6, and 6 illustrate that our proposed method outperforms the
benchmarks under all data missing rates. For example, when the degradation
image signals are complete, the median absolute prediction errors (and the
Interquartile Ranges, i.e., IQRs) of the proposed method and the two
benchmarks are 0.003 (0.003), 0.027 (0.035), and 0.025 (0.033), respectively;
when 10% images are missing, the median absolute prediction errors (and IQRs)
of the three methods are respectively 0.019 (0.017), 0.058 (0.067), and 0.053
(0.063); when 50% of images are missing, they are 0.052 (0.084), 0.302
(0.405), and 0.104 (0.168). We believe this is because our proposed method
applies historical TTFs to supervise the low-dimensional tensor dimension
reduction, and thus the extracted features are more effective for failure time
prediction. Unlike our method, the two baseline models use MPCA, an
unsupervised tensor dimension reduction method, for feature extraction. Since
the extracted features are only determined by the image streams, and no TTF
gets involved, they are not as effective as the features extracted by our
proposed method, and thus their failure time prediction accuracy and precision
are compromised.

The figures 4, 4, 6, and 6 also suggest that the performances of all the three
models deteriorate, and the superiority of our proposed method over the two
benchmarks decreases, with the increase of data missing rate. For example,
when data is complete, the median absolute prediction errors (and IQRs) of
“Proposed_CV” and “MPCA_CV” are 0.003 (0.003) and 0.025 (0.033), respectively;
when the missing rate increases to 90%, they are respectively 0.13 (0.21) and
0.16 (0.19), which are almost comparable. This is reasonable since the
performance of all the models are compromised more when more data are missing.
In addition, no model will perform well if a high percentage (say more than
90%) of data are missing since it implies that very limited useful degradation
information is available for modeling.

Figure 3: Prediction errors when data is complete in Numerical Study.

Figure 4: Prediction errors when 10% data is missing in Numerical Study.

Figure 5: Prediction errors when 50% data is missing in Numerical Study.

Figure 6: Prediction errors when 90% data is missing in Numerical Study.

The Figures 4-6 also demonstrate that “MPCA_CV” always outperforms “MPCA
(97%)”, and the superiority of “MPCA_CV” is augmented with the increase of
data missing rate. For instance, when 10% images are missing, the median
absolute prediction errors (and IQR) of “MPCA (97%)” and “MPCA_CV” are 0.058
(0.067) and 0.053 (0.063), respectively; when the missing rate is 50%, they
are 0.302 (0.405) and 0.104 (0.168). One of the possible reasons is that “MPCA
(97%)” determines the dimension of the tensor subspace by setting the “FVE” as
97%, which usually results in relatively high-dimensional features, although
the dimension is smaller than that of the original image tensor. Relatively
high-dimensional features implies an insufficient dimension reduction. In
addition, it means the number of parameters in the subsequent LLS-based
prognostic model is relatively large, which poses estimation challenges given
that the number of samples (assets) for model training is limited.

##  6 Case Study

In this section, we use degradation image streams obtained from a rotating
machinery test bed to validate the effectiveness of our proposed method. The
test bed is designed to perform accelerated degradation tests on rolling
element thrust bearings. Specifically, bearings were run from brand new to
failure. An FLIR T300 infrared camera was used to monitor the degradation
process and collect degradation images over time. In the meanwhile, an
accelerometer was mounted on the test bed to monitor the vibration of the
bearing, and the failure time is defined as the time point where the amplitude
of defective vibration frequencies crosses a threshold based on ISO standards
for machine vibration. The data set consists of $284$ degradation image
streams and their corresponding TTFs, and each image has $40\times 20$ pixels.
As an illustration, a sequence of images obtained at different (ordered) time
periods of one of the bearing are shown in Figure 7. More details about the
experimental setup and the data set can be found in Gebraeel et al. (2009) and
Fang et al. (2019).

Figure 7: An illustration of one infrared degradation image stream.

We use 5-fold cross validation to evaluate the performance of our proposed
model and the two benchmarks discussed in Section 5. Similar to the simulation
study, we conduct 10-fold cross validation to determine the optimal weight
parameter in criterion (4) and the most appropriate dimension of the tensor
subspace. In addition, we also consider four levels of data incompleteness:
(1) $0\%$ missing (i.e., complete), (2) $10\%$ missing, (3) $50\%$ missing,
and (4) $90\%$ missing. Figure 9 illustrates the absolute prediction errors
when degradation image streams are complete. Figure 9 shows prediction errors
when 10% of the images of each bearing are missing. Figures 11 and 11
demonstrate the absolute prediction errors when the missing rate are 50% and
90%, respectively.

Figure 8: Prediction errors when data is complete in Case Study.

Figure 9: Prediction errors when 10% data is missing in Case Study.

Figure 10: Prediction errors when 50% data is missing in Case Study.

Figure 11: Prediction errors when 90% data is missing in Case Study.

Similar to the discovery in the numerical study in Section 5, Figures 9, 9,
11, and 11 indicate that our proposed method constantly works better than the
two benchmarks under all the 4 data missing rates. For example, the median
absolute prediction errors (and IQRs) of our proposed method and the two
benchmarks are 0.03 (0.04), 0.3 (0.24), and 0.1 (0.16), respectively, when the
degradation image streams are complete. When 50% of the images are missing,
the median absolute prediction errors (and IQR) of are respectively 0.09
(0.17), 0.62 (0.65), and 0.15 (0.19). We believe this is because our proposed
model is a supervised dimension reduction-based method, which uses TTF
information to supervise the defection of the low-dimensional tensor subspace,
while the benchmarks are unsupervised dimension reduction-based methods
without TTF information involved. Since our method considers TTF information
when detecting the tensor subspace, the extracted features are more effective
for failure time prediction.

Figure 9, 9, 11, 11 also show that the prediction errors of all the 3 methods
increase with the increase of data missing rates. For example, when the
missing rates are 0%, 10%, 50%, and 90%, the median absolute prediction errors
(and IQRs) of “MPCA (97%)” are 0.3 (0.24), 0.49 (0.42), 0.62 (0.65), and 0.83
(0.78), respectively, while they are respectively 0.09 (0.11), 0.1 (0.16),
0.15 (0.19), and 0.31 (0.22) for “MPCA_CV”, and 0.03 (0.04), 0.05 (0.08), 0.09
(0.17), and 0.29 (0.21) for our proposed method. This is reasonable since a
higher data missing rate means less useful degradation information and thus a
worse model performance. In addition, we observe that the superiority of our
proposed method over the two benchmarks decrease with the increase of data
missing rates. For example, the prediction accuracy of our method and “MPCA
(97%)” are comparable when the data missing rate is 90%. We again believe this
is because not much useful information is available when data is highly
incomplete and none of the two models perform well with such limited data.

We also observe that “MPCA_CV” always outperforms “MPCA (97%)”, and the
superiority of “MPCA_CV” is augmented with the increase of data missing rate.
For example, when 10% images are missing, the median absolute prediction
errors (and IQR) of “MPCA (97%)” and MPCA_CV are 0.49 (0.42) and 0.1 (0.16),
respectively; when the missing rate is 50%, the they are 0.62 (0.65) and 0.15
(0.19). We again believe this is because “MPCA (97%)” determines the dimension
of the tensor subspace by setting the “FVE” as 97%, which results in
relatively high-dimensional features due to the insufficient dimension
reduction. Also, it results in a parameter estimation challenge since the
number of parameters to be estimated in the prognostic model is relatively
large comparing to the limited number of historical samples for model
training. This suggests that cross validation is a better method to determine
the dimension of the tensor subspace, especially when we do not have enough
number of samples for model training.

##  7 Conclusions

This paper proposed a supervised tensor dimension reduction-based prognostic
model for applications with incomplete degradation imaging data. This is
achieved by first developing a new supervised tensor dimension reduction
method that reduces the dimension of incomplete high-dimensional degradation
image streams and provides low-dimensional tensor features, which are then
used to build a prognostic model based on (log)-location-scale regression.

The supervised tensor dimension reduction method uses historical TTFs to
supervise the detection of a low-dimensional tensor subspace to reduce the
dimension of incomplete high-dimensional image streams. Mathematically, it is
formulated as an optimization criterion that combines a feature extraction
term and a regression term. The feature extraction term focuses on identifying
a tensor space to extract low-dimensional tensor features from high-
dimensional image streams. The regression term regresses failure times against
the features extracted by the first term using LLS regression. By jointly
optimizing the two terms, it is expected to detect an appropriate tensor
subspace such that the extracted features are effective for TTF prediction. To
estimate the parameters of the supervised dimension reduction method, we
developed a Block Updating Algorithm for applications where TTFs follow
distributions in the (log)-location-scale family. The algorithm works by
splitting the parameters into several blocks and cyclically optimizing one
block of parameters while keeping other blocks fixed until convergence. In
addition, we showed that if TTFs follow normal or lognormal distributions,
there is a closed form solution when optimizing each block of the parameters,
no matter the imaging data is complete or incomplete.

Simulated data as well as a data set from rotating machinery were used to
validate the effectiveness of our proposed method. The results showed that our
proposed prognostic method consistently outperforms the unsupervised tensor
reduction-based benchmarks under various data missing rates. This validated
the benefits and importance of using failure time information to supervise the
dimension reduction of high-dimensional degradation image streams when
building prognostic models.

Appendix

##  1 Proof of Proposition 1

The original optimization problem is

$\arg\min_{\bm{U}_{1}}~{}~{}\alpha||\mathcal{P}_{\Omega}(\mathcal{X}-\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{2}\bm{U}_{2}^{\top}\times_{3}\bm{U}_{3}^{\top})||_{F}^{2}+(1-\alpha)||\bm{y}-\bm{1}_{M}\cdot\tilde{\beta}_{0}-\bm{S}_{(4)}\cdot\tilde{\bm{\beta}}_{1}||_{F}^{2},$

which is equivalent to the following problem when data is complete:

$\arg\min_{\bm{U}_{1}}~{}~{}\alpha||\mathcal{X}-\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{2}\bm{U}_{2}^{\top}\times_{3}\bm{U}_{3}^{\top}||_{F}^{2}+(1-\alpha)||\bm{y}-\bm{1}_{M}\cdot\tilde{\beta}_{0}-\bm{S}_{(4)}\cdot\tilde{\bm{\beta}}_{1}||_{F}^{2},$

which is convex. Thus, it can be solved by setting the derivatives to be
zeros–that is $\frac{d\Psi}{d\bm{U}_{1}}=\bm{0}$, where
$\Psi=\alpha||\mathcal{X}-\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{2}\bm{U}_{2}^{\top}\times_{3}\bm{U}_{3}^{\top}||_{F}^{2}+(1-\alpha)||\bm{y}-\bm{1}_{M}\cdot\tilde{\beta}_{0}-\bm{S}_{(4)}\cdot\tilde{\bm{\beta}}_{1}||_{F}^{2}$.
This implies
$\frac{d}{d\bm{U}_{1}}(||\mathcal{X}-\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{2}\bm{U}_{2}^{\top}\times_{3}\bm{U}_{3}^{\top}||_{F}^{2})=\bm{0}$.
According to the communication law of tensor mode multiplication, we have
$\frac{d}{d\bm{U}_{1}}(||\mathcal{X}-(\mathcal{S}\times_{2}\bm{U}_{2}^{\top}\times_{3}\bm{U}_{3}^{\top})\times_{1}\bm{U}_{1}^{\top}||_{F}^{2})=\bm{0}$.
Thus,
$\frac{d}{d\bm{U}_{1}}(||\mathcal{X}-\mathcal{S}_{{U}_{1}}\times_{1}\bm{U}_{1}^{\top}||_{F}^{2})=\bm{0}$,
where
$\mathcal{S}_{{U}_{1}}=\mathcal{S}\times_{2}\bm{U}_{2}^{\top}\times_{3}\bm{U}_{3}^{\top}$.
Furthermore, we have
$\frac{d}{d\bm{U}_{1}}(||\bm{X}_{(1)}-\bm{U}_{1}^{\top}\cdot\bm{S}_{{U}_{1}(1)}||_{F}^{2})=\bm{0}$
due to the fact that $||\mathcal{S}||_{F}^{2}=||\bm{S}_{(n)}||_{F}^{2}$ and
the property of tensor mode multiplication
$\mathcal{S}\times_{n}\bm{U}=\bm{U}\cdot\bm{S}_{(n)}$. By taking the
derivative of the Frobenius norm, we have
$2(\bm{X}_{(1)}-\bm{U}_{1}^{\top}\cdot\bm{S}_{{U}_{1}(1)})\cdot(-\bm{S}_{{U}_{1}(1)}^{\top})=\bm{0}$.
Thus,
$\bm{U}_{1}^{\top}\cdot\bm{S}_{{U}_{1}(1)}\cdot\bm{S}_{{U}_{1}(1)}^{\top}=\bm{X}_{(1)}\cdot\bm{S}_{{U}_{1}(1)}^{\top}$,
which gives that
$\bm{U}_{1}^{\top}=\bm{X}_{(1)}\cdot\bm{S}_{{U}_{1}(1)}^{\top}\cdot(\bm{S}_{{U}_{1}(1)}\cdot\bm{S}_{{U}_{1}(1)}^{\top})^{-1}$.
Finally, we have
$\bm{U}_{1}=(\bm{X}_{(1)}\cdot\bm{S}_{{U}_{1}(1)}^{\top}\cdot(\bm{S}_{{U}_{1}(1)}\cdot\bm{S}_{{U}_{1}(1)}^{\top})^{-1})^{\top}$.

##  2 Proof of Proposition 2

The original optimization problem is

$\arg\min_{\bm{U}_{2}}~{}~{}\alpha||\mathcal{P}_{\Omega}(\mathcal{X}-\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{2}\bm{U}_{2}^{\top}\times_{3}\bm{U}_{3}^{\top})||_{F}^{2}+(1-\alpha)||\bm{y}-\bm{1}_{M}\cdot\tilde{\beta}_{0}-\bm{S}_{(4)}\cdot\tilde{\bm{\beta}}_{1}||_{F}^{2}\mbox{},$

which is equivalent to the following problem when data is complete:

$\arg\min_{\bm{U}_{2}}~{}~{}\alpha||\mathcal{X}-\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{2}\bm{U}_{2}^{\top}\times_{3}\bm{U}_{3}^{\top}||_{F}^{2}+(1-\alpha)||\bm{y}-\bm{1}_{M}\cdot\tilde{\beta}_{0}-\bm{S}_{(4)}\cdot\tilde{\bm{\beta}}_{1}||_{F}^{2}\mbox{}$

which is convex. Therefore, it can be solved by setting the derivatives to be
zeros–that is $\frac{d\Psi}{d\bm{U}_{2}}=\bm{0}$, where
$\Psi=\alpha||\mathcal{X}-\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{2}\bm{U}_{2}^{\top}\times_{3}\bm{U}_{3}^{\top}||_{F}^{2}+(1-\alpha)||\bm{y}-\bm{1}_{M}\cdot\tilde{\beta}_{0}-\bm{S}_{(4)}\cdot\tilde{\bm{\beta}}_{1}||_{F}^{2}$.
This implies
$\frac{d}{d\bm{U}_{2}}(||\mathcal{X}-\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{2}\bm{U}_{2}^{\top}\times_{3}\bm{U}_{3}^{\top}||_{F}^{2})=\bm{0}$.
According to the communication law of tensor mode multiplication, we have
$\frac{d}{d\bm{U}_{2}}(||\mathcal{X}-(\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{3}\bm{U}_{3}^{\top})\times_{2}\bm{U}_{2}^{\top}||_{F}^{2})=\bm{0}$.
Thus,
$\frac{d}{d\bm{U}_{2}}(||\mathcal{X}-\mathcal{S}_{{U}_{2}}\times_{2}\bm{U}_{2}^{\top}||_{F}^{2})=\bm{0}$,
where
$\mathcal{S}_{{U}_{2}}=\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{3}\bm{U}_{3}^{\top}$.
Furthermore, we have
$\frac{d}{d\bm{U}_{2}}(||\bm{X}_{(2)}-\bm{U}_{2}^{\top}\cdot\bm{S}_{{U}_{2}(2)}||_{F}^{2})=\bm{0}$
due to the fact that $||\mathcal{S}||_{F}^{2}=||\bm{S}_{(n)}||_{F}^{2}$ and
the property of tensor mode multiplicaion
$\mathcal{S}\times_{n}\bm{U}=\bm{U}\cdot\bm{S}_{(n)}$. By taking the
derivative of the Frobenius norm, we have
$2(\bm{X}_{(2)}-\bm{U}_{2}^{\top}\cdot\bm{S}_{{U}_{2}(2)})\cdot(-\bm{S}_{{U}_{2}(2)}^{\top})=\bm{0}$.
Thus,
$\bm{U}_{2}^{\top}\cdot\bm{S}_{{U}_{2}(2)}\cdot\bm{S}_{{U}_{2}(2)}^{\top}=\bm{X}_{(2)}\cdot\bm{S}_{{U}_{2}(2)}^{\top}$
which gives that
$\bm{U}_{2}^{\top}=\bm{X}_{(2)}\cdot\bm{S}_{{U}_{2}(2)}^{\top}\cdot(\bm{S}_{{U}_{2}(2)}\cdot\bm{S}_{{U}_{2}(2)}^{\top})^{-1}$.
Finally we have
$\bm{U}_{2}=(\bm{X}_{(2)}\cdot\bm{S}_{{U}_{2}(2)}^{\top}\cdot(\bm{S}_{{U}_{2}(2)}\cdot\bm{S}_{{U}_{2}(2)}^{\top})^{-1})^{\top}$.

##  3 Proof of Proposition 3

The original optimization problem is

$\arg\min_{\bm{U}_{3}}~{}~{}\alpha||\mathcal{P}_{\Omega}(\mathcal{X}-\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{2}\bm{U}_{2}^{\top}\times_{3}\bm{U}_{3}^{\top})||_{F}^{2}+(1-\alpha)||\bm{y}-\bm{1}_{M}\cdot\tilde{\beta}_{0}-\bm{S}_{(4)}\cdot\tilde{\bm{\beta}}_{1}||_{F}^{2},$

which is equivalent to the following problem when data is complete:

$\arg\min_{\bm{U}_{3}}~{}~{}\alpha||\mathcal{X}-\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{2}\bm{U}_{2}^{\top}\times_{3}\bm{U}_{3}^{\top}||_{F}^{2}+(1-\alpha)||\bm{y}-\bm{1}_{M}\cdot\tilde{\beta}_{0}-\bm{S}_{(4)}\cdot\tilde{\bm{\beta}}_{1}||_{F}^{2}$

which is convex. Therefore, it can be solved by setting the derivatives to be
zeros–that is $\frac{d\Psi}{d\bm{U}_{3}}=\bm{0}$, where
$\Psi=\alpha||\mathcal{X}-\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{2}\bm{U}_{2}^{\top}\times_{3}\bm{U}_{3}^{\top}||_{F}^{2}+(1-\alpha)||\bm{y}-\bm{1}_{M}\cdot\tilde{\beta}_{0}-\bm{S}_{(4)}\cdot\tilde{\bm{\beta}}_{1}||_{F}^{2}$.
This implies
$\frac{d}{d\bm{U}_{3}}(||\mathcal{X}-\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{2}\bm{U}_{2}^{\top}\times_{3}\bm{U}_{3}^{\top}||_{F}^{2})=\bm{0}$.
According to the communication law of tensor mode multiplication, we have
$\frac{d}{d\bm{U}_{3}}(||\mathcal{X}-(\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{2}\bm{U}_{2}^{\top})\times_{3}\bm{U}_{3}^{\top}||_{F}^{2})=\bm{0}$.
Thus,
$\frac{d}{d\bm{U}_{3}}(||\mathcal{X}-\mathcal{S}_{{U}_{3}}\times_{3}\bm{U}_{3}^{\top}||_{F}^{2})=\bm{0}$,
where
$\mathcal{S}_{{U}_{3}}=\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{2}\bm{U}_{2}^{\top}$.
Furthermore, we have
$\frac{d}{d\bm{U}_{3}}(||\bm{X}_{(3)}-\bm{U}_{3}^{\top}\cdot\bm{S}_{{U}_{3}(3)}||_{F}^{2})=\bm{0}$
due to the fact that $||\mathcal{S}||_{F}^{2}=||\bm{S}_{(n)}||_{F}^{2}$ and
the property of tensor mode multiplicaion
$\mathcal{S}\times_{n}\bm{U}=\bm{U}\cdot\bm{S}_{(n)}$. By taking the
derivative of the Frobenius norm, we have
$2(\bm{X}_{(3)}-\bm{U}_{3}^{\top}\cdot\bm{S}_{{U}_{3}(3)})\cdot(-\bm{S}_{{U}_{3}(3)}^{\top})=\bm{0}$.
Thus,
$\bm{U}_{3}^{\top}\cdot\bm{S}_{{U}_{3}(3)}\cdot\bm{S}_{{U}_{3}(3)}^{\top}=\bm{X}_{(3)}\cdot\bm{S}_{{U}_{3}(3)}^{\top}$,
which gives that
$\bm{U}_{3}^{\top}=\bm{X}_{(3)}\cdot\bm{S}_{{U}_{3}(3)}^{\top}\cdot(\bm{S}_{{U}_{3}(3)}\cdot\bm{S}_{{U}_{3}(3)}^{\top})^{-1}$.
Finally, we have
$\bm{U}_{3}=(\bm{X}_{(3)}\cdot\bm{S}_{{U}_{3}(3)}^{\top}\cdot(\bm{S}_{{U}_{3}(3)}\cdot\bm{S}_{{U}_{3}(3)}^{\top})^{-1})^{\top}$.

##  4 Proof of Proposition 4

The original optimization problem is

$\arg\min_{S}~{}~{}\alpha||\mathcal{P}_{\Omega}(\mathcal{X}-\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{2}\bm{U}_{2}^{\top}\times_{3}\bm{U}_{3}^{\top})||_{F}^{2}+(1-\alpha)||\bm{y}-\bm{1}_{M}\cdot\tilde{\beta}_{0}-\bm{S}_{(4)}\cdot\tilde{\bm{\beta}}_{1}||_{F}^{2},$

which is equivalent to the following problem when data is complete:

$\arg\min_{\hat{S}}~{}~{}\alpha||\mathcal{X}-\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{2}\bm{U}_{2}^{\top}\times_{3}\bm{U}_{3}^{\top}||_{F}^{2}+(1-\alpha)||\bm{y}-\bm{1}_{M}\cdot\tilde{\beta}_{0}-\bm{S}_{(4)}\cdot\tilde{\bm{\beta}}_{1}||_{F}^{2},$

which is convex. Thus, it can be solved by setting the derivatives to be
zeros–that is $\frac{d\Psi}{d{\bm{S}}}=\bm{0}$, where
$\Psi=\alpha||\mathcal{X}-\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{2}\bm{U}_{2}^{\top}\times_{3}\bm{U}_{3}^{\top}||_{F}^{2}+(1-\alpha)||\bm{y}-\bm{1}_{M}\cdot\tilde{\beta}_{0}-\bm{S}_{(4)}\cdot\tilde{\bm{\beta}}_{1}||_{F}^{2}$.
According to the connection between Kronecker product and tensor mode
multiplication (Kolda, 2006), we have
$\frac{d}{d{\bm{S}}}(\alpha||\bm{X}_{(4)}-\bm{S}_{(4)}\cdot(\bm{U}_{3}\otimes\bm{U}_{2}\otimes\bm{U}_{1})||_{F}^{2}+(1-\alpha)||\bm{y}-\bm{1}_{M}\cdot\tilde{\beta}_{0}-\bm{S}_{(4)}\cdot\tilde{\bm{\beta}}_{1}||_{F}^{2})=\bm{0}$.
By taking the derivative of the Frobenius norm, we have
$2\alpha\cdot[\bm{X}_{(4)}-\bm{S}_{(4)}\cdot(\bm{U}_{3}\otimes\bm{U}_{2}\otimes\bm{U}_{1})]\cdot[-(\bm{U}_{3}\otimes\bm{U}_{2}\otimes\bm{U}_{1})^{\top}]+2(1-\alpha)\cdot[\bm{y}-\bm{1}_{M}\cdot\tilde{\beta}_{0}-\bm{S}_{(4)}\cdot\tilde{\bm{\beta}}_{1}]\cdot(-\tilde{\bm{\beta}}_{1}^{\top})=\bm{0}$.
Thus,
$-2\alpha\cdot\bm{X}_{(4)}\cdot(\bm{U}_{3}\otimes\bm{U}_{2}\otimes\bm{U}_{1})^{\top}+2\alpha\cdot\bm{S}_{(4)}\cdot(\bm{U}_{3}\otimes\bm{U}_{2}\otimes\bm{U}_{1})\cdot(\bm{U}_{3}\otimes\bm{U}_{2}\otimes\bm{U}_{1})^{\top}+2(1-\alpha)(\bm{y}-\bm{1}_{M}\cdot\beta_{0})\cdot(-\tilde{\bm{\beta}}_{1}^{\top})+2(1-\alpha)\cdot(\bm{S}_{(4)}\cdot\tilde{\bm{\beta}}_{1}\cdot\tilde{\bm{\beta}}_{1}^{\top})=\bm{0}$.
Finally, we have
$\bm{S}_{(4)}=\big{[}\alpha\cdot\bm{X}_{(4)}\cdot(\bm{U}_{3}\otimes\bm{U}_{2}\otimes\bm{U}_{1})^{\top}+(1-\alpha)\cdot(\bm{y}-\bm{1}_{M}\cdot\tilde{\beta}_{0})\cdot\tilde{\bm{\beta}}_{1}^{\top}\big{]}\cdot\big{[}\alpha\cdot(\bm{U}_{3}\otimes\bm{U}_{2}\otimes\bm{U}_{1})\cdot(\bm{U}_{3}\otimes\bm{U}_{2}\otimes\bm{U}_{1})^{\top}+(1-\alpha)\cdot\tilde{\bm{\beta}}_{1}\cdot\tilde{\bm{\beta}}_{1}^{\top}\big{]}^{-1}$.

##  5 Proof of Lemma 1

Let $\bm{a}_{m}\in\mathbb{R}^{1\times N}$ denotes the $m$th row of matrix
$\bm{A}\in\mathbb{R}^{M\times N}$ and $\bm{b}_{m}\in\mathbb{R}^{1\times P}$
denotes the $m$th row of matrix $\bm{B}\in\mathbb{R}^{M\times P}$,
$m=1,\ldots,M$, then we have

$\bm{A}-\bm{B}\bm{C}=[(\bm{a}_{1}-\bm{b}_{1}\bm{C})^{\top},\ldots,(\bm{a}_{M}-\bm{b}_{M}\bm{C})^{\top}]^{\top}.$

Based on the definition of Frobenius norm, the original objective function in
Lemma 1 can be transformed as follows:

$\|\bm{A}-\bm{B}\bm{C}\|_{F}^{2}=\|[(\bm{a}_{1}-\bm{b}_{1}\bm{C})^{\top},\ldots,(\bm{a}_{M}-\bm{b}_{M}\bm{C})^{\top}]^{\top}\|_{F}^{2}=\sum_{m=1}^{M}\|\bm{a}_{m}-\bm{b}_{m}\bm{C}\|_{F}^{2}.$

Therefore, we have the following:

$\arg\min_{\bm{B}}\|\bm{A}-\bm{B}\bm{C}\|_{F}^{2}=\arg\min_{\\{\bm{b}_{m}\\}_{m=1}^{M}}\sum_{m=1}^{M}\|\bm{a}_{m}-\bm{b}_{m}\bm{C}\|_{F}^{2}.$

where ${\bm{B}=[\bm{b}_{1}^{\top},\ldots,\bm{b}_{M}^{\top}]^{\top}}$.
Therefore, to solve the original objective function, we can simply solve the
following $M$ sub problems:

$\arg\min_{\bm{b}_{m}}\|\bm{a}_{m}-\bm{b}_{m}\bm{C}\|_{F}^{2},\quad
m=1,\ldots,M.$

##  6 Proof of Proposition 5

The original optimization problem is

$\arg\min_{\bm{U}_{1}}~{}~{}\alpha||\mathcal{P}_{\Omega}(\mathcal{X}-\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{2}\bm{U}_{2}^{\top}\times_{3}\bm{U}_{3}^{\top})||_{F}^{2}+(1-\alpha)||\bm{y}-\bm{1}_{M}\cdot\tilde{\beta}_{0}-\bm{S}_{(4)}\cdot\tilde{\bm{\beta}}_{1}||_{F}^{2},$

which is equivalent to the following problem when data is missing:

$\arg\min_{\bm{U}_{1}}~{}~{}\alpha||\mathcal{X}-(\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{2}\bm{U}_{2}^{\top}\times_{3}\bm{U}_{3}^{\top})\odot\text{logic}(\mathcal{X})||_{F}^{2}+(1-\alpha)||\bm{y}-\bm{1}_{M}\cdot\tilde{\beta}_{0}-\bm{S}_{(4)}\cdot\tilde{\bm{\beta}}_{1}||_{F}^{2},$

where $\odot$ is the inner product, logic($\mathcal{X}$) denotes the logical
value of $\mathcal{X}$–that is–if an entry is observed, its logical value is
1; otherwise, it is 0. Since the problem is convex, it can be solved by
setting the derivatives to be zeros, i.e., $\frac{d\Psi}{d\bm{U}_{1}}=\bm{0}$,
where
$\Psi=\alpha||\mathcal{X}-(\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{2}\bm{U}_{2}^{\top}\times_{3}\bm{U}_{3}^{\top})\odot\text{logic}(\mathcal{X})||_{F}^{2}+(1-\alpha)||\bm{y}-\bm{1}_{M}\cdot\tilde{\beta}_{0}-\bm{S}_{(4)}\cdot\tilde{\bm{\beta}}_{1}||_{F}^{2}$.
This implies
$\frac{d}{d\bm{U}_{1}}(||\mathcal{X}-(\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{2}\bm{U}_{2}^{\top}\times_{3}\bm{U}_{3}^{\top})\odot\text{logic}(\mathcal{X})||_{F}^{2})=\bm{0}$.
According to the communication law of tensor mode multiplication, we have
$\frac{d}{d\bm{U}_{1}}(||\mathcal{X}-[(\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{2}\bm{U}_{2}^{\top})\times_{3}\bm{U}_{3}^{\top}]\odot\text{logic}(\mathcal{X})||_{F}^{2})=\bm{0}$.
Thus,
$\frac{d}{d\bm{U}_{1}}(||\mathcal{X}-(\mathcal{S}_{{U}_{1}}\times_{1}\bm{U}_{1}^{\top})\odot\text{logic}(\mathcal{X})||_{F}^{2})=\bm{0}$,
where
$\mathcal{S}_{{U}_{1}}=\mathcal{S}\times_{2}\bm{U}_{2}^{\top}\times_{3}\bm{U}_{3}^{\top}$.
Furthermore, we have
$\frac{d}{d\bm{U}_{1}}(||\bm{X}_{(1)}-(\bm{U}_{1}^{\top}\cdot\bm{S}_{{U}_{1}(1)})\odot\text{logic}(\bm{X}_{(1)})||_{F}^{2})=\bm{0}$
since $||\mathcal{S}||_{F}^{2}=||\bm{S}_{(n)}||_{F}^{2}$ and
$\mathcal{S}\times_{n}\bm{U}=\bm{U}\cdot\bm{S}_{(n)}$.

Figure 12 shows the pattern of mode-1 matricization of the 4D tensor
$\mathcal{X}$ when it has missing entries whose indices can be denoted by a
set $\Omega\subseteq\\{(i_{1},i_{2},i_{3},m),1\leq i_{1}\leq I_{1},1\leq
i_{2}\leq I_{2},1\leq i_{3}\leq I_{3},1\leq m\leq M\\}$. Based on Lemma 1, we
can sequentially optimize each column of $\bm{U}_{1}$.

Figure 12: An illustration of the data missing pattern in Proposition 5
(stripes representing available elements).

Specifically, we denote the $i_{1}$th row in $\bm{U}_{1}^{\top}$ as
$\bm{u}_{1}^{i_{1}\top}$ (blue solid row of $\bm{U}_{1}^{\top}$ in Figure 12).
The available elements in the $i_{1}$th row of $\bm{X}_{(1)}$ are denoted as
$\bm{x}_{(1)}^{i_{1},\pi_{i_{1}}}$ (blue striped squares of
$\bm{x}_{(1)}^{i_{1}}$ in Figure 12). In $\bm{S}_{U_{1}(1)}^{\pi_{i_{1}}}$, we
choose the columns whose indices are the same as those of the available
elements of $\bm{x}_{(1)}^{i_{1}}$ (blue striped columns of
$\bm{S}_{U_{1}(1)}^{\pi_{i_{1}}}$ in Figure 12). As a result, we have
$\frac{d}{d\bm{u}_{1}^{i_{1}}}(\Sigma_{{i_{1}}=1}^{I_{1}}||\bm{x}_{(1)}^{i_{1},\pi_{i_{1}}}-({\bm{u}_{1}^{i_{1}\top}}\cdot\bm{S}_{{U}_{1}(1)}^{\pi_{i_{1}}})||_{F}^{2})=\bm{0}$.
Because we only take the derivative of $\bm{u}_{1}^{i_{1}}$, we have
$\frac{d}{d\bm{u}_{1}^{i_{1}}}(||\bm{x}_{(1)}^{i_{1},\pi_{i_{1}}}-({\bm{u}_{1}^{i_{1}\top}}\cdot\bm{S}_{{U}_{1}(1)}^{\pi_{i_{1}}})||_{F}^{2})=\bm{0}$.
By taking the derivative of the Frobenius norm, we have
$2(\bm{x}_{(1)}^{i_{1},\pi_{i_{1}}}-\bm{u}_{1}^{i_{1}\top}\cdot\bm{S}_{{U}_{1}(1)}^{\pi_{i_{1}}})\cdot(-{\bm{S}_{{U}_{1}(1)}^{\pi_{i_{1}}}}^{\top})=\bm{0}$.
Thus,
$\bm{u}_{1}^{i_{1}\top}\cdot\bm{S}_{{U}_{1}(1)}^{\pi_{i_{1}}}\cdot{\bm{S}_{{U}_{1}(1)}^{\pi_{i_{1}}}}^{\top}=\bm{x}_{(1)}^{i_{1},\pi_{i_{1}}}\cdot{\bm{S}_{{U}_{1}(1)}^{\pi_{i_{1}}}}^{\top}$,
which gives that
$\bm{u}_{1}^{i_{1}}=(\bm{x}_{(1)}^{i_{1},\pi_{i_{1}}}\cdot{\bm{S}_{{U}_{1}(1)}^{\pi_{i_{1}}}}^{\top}\cdot(\bm{S}_{{U}_{1}(1)}^{\pi_{i_{1}}}\cdot{\bm{S}_{{U}_{1}(1)}^{\pi_{i_{1}}}}^{\top})^{-1})^{\top}$.

##  7 Proof of Proposition 6

The original optimization problem is

$\arg\min_{\bm{U}_{2}}~{}~{}\alpha||\mathcal{P}_{\Omega}(\mathcal{X}-\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{2}\bm{U}_{2}^{\top}\times_{3}\bm{U}_{3}^{\top})||_{F}^{2}+(1-\alpha)||\bm{y}-\bm{1}_{M}\cdot\tilde{\beta}_{0}-\bm{S}_{(4)}\cdot\tilde{\bm{\beta}}_{1}||_{F}^{2},$

which is equivalent to the following problem when data is missing:

$\arg\min_{\bm{U}_{2}}~{}~{}\alpha||\mathcal{X}-(\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{2}\bm{U}_{2}^{\top}\times_{3}\bm{U}_{3}^{\top})\odot\text{logic}(\mathcal{X})||_{F}^{2}+(1-\alpha)||\bm{y}-\bm{1}_{M}\cdot\tilde{\beta}_{0}-\bm{S}_{(4)}\cdot\tilde{\bm{\beta}}_{1}||_{F}^{2},$

where $\odot$ is the inner product, logic($\mathcal{X}$) denotes the logical
value of $\mathcal{X}$. Since the problem is convex, it can be solved by
setting the derivatives to be zeros, i.e., $\frac{d\Psi}{d\bm{U}_{2}}=\bm{0}$,
where
$\Psi=\alpha||\mathcal{X}-(\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{2}\bm{U}_{2}^{\top}\times_{3}\bm{U}_{3}^{\top})\odot\text{logic}(\mathcal{X})||_{F}^{2}+(1-\alpha)||\bm{y}-\bm{1}_{M}\cdot\tilde{\beta}_{0}-\bm{S}_{(4)}\cdot\tilde{\bm{\beta}}_{1}||_{F}^{2}$.
This implies
$\frac{d}{d\bm{U}_{2}}(||\mathcal{X}-(\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{2}\bm{U}_{2}^{\top}\times_{3}\bm{U}_{3}^{\top})\odot\text{logic}(\mathcal{X})||_{F}^{2})=\bm{0}$.
According to the communication law of tensor mode multiplication, we have
$\frac{d}{d\bm{U}_{2}}(||\mathcal{X}-[(\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{3}\bm{U}_{3}^{\top})\times_{2}\bm{U}_{2}^{\top}]\odot\text{logic}(\mathcal{X})||_{F}^{2})=\bm{0}$.
Thus,
$\frac{d}{d\bm{U}_{2}}(||\mathcal{X}-(\mathcal{S}_{{U}_{2}}\times_{2}\bm{U}_{2}^{\top})\odot\text{logic}(\mathcal{X})||_{F}^{2})=\bm{0}$,
where
$\mathcal{S}_{{U}_{2}}=\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{3}\bm{U}_{3}^{\top}$.
Furthermore, we have
$\frac{d}{d\bm{U}_{2}}(||\bm{X}_{(2)}-(\bm{U}_{2}^{\top}\cdot\bm{S}_{{U}_{2}(2)})\odot\text{logic}(\bm{X}_{(2)})||_{F}^{2})=\bm{0}$
since $||\mathcal{S}||_{F}^{2}=||\bm{S}_{(n)}||_{F}^{2}$ and
$\mathcal{S}\times_{n}\bm{U}=\bm{U}\cdot\bm{S}_{(n)}$.

Figure 13: An illustration of the data missing pattern in Proposition 6
(stripes representing available elements).

Figure 13 shows the pattern of mode-2 matricization of the 4D tensor
$\mathcal{X}$ when it has missing entries whose indices can be denoted by a
set $\Omega\subseteq\\{(i_{1},i_{2},i_{3},m),1\leq i_{1}\leq I_{1},1\leq
i_{2}\leq I_{2},1\leq i_{3}\leq I_{3},1\leq m\leq M\\}$. Based on Lemma 1, we
can sequentially optimize each column of $\bm{U}_{2}$.

Specifically, we denote the $i_{2}$th row in $\bm{U}_{2}^{\top}$ as
$\bm{u}_{2}^{i_{2}\top}$ (blue solid row of $\bm{U}_{2}^{\top}$ in Figure 13).
The available elements in the $i_{2}$th row of $\bm{X}_{(2)}$ are denoted as
$\bm{x}_{(2)}^{i_{2},\pi_{i_{2}}}$ (blue striped squares of
$\bm{x}_{(2)}^{i_{2}}$ in Figure 13). In $\bm{S}_{U_{2}(2)}^{\pi_{i_{2}}}$, we
choose the columns whose indices are the same as those of the available
elements in $\bm{x}_{(2)}^{i_{2}}$ (blue striped columns of
$\bm{S}_{U_{2}(2)}^{\pi_{i_{2}}}$ in Figure 13). As a result, we have
$\frac{d}{d\bm{u}_{2}^{i_{2}}}(\Sigma_{{i_{2}}=1}^{I_{2}}||\bm{x}_{(2)}^{i_{2},\pi_{i_{2}}}-({\bm{u}_{2}^{i_{2}\top}}\cdot\bm{S}_{{U}_{2}(2)}^{\pi_{i_{2}}})||_{F}^{2})=\bm{0}$.
Since on the derivative of $\bm{u}_{2}^{i_{2}}$ is taken, we have
$\frac{d}{d\bm{u}_{2}^{i_{2}}}(||\bm{x}_{(2)}^{i_{2},\pi_{i_{2}}}-({\bm{u}_{2}^{i_{2}\top}}\cdot\bm{S}_{{U}_{2}(2)}^{\pi_{i_{2}}})||_{F}^{2})=\bm{0}$.
By taking the derivative of the Frobenius norm, we have
$2(\bm{x}_{(2)}^{i_{2},\pi_{i_{2}}}-\bm{u}_{2}^{i_{2}\top}\cdot\bm{S}_{{U}_{2}(2)}^{\pi_{i_{2}}})\cdot(-{\bm{S}_{{U}_{2}(2)}^{\pi_{i_{2}}}}^{\top})=\bm{0}$.
Thus,
$\bm{u}_{2}^{i_{2}\top}\cdot\bm{S}_{{U}_{2}(2)}^{\pi_{i_{2}}}\cdot{\bm{S}_{{U}_{2}(2)}^{\pi_{i_{2}}}}^{\top}=\bm{x}_{(2)}^{i_{2},\pi_{i_{2}}}\cdot{\bm{S}_{{U}_{2}(2)}^{\pi_{i_{2}}}}^{\top}$,
which gives that
$\bm{u}_{2}^{i_{2}}=(\bm{x}_{(2)}^{i_{2},\pi_{i_{2}}}\cdot{\bm{S}_{{U}_{2}(2)}^{\pi_{i_{2}}}}^{\top}\cdot(\bm{S}_{{U}_{2}(2)}^{\pi_{i_{2}}}\cdot{\bm{S}_{{U}_{2}(2)}^{\pi_{i_{2}}}}^{\top})^{-1})^{\top}$.

##  8 Proof of Proposition 7

The original optimization problem is

$\arg\min_{\bm{U}_{3}}~{}~{}\alpha||\mathcal{P}_{\Omega}(\mathcal{X}-\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{2}\bm{U}_{2}^{\top}\times_{3}\bm{U}_{3}^{\top})||_{F}^{2}+(1-\alpha)||\bm{y}-\bm{1}_{M}\cdot\tilde{\beta}_{0}-\bm{S}_{(4)}\cdot\tilde{\bm{\beta}}_{1}||_{F}^{2},$

which is equivalent to the following problem when data is missing:

$\arg\min_{\bm{U}_{3}}~{}~{}\alpha||\mathcal{X}-(\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{2}\bm{U}_{2}^{\top}\times_{3}\bm{U}_{3}^{\top})\odot\text{logic}(\mathcal{X})||_{F}^{2}+(1-\alpha)||\bm{y}-\bm{1}_{M}\cdot\tilde{\beta}_{0}-\bm{S}_{(4)}\cdot\tilde{\bm{\beta}}_{1}||_{F}^{2},$

where $\odot$ is the inner product, logic($\mathcal{X}$) denotes the logical
value of $\mathcal{X}$. Since the problem is convex, it can be solved by
setting the derivatives to be zeros, i.e., $\frac{d\Psi}{d\bm{U}_{3}}=\bm{0}$,
where
$\Psi=\alpha||\mathcal{X}-(\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{2}\bm{U}_{2}^{\top}\times_{3}\bm{U}_{3}^{\top})\odot\text{logic}(\mathcal{X})||_{F}^{2}+(1-\alpha)||\bm{y}-\bm{1}_{M}\cdot\tilde{\beta}_{0}-\bm{S}_{(4)}\cdot\tilde{\bm{\beta}}_{1}||_{F}^{2}$.
This implies
$\frac{d}{d\bm{U}_{3}}(||\mathcal{X}-(\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{2}\bm{U}_{2}^{\top}\times_{3}\bm{U}_{3}^{\top})\odot\text{logic}(\mathcal{X})||_{F}^{2})=\bm{0}$.
According to the communication law of tensor mode multiplication, we have
$\frac{d}{d\bm{U}_{3}}(||\mathcal{X}-[(\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{2}\bm{U}_{2}^{\top})\times_{3}\bm{U}_{3}^{\top}]\odot\text{logic}(\mathcal{X})||_{F}^{2})=\bm{0}$.
Thus,
$\frac{d}{d\bm{U}_{3}}(||\mathcal{X}-(\mathcal{S}_{{U}_{3}}\times_{3}\bm{U}_{3}^{\top})\odot\text{logic}(\mathcal{X})||_{F}^{2})=\bm{0}$,
where
$\mathcal{S}_{{U}_{3}}=\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{2}\bm{U}_{2}^{\top}$.
Furthermore, we have
$\frac{d}{d\bm{U}_{3}}(||\bm{X}_{(3)}-(\bm{U}_{3}^{\top}\cdot\bm{S}_{{U}_{3}(3)})\odot\text{logic}(\bm{X}_{(3)})||_{F}^{2})=\bm{0}$
since $||\mathcal{S}||_{F}^{2}=||\bm{S}_{(n)}||_{F}^{2}$ and
$\mathcal{S}\times_{n}\bm{U}=\bm{U}\cdot\bm{S}_{(n)}$.

Figure 14: An illustration of the data missing pattern in Proposition 7
(stripes representing available elements).

Figure 14 shows the pattern of mode-3 matricization of the 4D tensor
$\mathcal{X}$ when it has missing entries whose indices can be denoted by a
set $\Omega\subseteq\\{(i_{1},i_{2},i_{3},m),1\leq i_{1}\leq I_{1},1\leq
i_{2}\leq I_{2},1\leq i_{3}\leq I_{3},1\leq m\leq M\\}$. Based on Lemma 1, we
can sequentially optimize each column of $\bm{U}_{3}$. The $i_{3}$th row in
$\bm{U}_{3}^{\top}$ is denoted as $\bm{u}_{3}^{i_{3}\top}$ (blue solid row of
$\bm{U}_{3}^{\top}$ in Figure 14). The available elements in the $i_{3}$th row
of $\bm{X}_{(3)}$ are denoted as $\bm{x}_{(3)}^{i_{3},\pi_{i_{3}}}$ (blue
striped squares of $\bm{x}_{(3)}^{i_{3}}$ in Figure 14). In
$\bm{S}_{U_{3}(3)}^{\pi_{i_{3}}}$, we choose the columns whose indices are the
same as those of the available elements of $\bm{x}_{(3)}^{i_{3}}$ (blue
striped columns of $\bm{S}_{U_{3}(3)}^{\pi_{i_{3}}}$ in Figure 14). Thus, we
have
$\frac{d}{d\bm{u}_{3}^{i_{3}}}(\Sigma_{{i_{3}}=1}^{I_{3}}||\bm{x}_{(3)}^{i_{3},\pi_{i_{3}}}-({\bm{u}_{3}^{i_{3}\top}}\cdot\bm{S}_{{U}_{3}(3)}^{\pi_{i_{3}}})||_{F}^{2})=\bm{0}$,
which yields
$\frac{d}{d\bm{u}_{3}^{i_{3}}}(||\bm{x}_{(3)}^{i_{3},\pi_{i_{3}}}-({\bm{u}_{3}^{i_{3}\top}}\cdot\bm{S}_{{U}_{3}(3)}^{\pi_{i_{3}}})||_{F}^{2})=\bm{0}$.
By taking the derivative of the Frobenius norm, we have
$2(\bm{x}_{(3)}^{i_{3},\pi_{i_{3}}}-\bm{u}_{3}^{i_{3}\top}\cdot\bm{S}_{{U}_{3}(3)}^{\pi_{i_{3}}})\cdot(-{\bm{S}_{{U}_{3}(3)}^{\pi_{i_{3}}}}^{\top})=\bm{0}$.
Thus,
$\bm{u}_{3}^{i_{3}\top}\cdot\bm{S}_{{U}_{3}(3)}^{\pi_{i_{3}}}\cdot{\bm{S}_{{U}_{3}(3)}^{\pi_{i_{3}}}}^{\top}=\bm{x}_{(3)}^{i_{3},\pi_{i_{3}}}\cdot{\bm{S}_{{U}_{3}(3)}^{\pi_{i_{3}}}}^{\top}$,
which gives that
$\bm{u}_{3}^{i_{3}}=(\bm{x}_{(3)}^{i_{3},\pi_{i_{3}}}\cdot{\bm{S}_{{U}_{3}(3)}^{\pi_{i_{3}}}}^{\top}\cdot(\bm{S}_{{U}_{3}(3)}^{\pi_{i_{3}}}\cdot{\bm{S}_{{U}_{3}(3)}^{\pi_{i_{3}}}}^{\top})^{-1})^{\top}$.

##  9 Proof of Proposition 8

The original optimization problem is

$\arg\min_{S}~{}~{}\alpha||\mathcal{P}_{\Omega}(\mathcal{X}-\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{2}\bm{U}_{2}^{\top}\times_{3}\bm{U}_{3}^{\top})||_{F}^{2}+(1-\alpha)||\bm{y}-\bm{1}_{M}\cdot\tilde{\beta}_{0}-\bm{S}_{(4)}\cdot\tilde{\bm{\beta}}_{1}||_{F}^{2}$

which is equivalent to the following problem when data is missing:

$\arg\min_{S}~{}~{}\alpha||\mathcal{X}-(\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{2}\bm{U}_{2}^{\top}\times_{3}\bm{U}_{3}^{\top})\odot\text{logic}(\mathcal{X})||_{F}^{2}+(1-\alpha)||\bm{y}-\bm{1}_{M}\cdot\tilde{\beta}_{0}-\bm{S}_{(4)}\cdot\tilde{\bm{\beta}}_{1}||_{F}^{2}$

where $\odot$ is the inner product, logic($\mathcal{X}$) denotes the logical
value of $\mathcal{X}$. Since the problem is convex, it can be solved by
setting the derivative to be zeros–that is $\frac{d\Psi}{d\bm{S}}=\bm{0}$,
where
$\Psi=\alpha||\mathcal{X}-(\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{2}\bm{U}_{2}^{\top}\times_{3}\bm{U}_{3}^{\top})\odot\text{logic}(\mathcal{X})||_{F}^{2}+(1-\alpha)||\bm{y}-\bm{1}_{M}\cdot\tilde{\beta}_{0}-\bm{S}_{(4)}\cdot\tilde{\bm{\beta}}_{1}||_{F}^{2}$.
According to the connection between Kronecker product and tensor mode
multiplication, we have
$\frac{d}{d\bm{S}}(\alpha||\bm{X}_{(4)}-[\bm{S}_{(4)}\cdot(\bm{U}_{3}\otimes\bm{U}_{2}\otimes\bm{U}_{1})]\odot\text{logic}(\bm{X}_{(4)})||_{F}^{2}+(1-\alpha)||\bm{y}-\bm{1}_{M}\odot\tilde{\beta}_{0}-\bm{S}_{(4)}\odot\tilde{\bm{\beta}}_{1}||_{F}^{2})=\bm{0}$.

Figure 15: An illustration of the data missing pattern in Proposition 8
(stripes representing available elements).

Figure 15 shows the pattern of mode-4 matricization of the 4D tensor
$\mathcal{X}$ when it has missing entries whose indices can be denoted by a
set $\Omega\subseteq\\{(i_{1},i_{2},i_{3},m),1\leq i_{1}\leq I_{1},1\leq
i_{2}\leq I_{2},1\leq i_{3}\leq I_{3},1\leq m\leq M\\}$. Based on Lemma 1, we
can sequentially optimize each column of $\bm{U}_{3}$.

The $m$th row in $\bm{S}_{(4)}$ is denoted as $\bm{s}_{(4)}^{m}$ (blue solid
row of $\bm{S}_{(4)}$ in Figure 15). The available elements in the $m$th row
of $\bm{X}_{(4)}$ are denoted as $\bm{x}_{(4)}^{m,\pi_{m}}$ (blue striped
squares of $\bm{x}_{(4)}^{m}$ in Figure 15). In
$(\bm{U}_{3}\otimes\bm{U}_{2}\otimes\bm{U}_{1})$, we choose the columns whose
indices are the same as those of the available elements of $\bm{x}_{(4)}^{m}$
(blue striped columns of
$(\bm{U}_{3}\otimes\bm{U}_{2}\otimes\bm{U}_{1})^{\pi_{m}}$ in Figure 15).
Thus, we have
$\frac{d}{d\bm{s}_{(4)}^{m}}(\alpha\\{\Sigma_{m=1}^{M}||\bm{x}_{(4)}^{m,\pi_{m}}-[\bm{s}_{(4)}^{m}\cdot(\bm{U}_{3}\otimes\bm{U}_{2}\otimes\bm{U}_{1})^{\pi_{m}}]\\}||_{F}^{2}+(1-\alpha)\Sigma_{m=1}^{M}||y_{m}-\tilde{\beta}_{0}-\bm{s}_{(4)}^{m}\cdot\tilde{\bm{\beta}}_{1}||_{F}^{2})=\bm{0}$,
which yields
$\frac{d}{d\bm{s}_{(4)}^{m}}(\alpha\\{||\bm{x}_{(4)}^{m,\pi_{m}}-[\bm{s}_{(4)}^{m}\cdot(\bm{U}_{3}\otimes\bm{U}_{2}\otimes\bm{U}_{1})^{\pi_{m}}]\\}||_{F}^{2}+(1-\alpha)||y_{m}-\tilde{\beta}_{0}-\bm{s}_{(4)}^{m}\cdot\tilde{\bm{\beta}}_{1}||_{F}^{2})=\bm{0}$.
By taking the derivative of Frobenius norm, we have
$2\alpha\cdot[\bm{x}_{(4)}^{m,\pi_{m}}-\bm{s}_{(4)}^{m}\cdot(\bm{U}_{3}\otimes\bm{U}_{2}\otimes\bm{U}_{1})^{\pi_{m}}]\cdot[-(\bm{U}_{3}\otimes\bm{U}_{2}\otimes\bm{U}_{1})^{\pi_{m}\top}]+2(1-\alpha)\cdot[y_{m}-\tilde{\beta}_{0}-\bm{s}_{(4)}^{m}\cdot\tilde{\bm{\beta}}_{1}]\cdot(-\tilde{\bm{\beta}}_{1}^{\top})=\bm{0}$.
Thus,
$-2\alpha\cdot\bm{x}_{(4)}^{m,\pi_{m}}\cdot(\bm{U}_{3}\otimes\bm{U}_{2}\otimes\bm{U}_{1})^{\pi_{m}\top}+2\alpha\cdot\bm{s}_{(4)}^{m}\cdot(\bm{U}_{3}\otimes\bm{U}_{2}\otimes\bm{U}_{1})^{\pi_{m}}\cdot(\bm{U}_{3}\otimes\bm{U}_{2}\otimes\bm{U}_{1})^{\pi_{m}\top}+2(1-\alpha)(y_{m}-\tilde{\beta}_{0})\cdot(-\tilde{\bm{\beta}}_{1}^{\top})+2(1-\alpha)\cdot(\bm{s}_{(4)}^{m}\cdot\tilde{\bm{\beta}}_{1}\cdot\tilde{\bm{\beta}}_{1}^{\top})=\bm{0}$,
which gives that
$\bm{s}_{(4)}^{m}=\big{[}\alpha\cdot\bm{x}_{(4)}^{m,\pi_{m}}\cdot(\bm{U}_{3}\otimes\bm{U}_{2}\otimes\bm{U}_{1})^{\pi_{m}\top}+(1-\alpha)\cdot(y_{m}-\tilde{\beta}_{0})\cdot\tilde{\bm{\beta}}_{1}^{\top}\big{]}\cdot\big{[}\alpha\cdot(\bm{U}_{3}\otimes\bm{U}_{2}\otimes\bm{U}_{1})^{\pi_{m}}\cdot(\bm{U}_{3}\otimes\bm{U}_{2}\otimes\bm{U}_{1})^{\pi_{m}\top}+(1-\alpha)\cdot\tilde{\bm{\beta}}_{1}\cdot\tilde{\bm{\beta}}_{1}^{\top}\big{]}^{-1}$.

##  10 Proof of Proposition 9

The original optimization problem is

$\arg\min_{\bm{U}_{1}}~{}~{}\alpha||\mathcal{P}_{\Omega}(\mathcal{X}-\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{2}\bm{U}_{2}^{\top}\times_{3}\bm{U}_{3}^{\top})||_{F}^{2}+(1-\alpha)||\bm{y}-\bm{1}_{M}\cdot\tilde{\beta}_{0}-\bm{S}_{(4)}\cdot\tilde{\bm{\beta}}_{1}||_{F}^{2},$

which is equivalent to the following problem when data is incomplete:

$\arg\min_{\bm{U}_{1}}~{}~{}\alpha||\mathcal{X}-(\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{2}\bm{U}_{2}^{\top}\times_{3}\bm{U}_{3}^{\top})\odot\text{logic}(\mathcal{X})||_{F}^{2}+(1-\alpha)||\bm{y}-\bm{1}_{M}\cdot\tilde{\beta}_{0}-\bm{S}_{(4)}\cdot\tilde{\bm{\beta}}_{1}||_{F}^{2},$

where $\odot$ is the inner product, logic($\mathcal{X}$) denotes the logical
value of $\mathcal{X}$. Since the problem is convex, it can be solved by
setting the derivatives to be zeros, i.e.,
$\frac{d\Psi}{d\bm{U}_{1}}=\mathbf{0}$, where
$\Psi=\alpha||\mathcal{X}-(\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{2}\bm{U}_{2}^{\top}\times_{3}\bm{U}_{3}^{\top})\odot\text{logic}(\mathcal{X})||_{F}^{2}+(1-\alpha)||\bm{y}-\bm{1}_{M}\cdot\tilde{\beta}_{0}-\bm{S}_{(4)}\cdot\tilde{\bm{\beta}}_{1}||_{F}^{2}$.
This implies
$\frac{d}{d\bm{U}_{1}}(||\mathcal{X}-(\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{2}\bm{U}_{2}^{\top}\times_{3}\bm{U}_{3}^{\top})\odot\text{logic}(\mathcal{X})||_{F}^{2})=\bm{0}$.
According to the communication law of tensor mode multiplication, we have
$\frac{d}{d\bm{U}_{1}}(||\mathcal{X}-[(\mathcal{S}\times_{2}\bm{U}_{2}^{\top}\times_{3}\bm{U}_{3}^{\top})\times_{1}\bm{U}_{1}^{\top}]\odot\text{logic}(\mathcal{X})||_{F}^{2})=\bm{0}$.
Thus,
$\frac{d}{d\bm{U}_{1}}(||\mathcal{X}-(\mathcal{S}_{{U}_{1}}\times_{1}\bm{U}_{1}^{\top})\odot\text{logic}(\mathcal{X})||_{F}^{2})=\bm{0}$,
where
$\mathcal{S}_{{U}_{1}}=\mathcal{S}\times_{2}\bm{U}_{2}^{\top}\times_{3}\bm{U}_{3}^{\top}$.
Furthermore, we have
$\frac{d}{d\bm{U}_{1}}(||\bm{X}_{(1)}-(\bm{U}_{1}^{\top}\cdot\bm{S}_{{U}_{1}(1)})\odot\text{logic}(\bm{X}_{(1)})||_{F}^{2})=\bm{0}$
due to the fact that $||\mathcal{S}||_{F}^{2}=||\bm{S}_{(n)}||_{F}^{2}$ and
$\mathcal{S}\times_{n}\bm{U}=\bm{U}\cdot\bm{S}_{(n)}$ (a property of tensor
mode multiplication).

Figure 16: An illustration of the data missing pattern in Proposition 9
(stripes representing available columns).

As discussed earlier, for applications with missing images, the indices of
tensor $\mathcal{X}$’s missing entries can be denoted as
$\Omega\subseteq\\{(:,:,i_{3},m)$, $1\leq i_{3}\leq I_{3},1\leq m\leq M\\}$,
where “$:$” denotes all the indices in a dimension. As a result, it can be
easily shown that $\mathcal{X}$’s mode-1 matricization $\bm{X}_{(1)}$ has
missing columns (see Figure 16 for an illustration). Let $\pi$ be the set
consisting of the indices of available columns in $\bm{X}_{(1)}$, then we need
to solve
$\frac{d}{d\bm{U}_{1}}(||\bm{X}_{(1)}^{\pi}-\bm{U}_{1}^{\top}\cdot\bm{S}_{{U}_{1}(1)}^{\pi}||_{F}^{2})=\bm{0}$,
where $\bm{S}_{{U}_{1}(1)}^{\pi}$ denotes a matrix constituting the $\pi$
columns of $\bm{S}_{{U}_{1}(1)}$. Thus, we have
$2(\bm{X}_{(1)}^{\pi}-\bm{U}_{1}^{\top}\cdot\bm{S}_{{U}_{1}(1)}^{\pi})\cdot(-{\bm{S}_{{U}_{1}(1)}^{\pi}}^{\top})=\bm{0}$.
Thus,
$\bm{U}_{1}^{\top}\cdot\bm{S}_{{U}_{1}(1)}^{\pi}\cdot{\bm{S}_{{U}_{1}(1)}^{\pi}}^{\top}=\bm{X}_{(1)}^{\pi}\cdot{\bm{S}_{{U}_{1}(1)}^{\pi}}^{\top}$,
which gives that
$\bm{U}_{1}^{\top}=\bm{X}_{(1)}^{\pi}\cdot{\bm{S}_{{U}_{1}(1)}^{\pi}}^{\top}\cdot(\bm{S}_{{U}_{1}(1)}^{\pi}\cdot{\bm{S}_{{U}_{1}(1)}^{\pi}}^{\top})^{-1}$.
This yields the solution
$\bm{U}_{1}=(\bm{X}_{(1)}^{\pi}\cdot{\bm{S}_{{U}_{1}(1)}^{\pi}}^{\top}\cdot(\bm{S}_{{U}_{1}(1)}^{\pi}\cdot{\bm{S}_{{U}_{1}(1)}^{\pi}}^{\top})^{-1})^{\top}$.

##  11 Proof of Proposition 10

The original optimization problem is

$\arg\min_{\bm{U}_{2}}~{}~{}\alpha||\mathcal{P}_{\Omega}(\mathcal{X}-\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{2}\bm{U}_{2}^{\top}\times_{3}\bm{U}_{3}^{\top})||_{F}^{2}+(1-\alpha)||\bm{y}-\bm{1}_{M}\cdot\tilde{\beta}_{0}-\bm{S}_{(4)}\cdot\tilde{\bm{\beta}}_{1}||_{F}^{2},$

which is equivalent to the following problem when data is incomplete:

$\arg\min_{\bm{U}_{2}}~{}~{}\alpha||\mathcal{X}-(\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{2}\bm{U}_{2}^{\top}\times_{3}\bm{U}_{3}^{\top})\odot\text{logic}(\mathcal{X})||_{F}^{2}+(1-\alpha)||\bm{y}-\bm{1}_{M}\cdot\tilde{\beta}_{0}-\bm{S}_{(4)}\cdot\tilde{\bm{\beta}}_{1}||_{F}^{2},$

where $\odot$ is the inner product, logic($\mathcal{X}$) denotes the logical
value of $\mathcal{X}$. Since the optimization criterion is convex, it can be
solved by setting the derivatives to be zeros, i.e.,
$\frac{d\Psi}{d\bm{U}_{2}}=\mathbf{0}$, where
$\Psi=\alpha||\mathcal{X}-(\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{2}\bm{U}_{2}^{\top}\times_{3}\bm{U}_{3}^{\top})\odot\text{logic}(\mathcal{X})||_{F}^{2}+(1-\alpha)||\bm{y}-\bm{1}_{M}\cdot\tilde{\beta}_{0}-\bm{S}_{(4)}\cdot\tilde{\bm{\beta}}_{1}||_{F}^{2}$.
This implies
$\frac{d}{d\bm{U}_{2}}(||\mathcal{X}-(\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{2}\bm{U}_{2}^{\top}\times_{3}\bm{U}_{3}^{\top})\odot\text{logic}(\mathcal{X})||_{F}^{2})=\bm{0}$.
According to the communication law of tensor mode multiplication, we have
$\frac{d}{d\bm{U}_{2}}(||\mathcal{X}-[(\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{3}\bm{U}_{3}^{\top})\times_{2}\bm{U}_{2}^{\top}]\odot\text{logic}(\mathcal{X})||_{F}^{2})=\bm{0}$.
As a result,
$\frac{d}{d\bm{U}_{2}}(||\mathcal{X}-(\mathcal{S}_{{U}_{2}}\times_{2}\bm{U}_{2}^{\top})\odot\text{logic}(\mathcal{X})||_{F}^{2})=\bm{0}$,
where
$\mathcal{S}_{{U}_{2}}=\mathcal{S}\times_{1}\bm{U}_{1}^{\top}\times_{3}\bm{U}_{3}^{\top}$.
Furthermore, we have
$\frac{d}{d\bm{U}_{2}}(||\bm{X}_{(2)}-(\bm{U}_{2}^{\top}\cdot\bm{S}_{{U}_{2}(2)})\odot\text{logic}(\bm{X}_{(2)})||_{F}^{2})=\bm{0}$
since $||\mathcal{S}||_{F}^{2}=||\bm{S}_{(n)}||_{F}^{2}$ and
$\mathcal{S}\times_{n}\bm{U}=\bm{U}\cdot\bm{S}_{(n)}$.

Figure 17: An illustration of the data missing pattern in Proposition 10
(stripes representing available columns).

It can be easily shown that $\mathcal{X}$’s mode-2 matricization
$\bm{X}_{(2)}$ has missing columns as well (see Figure 17 for an
illustration). Therefore, similar to the proof of Proposition 9, we have
$\frac{d}{d\bm{U}_{2}}(||\bm{X}_{(2)}^{\pi}-\bm{U}_{2}^{\top}\cdot\bm{S}_{{U}_{2}(2)}^{\pi}||_{F}^{2})=\bm{0}$,
where $\pi$ denotes the indices of available columns in $\bm{X}_{(2)}$,
$\bm{S}_{{U}_{2}(2)}^{\pi}$ denotes a matrix constituting the $\pi$ columns of
$\bm{S}_{{U}_{2}(2)}$. As a result, we have
$2(\bm{X}_{(2)}^{\pi}-\bm{U}_{2}^{\top}\cdot\bm{S}_{{U}_{2}(2)}^{\pi})\cdot(-{\bm{S}_{{U}_{2}(2)}^{\pi}}^{\top})=\bm{0}$.
Thus,
$\bm{U}_{2}^{\top}\cdot\bm{S}_{{U}_{2}(2)}^{\pi}\cdot{\bm{S}_{{U}_{2}(2)}^{\pi}}^{\top}=\bm{X}_{(2)}^{\pi}\cdot{\bm{S}_{{U}_{2}(2)}^{\pi}}^{\top}$,
which gives that
$\bm{U}_{2}^{\top}=\bm{X}_{(2)}^{\pi}\cdot{\bm{S}_{{U}_{2}(2)}^{\pi}}^{\top}\cdot(\bm{S}_{{U}_{2}(2)}^{\pi}\cdot{\bm{S}_{{U}_{2}(2)}^{\pi}}^{\top})^{-1}$.
This yields the analytical solution
$\bm{U}_{2}=(\bm{X}_{(2)}^{\pi}\cdot{\bm{S}_{{U}_{2}(2)}^{\pi}}^{\top}\cdot(\bm{S}_{{U}_{2}(2)}^{\pi}\cdot{\bm{S}_{{U}_{2}(2)}^{\pi}}^{\top})^{-1})^{\top}$.

## References

  * Abdi and Williams (2010) Abdi, H. and L. J. Williams (2010).  Principal component analysis.  Wiley interdisciplinary reviews: computational statistics 2(4), 433–459.
  * Aydemir and Paynabar (2019) Aydemir, G. and K. Paynabar (2019).  Image-based prognostics using deep learning approach.  IEEE Transactions on Industrial Informatics 16(9), 5956–5964.
  * Bogdanoff and Kozin (1985) Bogdanoff, J. L. and F. Kozin (1985).  Probabilistic models of cumulative damage(book).  New York, Wiley-Interscience, 1985, 350 p.
  * Dong et al. (2021) Dong, Y., T. Xia, D. Wang, X. Fang, and L. Xi (2021).  Infrared image stream based regressors for contactless machine prognostics.  Mechanical Systems and Signal Processing 154, 107592.
  * Doray (1994) Doray, L. (1994).  Ibnr reserve under a loglinear location-scale regression model.  In Casualty Actuarial Society Forum, Volume 2, pp. 607–652. Citeseer.
  * Fang et al. (2019) Fang, X., K. Paynabar, and N. Gebraeel (2019).  Image-based prognostics using penalized tensor regression.  Technometrics 61(3), 369–384.
  * Filipović and Jukić (2015) Filipović, M. and A. Jukić (2015).  Tucker factorization with missing data with application to low-n-rank tensor completion.  Multidimensional systems and signal processing 26(3), 677–692.
  * Gebraeel et al. (2009) Gebraeel, N., A. Elwany, and J. Pan (2009).  Residual life predictions in the absence of prior degradation knowledge.  IEEE Transactions on Reliability 58(1), 106–117.
  * Gebraeel et al. (2005) Gebraeel, N. Z., M. A. Lawley, R. Li, and J. K. Ryan (2005).  Residual-life distributions from component degradation signals: A bayesian approach.  IiE Transactions 37(6), 543–557.
  * Hong and Meeker (2010) Hong, Y. and W. Q. Meeker (2010).  Field-failure and warranty prediction based on auxiliary use-rate information.  Technometrics 52(2), 148–159.
  * Hong and Meeker (2013) Hong, Y. and W. Q. Meeker (2013).  Field-failure predictions based on failure-time data with dynamic covariate information.  Technometrics 55(2), 135–149.
  * Jiang et al. (2022) Jiang, Y., T. Xia, D. Wang, X. Fang, and L. Xi (2022).  Adversarial regressive domain adaptation framework for infrared thermography-based unsupervised remaining useful life prediction.  IEEE Transactions on Industrial Informatics.
  * Kolda (2006) Kolda, T. G. (2006).  Multilinear operators for higher-order decompositions.  Technical report, Citeseer.
  * Kolda and Bader (2009) Kolda, T. G. and B. W. Bader (2009).  Tensor decompositions and applications.  SIAM review 51(3), 455–500.
  * Liu et al. (2012) Liu, J., P. Musialski, P. Wonka, and J. Ye (2012).  Tensor completion for estimating missing values in visual data.  IEEE transactions on pattern analysis and machine intelligence 35(1), 208–220.
  * Liu et al. (2013) Liu, K., N. Z. Gebraeel, and J. Shi (2013).  A data-level fusion model for developing composite health indices for degradation modeling and prognostic analysis.  IEEE Transactions on Automation Science and Engineering 10(3), 652–664.
  * Lu et al. (2008) Lu, H., K. N. Plataniotis, and A. N. Venetsanopoulos (2008).  Mpca: Multilinear principal component analysis of tensor objects.  IEEE transactions on Neural Networks 19(1), 18–39.
  * Prautzsch et al. (2002) Prautzsch, H., W. Boehm, and M. Paluszny (2002).  Bézier and B-spline techniques, Volume 6.  Springer.
  * Ramsay and Silverman (2005) Ramsay, J. and B. Silverman (2005).  Principal components analysis for functional data.  Functional data analysis, 147–172.
  * Shu et al. (2015) Shu, Y., Q. Feng, and D. W. Coit (2015).  Life distribution analysis based on l évy subordinators for degradation with random jumps.  Naval Research Logistics (NRL) 62(6), 483–492.
  * Tang et al. (2021) Tang, D., M. Gong, J. Yu, L. Guo, and J. Di (2021).  System-level performance prediction for infrared systems based on energy redistribution in infrared images.  IEEE Transactions on Industrial Electronics 69(2), 2000–2011.
  * Xu et al. (2013) Xu, Y., R. Hao, W. Yin, and Z. Su (2013).  Parallel matrix factorization for low-rank tensor completion.  arXiv preprint arXiv:1312.1254.
  * Xu and Yin (2013) Xu, Y. and W. Yin (2013).  A block coordinate descent method for regularized multiconvex optimization with applications to nonnegative tensor factorization and completion.  SIAM Journal on imaging sciences 6(3), 1758–1789.
  * Yang et al. (2021) Yang, Z., P. Baraldi, and E. Zio (2021).  A multi-branch deep neural network model for failure prognostics based on multimodal data.  Journal of Manufacturing Systems 59, 42–50.
  * Ye and Chen (2014) Ye, Z.-S. and N. Chen (2014).  The inverse gaussian process as a degradation model.  Technometrics 56(3), 302–311.
  * Ye et al. (2014) Ye, Z.-S., M. Xie, L.-C. Tang, and N. Chen (2014).  Semiparametric estimation of gamma processes for deteriorating products.  Technometrics 56(4), 504–513.