IJDS-0001-1922.65

Low-Rank Robust Subspace Tensor Clustering

Low-Rank Robust Subspace Tensor Clustering for Metro Passenger Flow Modeling

Nurretin Dorukhan Sergin School of Computing and Augmented Intelligence,
Arizona State University, Tempe, AZ 85281 Jiuyun Hu School of Computing and
Augmented Intelligence, Arizona State University, Tempe, AZ 85281

Ziyue Li Cologne Institute for Information Systems, University of Cologne,
Cologne, Germany EWI gGmbH, University of Cologne, Cologne, Germany,
<EMAIL_ADDRESS>Chen Zhang Department of Industrial Engineering,
Tsinghua University, Beijing, China<EMAIL_ADDRESS>
Fugee Tsung Department of Industrial Engineering and Decision Analytics, Hong
Kong University of Science and Technology, Hong Kong Information Hub, Hong
Kong University of Science and Technology (Guangzhou), China<EMAIL_ADDRESS>
Hao Yan School of Computing and Augmented Intelligence, Arizona State
University, Tempe, AZ 85281<EMAIL_ADDRESS>
Tensor clustering has become an important topic, specifically in
spatiotemporal modeling, due to its ability to cluster spatial modes (e.g.,
stations or road segments) and temporal modes (e.g., time of the day or day of
the week). Our motivating example is from subway passenger flow modeling,
where similarities between stations are commonly found. However, the
challenges lie in the innate high-dimensionality of tensors and also the
potential existence of anomalies. This is because the three tasks, i.e.,
dimension reduction, clustering, and anomaly decomposition, are inter-
correlated to each other, and treating them in a separate manner will render a
suboptimal performance. Thus, in this work, we design a tensor-based subspace
clustering and anomaly decomposition technique for simultaneously outlier-
robust dimension reduction and clustering for high-dimensional tensors. To
achieve this, a novel low-rank robust subspace clustering decomposition model
is proposed by combining Tucker decomposition, sparse anomaly decomposition,
and subspace clustering. An effective algorithm based on Block Coordinate
Descent is proposed to update the parameters. Prudent experiments prove the
effectiveness of the proposed framework via the simulation study, with a gain
of +25% clustering accuracy than benchmark methods in a hard case. The
interrelations of the three tasks are also analyzed via ablation studies,
validating the interrelation assumption. Moreover, a case study in the station
clustering based on real passenger flow data is conducted, with quite valuable
insights discovered.

Subspace Clustering, Tensor Decomposition, Anomaly Detection, Spatiotemporal
Analysis This paper was first submitted on Oct 8th, 2022

##  1 Introduction

Higher-order tensors have been actively used in research since they have an
inclination to successfully preserve the complicated correlation structures of
data. A tensor can be defined mathematically as multi-dimensional arrays
(Kolda and Bader 2009). The order of a tensor is the number of dimensions,
also known as modes. Tensor clustering is a recent generalization of the basic
one-dimensional clustering problem to a high-dimensional version, and it seeks
to partition an order-$K$ input tensor into coherent sub-tensors (e.g., slices
on one mode) while optimizing some cluster-related measuring criteria (Jegelka
et al. 2009).

Tensor clustering has also become an important topic in spatiotemporal
modeling (Sun and Li 2019, Bahadori et al. 2014, Mao et al. 2022), due to its
ability to cluster complicated spatial modes (e.g., stations, road segments,
etc.) and temporal modes (e.g., time of day or day of week). For example, in
the modeling of passenger flow, similarities between spatial elements are
commonly found in spatio-temporal data (Li et al. 2020a). These clusters may
be due to natural geographical locations (e.g., neighboring stations) or
contextual information (e.g., points of interest) in the context of where the
data are collected. For example, in residential areas, stations often have a
large number of in-flow passengers in the morning on weekdays; instead, in the
business areas, in-flow peaks are expected in the afternoon on weekdays.
Exploring metro station clusters can help public transportation management,
like operational efficiency, Strategic Planning, Anomaly Detection, Policy
Making, and Land Use Planning. Overall, the following two challenges present
for data collected from such complex spatiotemporal systems.

  1. 1.

Complex spatio-temporal structure in high-dimensional data. The dimensionality
of the tensor data is often high (especially in the temporal domain), and
traditional statistical methods may suffer from the “curse of dimensionality”.
Furthermore, data points in the high-dimensional tensor often present
complicated and high-order variations. For the temporal dependencies,
passenger flow patterns show a strong periodicity with a period of 7 days,
which denotes the weekly transit patterns (Li et al. 2020b). This can also be
seen in Figure 1 (right). Furthermore, on the same day, there are also
complicated dependencies. For example, there are often two in-flow peaks for
the multi-functional regions in the daily profile during the weekdays, which
represent the morning and afternoon time for the daily commute, whereas a
dominant morning in-flow peak is usually observed in residential regions and
an evening in-flow peak commonly in business regions. For spatial structures,
we create a correlation heatmap of passenger flows in different stations on
day 1, as shown in Figure 1 (left). Figure 1 (left) demonstrates two kinds of
passenger flow similarities. The geographical similarity is based on the
distance between the stations; while the contextual similarity is based on the
function of the stations.

  2. 2.

Existence of sparse anomalies. The tensor may contain sparse outliers due to
measurement errors or rare events in urban passenger flow, such as weather,
special occasions, etc. The sparse outliers may refer to the sudden increase
of demand; this may be related to some unprecedented events happening, such as
concerts and festivals, or this could even be due to maintenance or breakdowns
in other stations whose spillover effect flows into the station in question.
Identifying these events can also help prepare for the sudden increase in
demand. To better understand the sparse anomaly feature, we create a
correlation heatmap of passenger flows at station 1 on different days, as
shown in Figure 1 (right). Figure 1 (right) shows that in most days, the
passenger flow follows the temporal structure of the similarities of the
weekday and weekend, except for 4 days. These days serve as the sparse anomaly
in the data. Finally, to further illustrate the existence of "sparse
anomalies", we have included another plot on the sparse anomalies, as shown in
Figure 1 (bottom). A good example is illustrated here of an unexpected peak in
the afternoon of January 30th, 2017, at the Fo Tan station. The peak is due to
people leaving home from the Lunar New Year horse racing event, which takes
place on the race course accessible from the Fo Tan station. Therefore, in
this research, we define the “anomaly" study in this research as the entry-
level pieces of the tensor data that do not fit the underlying distribution of
the main features of the data.

Figure 1: Correlation heatmaps showing (Left) geographical and contextual
similarity between stations; (Right) temporal structure and sparse anomaly.
(Bottom): Existence of sparse anomalies

In the literature, most of the methods implement a three-step approach for
tensor clustering. First, tensor decomposition can be used to reduce the
tensor dimensionality. Then, sparse outliers can be removed by robust tensor
decomposition, and lastly, clusters can be identified by existing clustering
algorithms such as K-Means or subspace clustering. In the first step of tensor
decomposition, related to dimension reduction techniques, non-negative tensor
decomposition (Lin et al. 2018), Tucker decomposition (Kolda and Bader 2009),
and CP decomposition (Kolda and Bader 2009) are widely used. In the step of
anomaly detection, to handle sparse outliers, Hu and Work (2021) proposed a
robust tensor decomposition method that can handle outliers data while
achieving dimension reduction. Unlike standard robust tensor decomposition
applications, they used the $L_{2,1}$ norm as they assumed fiber-sparse
anomalies. Finally, in the step of clustering, it is quite common to apply
clustering algorithms to the reduced-dimensional space. For example, Yang et
al. (2019) proposed a hybrid approach that combines tensor decomposition with
spectral clustering in a low-dimensional space.

However, none of these aforementioned works considers the interrelated nature
of clustering, outlier removal, and dimension reduction (Aggarwal 2015), with
the interrelation threefold as follows:

  1. 1.

Dimension Reduction and Clustering: Clustering directly in the high dimensions
will suffer from the “curse of dimensionality”, given that data becomes much
more sparse and Euclidean distance is not a good distance measure in high-
dimensional space. Therefore, dimension reduction techniques are often needed
to project data in a low-dimensional space before clustering methods are
applied (Yang et al. 2019).

  2. 2.

Outlier Decomposition and Clustering: On the other hand, if outliers on
stations are not correctly identified or decomposed, the results of clustering
for the overall data from the normal patterns will also be compromised.

  3. 3.

Dimension Reduction and Outlier Decomposition: Furthermore, a direct dimension
reduction on the outliers-corrupted data will also lead to the introduction of
erroneous information into the lower-dimensional products, i.e., decomposed
latent matrices, which will hurt the downstream task performance (i.e.,
clustering) using the low-dimensional space.

In conclusion, we claim that the tensor clustering problem, dimension
reduction, and outlier decomposition problem are inter-correlated tasks and
tied to each other, and isolated treatment of them may fail to provide
satisfactory performance.

In this paper, we propose a Low-Rank Robust Tensor Subspace Decomposition
(LRTSD) technique to achieve dimensionality reduction, spatial clustering, and
sparse anomaly decomposition simultaneously. The proposed method is inspired
by the recent development in the following three fields: Tucker decomposition
(Kolda and Bader 2009), tensor subspace clustering (Vidal 2011, Zhang et al.
2015), and sparse anomaly decomposition (Yan et al. 2017, Li et al. 2022a) (In
this paper, we interchangeably use the term outlier and anomaly to represent
the data which behaves outstandingly). Tucker decomposition is used as a
dimensionality reduction technique, which reduces the original high-
dimensional tensor data into a lower-dimensional core tensor while fitting
orthonormal bases (Tucker 1966). Furthermore, the proposed method is inspired
by the sparse anomaly decomposition methods (Yan et al. 2017, Li et al.
2022a), which aims to decompose the sparse anomaly component from the
background component, with the $L_{1}$ norm applied to encourage the sparsity
of detected anomalies. Finally, the proposed method is also inspired by the
subspace clustering method (Vidal 2011), which has achieved especially good
accuracy in high-dimensional clustering problems (Parsons et al. 2004).

In this paper, we innovatively apply subspace clustering to the core tensor in
the Tucker Decomposition, which is inherently more resilient to
corruption/anomaly. To our knowledge, this is the first instance of such
integration, and our simulation studies substantiate its superior performance.
The amalgamation of these three components significantly increases the
complexity of the model. To address this, we have developed an efficient
optimization algorithm to overcome this challenge, ensuring the
generalizability, computational feasibility, and enhanced accuracy of the
tensor clustering model.

The remainder of this paper is organized as follows. Section 2 gives a brief
overview of the basic tensor notations and multilinear algebra operation and
subspace clustering. Section 3 provides a detailed introduction to the
proposed subspace Tucker decomposition methods as well as an efficient
optimization algorithm. Section 4 conducts a simulation study, and Section 5
apply the proposed method to a real dataset based on Hong Kong Metro data.
Section 6 concludes and presents some future works.

##  2 Literature Review

In this section, we will briefly review the methodology related to tensor
decomposition and subspace clustering methods for tensor data.

###  2.1 Tensor Decomposition

Related to robust tensor decomposition methods, most existing works apply the
dimension reduction approach and tensor decomposition method for feature
extraction.

Historically, all data is vectorized to the sample dimension, and the ordinary
principal component analysis (PCA) is applied to the vectorized data to
extract and monitor features, which is typically termed vectorized PCA (VPCA)
(Nomikos and MacGregor 1994). However, such a method sacrifices the detection
power, since the vectorization on high-dimensional data destroys the innate
cross-dimensional relations.

Tensor is introduced as an efficient data structure to preserve inter-
dimensional correlations well, such as complex spatiotemporal correlations
among traffic data, as mentioned before.  Tensor decomposition is developed
accordingly to extract the features from the tensor data. Various tensor
decomposition methods are developed, namely: (1) CANDECOMP/PARAFAC (CP)
Decomposition (Hitchcock 1927), which represents a tensor
$\mathcal{X}\in\mathbb{R}^{I_{1}\times I_{2}\times\dots\times I_{K}}$ as the
weighted summation of a set of rank-one tensors; (2) Tucker Decomposition
(Tucker 1966), which decomposes a tensor into a core tensor
$\mathcal{C}\in\mathbb{R}^{J_{1}\times J_{2}\times\dots\times J_{K}}$
multiplied by a mode matrix $\bm{U}^{(k)}\in\mathbb{R}^{I_{k}\times J_{k}}$
along each dimension. The proposed method is based on the Tucker
decomposition, and detailed operations of the Tucker decomposition will be
introduced in Section 3.1.

Yan et al. (2015) proposed to utilize several anomaly detection and monitoring
schemes based on Tucker Decomposition (Tucker 1966), Multi-linear PCA (Lu et
al. 2008a), CP Decomposition (Hitchcock 1927), and Uncorrelated multi-linear
PCA (Lu et al. 2008b) with the application for image-based anomaly detection.
However, such approaches do not consider the properties of sparse outliers for
anomaly detection. To address this problem, tensor-based sparse anomaly
decomposition techniques have been proposed (Fanaee-T and Gama 2016) and
applied into video-imaging data (Yan et al. 2018), urban traffic flow
prediction (Sofuoglu and Aviyente 2021), crime rate monitoring (Zhao et al.
2021), pandemic monitoring (Zhao et al. 2020), etc. However, the major
limitation of these existing tensor-based anomaly decomposition methods is
that they assume that the tensor follows a joint smooth and low-rank
representation, which is not feasible for high-dimensional tensor data with
multiple latent clusters behind the data distribution (Li et al. 2020a).
However, a direct application of the decomposition method to the tensor data
with outliers for clustering will render an oversimplified model.

###  2.2 Tensor Subspace Clustering

On the other hand, there are some recent literature focusing on clustering
methods for tensor data. (Drakopoulos et al. 2019) gave an extensive
methodology review of the tensor clustering. Some methods such as topic-model-
based clustering (Li 2021, Li et al. 2022b), dynamic tensor clustering (Sun
and Li 2019), multi-view subspace clustering (Zhang et al. 2015), dynamic
subspace clustering (Zhang et al. 2020), joint K-means and high-order SVD
(Huang et al. 2008) are proposed. Among them, subspace clustering has been a
popular method for high-dimensional clustering (Parsons et al. 2004), which
learns data representation in certain low-dimensional subspaces and clusters
of the data points. The subspace clustering is commonly formulated based on
the data’s self-expression property (Gao et al. 2015), representing the
original input data $\bm{X}$ by itself: $\bm{X}=\bm{X}\bm{Z}+\bm{E}$, where
$\bm{Z}\in\mathbb{R}^{n\times n}$ is the subspace representation matrix, and
the nonzero elements in $\bm{Z}$ correspond to the data points from the same
subspace. Different properties could be achieved by introducing regularization
terms into the representation matrix. More details will be explained in
Section 3.2. However, like the other tensor clustering methods mentioned
above, they are typically operated on the high-dimensional tensor directly and
may suffer from the curse of dimensionality. For example, similar as tensor
regression (Lock 2018, Gahrooei et al. 2021), directly applying self-
expression on high-dimensional data leads to an even higher-dimensional self-
expressiveness term.

Recently, there have been some works focusing on combining tensor low-rank
decomposition and subspace clustering (Fu et al. 2014, 2016). Specifically, Fu
et al. (2016) proposed a subspace clustering that incorporates sparse
dictionary learning into tucker decomposition, where the tucker core tensor
could be the input tensor itself (Fu et al. 2014) or the inverse matricization
product of the dictionary and sparse representation (Fu et al. 2016). However,
these subspace clustering methods are directly targeted at high-dimensional
tensors, which are not suitable for tensors with sparse anomalies. Recently,
there have been some efforts to combine subspace clustering methods with
anomaly detection to handle the challenges of anomalies in multi-clustered
datasets.

Recently, there have been a few developed tensor clustering methods based on
deep learning and neural networks. (Wang et al. 2023) proposed the Tensorized
Hypergraph Neural Networks (THNN) to learn the hypergraph structure on higher
dimensions. In two-dimensional cases, further methods, such as spectral
clustering, can be performed to cluster the nodes. However, this method only
considered the hypergraph structure and did not take subspace learning and
anomaly into consideration. (Zhao et al. 2022a) proposed the Reinforced Tensor
Graph Neural Network (RTGNN) to learn multi-view graph data. This method
contains three modules, which makes it time-consuming. Tucker decomposition
only takes part in the initialization, and the anomaly is neglected in this
method. (Huang et al. 2023) extensively explored various methods, specifically
on hyperspectral image clustering. The improvement achieved by neural network-
based methods compared to tensor decomposition and self-regression-based
methods is data-dependent. However, compared to the tensor-based approach,
neural networks tend to become time-consuming and require many training
samples.

###  2.3 Sparse Anomaly Decomposition

Anomaly is usually detected based on the assumption that anomalous data do not
conform to expected behavior with outstandingly different features from the
homogeneous background (Chandola et al. 2009). Traditional methods detect the
anomaly based on distance and density, which assume that anomalies lie far
away from the background or locate in a less dense area. Various definitions
of distance or density have been introduced (Du and Zhang 2014). For example,
(Sun et al. 2022) enhanced the CP decomposition by a deep CP decomposition
neural network for feature extraction with an anomaly detection step.
Nonetheless, to further achieve the separation of the anomalies, a further
step such as clustering is still necessary. In this case, a one-step
decomposition is favored in practice.

This paper is also motivated by the recent development of the decomposition-
based methods (Sofuoglu and Aviyente 2021, Yan et al. 2017, Li et al. 2022a),
which try to decompose the original input data into several components
directly. For example, Yan et al. (2017) decomposed input data into the smooth
background and sparse anomaly, with smoothness regularization on the
background and sparsity penalty on the anomaly. Later, the extensions to
spatiotemporal data (Yan et al. 2018), deep learning-based decomposition
models (Zhao et al. 2022b) are proposed. For example, (Li et al. 2022a)
considered more components, namely, smooth background, sparse anomaly, sample-
specific deviation, and random noise, and they introduced a transfer learning
mechanism (Pan and Yang 2009) to solve the cold-start problem for anomaly
decomposition. However, these methods are not developed for tensor data, and
the decomposed components are pre-defined in a fixed manner. Recently, Shen et
al. (2022) proposed a tensor-based decomposition method to decompose the
static tensor background and smooth foreground. However, the method is
designed for static background and cannot be directly used to learn the tensor
data from different subspaces. Our method instead is to combine the sparse
decomposition method with the tensor decomposition approach, and it learns the
hidden components in a data-driven and unsupervised manner via tensor
subspace.

To summarize, this work is the first that achi†eves dimensionality reduction,
spatial clustering, and anomaly decomposition simultaneously.

##  3 Methodology

In this section, we will first review the basic tensor algebra in Section 3.1.
Then, we will give a brief review of the recent development of subspace
clustering in Section 3.2. More specifically, we will discuss the proposed
formulation combining tensor subspace clustering in Section 3.3. We will then
discuss how to solve the formulation efficiently in Section 3.4.

###  3.1 Basic Tensor Notation and Multilinear Algebra

In this section, we introduce basic notations, definitions, and operators in
multilinear (tensor) algebra that is used in this work. Throughout the
chapter, scalars are denoted by lowercase italic letters, e.g., $a$, vectors
by lowercase boldface letters, e.g., $\bm{a}$, matrices by uppercase boldface
letter, e.g., $\bm{A}$, and tensors by calligraphic letters, e.g.,
$\mathcal{A}$. For example, an order-$K$ tensor is represented by
$\mathcal{A}\in\mathbb{R}^{I_{1}\times\cdots\times I_{K}}$, where $I_{k}$
represents the mode-$k$ dimension of $\mathcal{A}$. The mode-$k$ product of a
tensor $\mathcal{A}$ by a matrix $\bm{V}\in\mathbb{R}^{P_{k}\times I_{k}}$ is
defined as follows:

$(\mathcal{A}\times_{k}\bm{V})(i_{1},\cdots,i_{k-1},j_{k},i_{k+1},\cdots,i_{K})=\sum_{i_{k}}A(i_{1},\cdots,i_{k},\cdots,i_{K})V(j_{k},i_{k}).$

The Frobenius norm of a tensor $\mathcal{A}$ can be defined as follows:

$\|\mathcal{A}\|_{F}^{2}=\sum_{i_{1},\cdots,i_{K}}A(i_{1},\cdots,i_{k},\cdots,i_{K})^{2}.$

The $n$-mode unfold operator maps the tensor $\mathcal{A}$ into matrix
$\bm{A}_{(n)}$, where the columns of $\bm{A}_{(n)}$ are the $n$-mode vectors
of $\mathcal{A}$.

Tucker decomposition decomposes a tensor into a core tensor multiplied by a
matrix along each mode, formulated as follows:

$\mathcal{A}=\mathcal{C}\times_{1}\bm{U}^{(1)}\times_{2}\bm{U}^{(2)}\cdots\times_{K}\bm{U}^{(K)},$

where $\mathcal{C}\in\mathbb{R}^{J_{1}\times J_{2}\dots\times J_{K}}$ is the
core tensor, $\bm{U}^{(k)}$ is an orthogonal $I_{k}\times J_{k}$ matrix and is
a principal component mode-$k$. The definition of Kronecker product is as
follows: Suppose $\bm{A}\in\mathbb{R}^{m\times n}$ and
$\bm{B}\in\mathbb{R}^{p\times q}$ are matrices, the Kronecker product of these
matrices, denoted by $\bm{A}\otimes\bm{B}$, is an $mp\times nq$ block matrix,
which can be formalized as follows:

$\bm{A}\otimes\bm{B}=\left[\begin{array}[]{ccc}a_{11}\bm{B}&\cdots&a_{1n}\bm{B}\\\
\vdots&\ddots&\vdots\\\ a_{m1}\bm{B}&\cdots&a_{mn}\bm{B}\end{array}\right]$

###  3.2 Subspace Clustering

An important but often unnoticed assumption of PCA is that the data lives in a
_single_ low-dimensional subspace. This assumption may be too crude for many
applications, including spatiotemporal problems such as motion segmentation
(Vidal et al. 2008). It is more likely in such cases that the data lives in a
mixture of subspaces instead of a single one (see Figure 2). The field of
subspace clustering deals with problems of this nature and has attracted
growing interest in recent years (Liu et al. 2019, Guo et al. 2019).

Figure 2: Data were drawn from two one-dimensional subspaces (red and blue)
and one two-dimensional subspace (green). Adapted from (Vidal 2011).

The central problem of subspace clustering is the grouping of points and their
respective clusters are both unknown in the beginning. Various frameworks have
been proposed in recent years, ranging from direct algebraic solutions to
probabilistic approaches (Vidal 2011), whereas spectral clustering (von
Luxburg 2007) dominates the approaches proposed in recent years due to their
superior performance over other frameworks. In short, spectral clustering aims
to build a similarity graph where nodes are individual samples and the edges
are similarity links. Algorithms such as graph cut can be used later to
identify the clusters.

An effective way to obtain similarity is to use the self-expressiveness
property (Parsons et al. 2004). Let $\bm{D}\in\mathbb{R}^{n\times p}$ be a
data matrix. The self-expressiveness can be defined with the equality
$\bm{D}=\bm{Z}\bm{D}$ where the self-expression matrix
$\bm{Z}\in\mathbb{R}^{n\times n}$ helps to express each point as a linear
combination of other points in the data. This matrix can later be used to find
clusters, as explained above.

Spectral subspace clustering methods differ in how they regulate the affinity
matrix to avoid the trivial solution of identity matrix $\bm{Z}=\bm{I}$.
Sparse subspace clustering (SSC) imposes sparsity on the affinity matrix (Yang
et al. 2015), which is formulated as follows:

$\text{$\underset{\bm{Z}}{\text{\text{min}}}$}\quad\,\|\,\bm{Z}\,\|\,_{1}+\frac{\lambda}{2}\,\|\,\bm{D}-\bm{Z}\bm{D}\,\|\,_{F}^{2}\quad
s.t.\quad\text{\text{diag}}(\bm{Z})=0,$

where the $\|.\|_{1}$ represent the L1 penalty of the self-expression matrix.
By applying this penalty, the self-expression matrix $\bm{Z}$ is also assumed
to be low-rank. Low-rank representation method was proposed by Liu et al.
(2013). Similarly, another method of low-rank subspace clustering was
developed by Vidal and Favaro (2014) where the core formulation uses the
nuclear-norm of the self-expression matrix $\bm{Z}$ defined by the sum of its
singular values as the best convex approximation of the rank of $\bm{Z}$.

$\underset{\bm{Z},\bm{D},\bm{A}}{\text{min}}\quad\,\|\,\bm{Z}\,\|\,_{*}+\frac{\tau}{2}\,\|\,\bm{A}-\bm{A}\bm{Z}\,\|\,_{F}^{2}+\frac{\alpha}{2}\,\|\,\bm{D}-\bm{A}\,\|\,_{F}^{2}\quad
s.t.\quad\bm{Z}^{\top}=\bm{Z}$

Thresholding ridge regression (TRR) (Peng et al. 2015) is an improvement on
SSC in the sense that it relaxes the requirement to assume the structure of
the affinity matrix, whether it is low-rank or sparse. The formulation is also
much simpler.

$\underset{\bm{Z}}{\text{min}}\quad\lambda\,\|\,\bm{Z}\,\|\,^{2}_{F}+\,\|\,\bm{D}-\bm{D}\bm{Z}\,\|\,_{F}^{2}$

The regularization parameter $\lambda$ represents the trade-off between the
reconstruction fidelity and clustering ability, which may be set with respect
to the use case. However, all these variations focus on the matrix formulation
and cannot be applied to the tensor formulation.

###  3.3 Proposed Low-Rank Robust Tensor Subspace Decomposition (LRTSD) Model

Without loss of generality, we will illustrate the proposed method using a
three-dimensional tensor. The methodology can be easily extended into higher-
order tensors. For a three-dimensional tensor, we assume that
$\mathcal{X}\in\mathbb{R}^{N\times I_{2}\times I_{3}}$ to represent the
original data, where the first dimension denotes the sample dimension that
requires the clustering. Furthermore, we assume that this tensor $\mathcal{X}$
can be decomposed into the normal variations data $\mathcal{L}$ and the
outlier tensor $\mathcal{A}$ (i.e., represents anomalous inflows in the
passenger flow data), which is assumed to be sparse. We assumed the normal
variations data $\mathcal{L}$ has a low-rank structure, which is a common
assumption for high-dimensional tensors for dimension reductions. Furthermore,
without loss of generality, we assume the first mode of tensor $\mathcal{X}$
is the dimension that needs further clustering. We name the first mode of the
tensor (mode 1) as the clustering mode and the other modes (modes 2 and 3) as
the non-clustering modes.

For example, the passenger flow data can be represented as a tensor
$\mathcal{X}\in\mathbb{R}^{N\times I_{2}\times I_{3}}$, with the three modes
representing the station (with $N$ stations), the day of the week (with
$I_{2}=7$ days), and the 5-minute time intervals (with $I_{3}=247$ intervals
from 5 am to 1:30 am) within each day for the period where the subway station
is open for, respectively. Each element of the tensor $X_{i,d,t}$ represents
the inflow into stations $i$ at day $d$ and time interval $t$. Here, the
station $i$ is the clustering mode, and $d$ and $t$ are the non-clustering
modes, given that the station $s$ is often heterogeneous and may exhibit
multiple clusters.

In this case, we propose to use a partial Tucker decomposition on the non-
clustering modes as:

$\mathcal{L}=\mathcal{C}\times_{2}\bm{U}_{2}\times_{3}\bm{U}_{3},$ (1)

where a core tensor $\mathcal{C}\in\mathbb{R}^{N\times P_{2}\times P_{3}}$
reduces the dimensionality of the original tensor. $\bm{U}_{2}$ and
$\bm{U}_{3}$ are respective orthogonal matrices for the non-clustering modes.
We propose to find subspaces in the latent space on the core tensor
$\mathcal{C}$ as opposed to original space, as this is computationally more
efficient. Here, the notations that will be used in the paper are given in
Table 1.

Table 1: Notations for Proposed LRTSD Model. Notation |  Explanation
---|---
$\mathcal{X}\in\mathbb{R}^{N\times I_{2}\times I_{3}}$ |  $N$: number of samples, $I_{2},I_{3}:$ feature dimensions
$\mathcal{C}\in\mathbb{R}^{N\times P_{2}\times P_{3}}$ |  Core tensor, with reduced dimensions $P_{2},P_{3}$
$\bm{U}_{l}\in\mathbb{R}^{I_{l},P_{l}},l=2,3$ |  orthogonal matrices for dimension reduction
$\mathcal{A}\in\mathbb{R}^{N\times I_{2}\times I_{3}}$ |  Anomaly tensor, same size as $\mathcal{X}$
$\bm{Z}\in\mathbb{R}^{N\times N}$ |  Self-expression matrix

Using the thresholding ridge regression framework (Peng et al. 2015), we
utilize the self-expressive expression
$\mathcal{C}-\mathcal{C}\times_{1}\bm{Z}$ by regularizing the self-expression
matrix $\bm{Z}\in\mathbb{R}^{N\times N}$ using Frobenius norm. Formally, the
model is described as an optimization problem:

$\displaystyle\underset{\bm{Z},\mathcal{C},\bm{U}_{2},\bm{U}_{3},\mathcal{A}}{\text{min}}$
$\displaystyle\begin{aligned}
\frac{1}{2}\,\|\,\bm{Z}\,\|\,_{F}^{2}+\frac{\lambda_{z}}{2}\,\|\,\mathcal{C}-\mathcal{C}\times_{1}\bm{Z}\,\|\,_{F}^{2}+\lambda_{a}\,\|\,\mathcal{A}\,\|\,_{1}\\\
+\frac{\lambda_{e}}{2}\,\|\,\mathcal{X}-\mathcal{A}-\mathcal{C}\times_{2}\bm{U}_{2}\times_{3}\bm{U}_{3}\,\|\,_{F}^{2}\end{aligned}$
(2) $\displaystyle\bm{U}_{2}^{\top}\bm{U}_{2}=\bm{I}$ (3)
$\displaystyle\bm{U}_{3}^{\top}\bm{U}_{3}=\bm{I}$ (4)

In summary, each term in the Equation (2) achieves following purpose:

  * •

$\,\|\,\mathcal{X}-\mathcal{A}-\mathcal{C}\times_{2}\bm{U}_{2}\times_{3}\bm{U}_{3}\,\|\,_{F}^{2}$
is the error term between the original input data $\mathcal{X}$ and the two
decomposed components, i.e., anomaly $\mathcal{A}$ and low rank structure
$\mathcal{L}=\mathcal{C}\times_{2}\bm{U}_{2}\times_{3}\bm{U}_{3}$;

  * •

$\,\|\,\mathcal{C}-\mathcal{C}\times_{1}\bm{Z}\,\|\,_{F}^{2}$ is the self-
expression term in tensor subspace clustering and is imposed on the decomposed
core tensor. $\,\|\,\bm{Z}\,\|\,_{F}^{2}$ is the ridge regularization;

  * •

$\,\|\,\mathcal{A}\,\|\,_{1}$ is the sparsity regularization term for the
anomaly.

To this end, this model achieves what we promised “a simultaneous framework
for dimension reduction, spatial clustering, and anomaly detection”.

###  3.4 Optimization Procedure

In this section, we propose an efficient algorithm to solve the proposed
optimization problem Equation 2 \- Equation 4. We will utilize the Block
Coordinate Descent procedure to obtain each one of the model parameters blocks
$\bm{Z},\mathcal{C},\bm{U}_{2},\bm{U}_{3}$ and $\mathcal{A}$ iteratively while
keeping the others fixed from the last iteration. When Algorithm 1 terminates
through convergence, graph cut is applied to $\bm{L}=\bm{Z}^{\top}+\bm{Z}$ to
achieve clusters, and $\mathcal{A}$ can be analyzed for detected point
anomalies. The rationale of the proofs of the propositions are relegated to
the Appendix.

Data: $\mathcal{X}$, $N$, $P_{2}$, $P_{3}$

Result: $\mathcal{A}$, $\mathcal{C}$, $\bm{U}_{2}$, $\bm{U}_{3}$, $\bm{Z}$

Initialize $\mathcal{A}=\mathbf{0}$, $\bm{Z}=\bm{I}$. Initialize
$\mathcal{C}$, $\bm{U}_{2}$, $\bm{U}_{3}$ by applying Tucker decomposition on
$\mathcal{X}$;

while  _not converged_ do

Update $\bm{Z}$ as in Equation 6 ;

Update $\mathcal{A}$ as in Equation 8 ;

Update $\mathcal{C}$ as in Equation 10 ;

Update $\bm{U}_{2}$ and $\bm{U}_{3}$ as in Equation 12 ;

end while

Algorithm 1 Low-rank and Sparse Tensor Decomposition with Ridge Regularized
Subspace Clustering

####  3.4.1 Calculating the Self-expression Matrix

The subproblem for optimizing $\bm{Z}$ when all other parameters are fixed is

$\displaystyle\underset{\bm{Z}}{\text{min}}$
$\displaystyle\frac{1}{2}\,\|\,\bm{Z}\,\|\,_{F}^{2}+\frac{\lambda_{z}}{2}\,\|\,\bm{C}_{(1)}-\bm{Z}\bm{C}_{(1)}\,\|\,_{F}^{2}$
(5)

###### Proposition 3.1

The optimization problem in Equation 5 admits a closed-form solution,

$\displaystyle\bm{Z}^{*}=\lambda_{z}(\bm{I}+\lambda_{z}\bm{C}_{(1)}\bm{C}_{(1)}^{\top})^{-1}(\bm{C}_{(1)}\bm{C}_{(1)}^{\top})$
(6)

Equation (5) suggests that there are only two terms in the problem related to
the self expression matrix $\bm{Z}$. Proposition 3.1 derived a closed-form
solution.

####  3.4.2 Calculating the Sparse Anomaly Tensor

The subproblem for optimizing $\mathcal{A}$ when all other parameters are
fixed is

$\displaystyle\underset{\mathcal{A}}{\text{min}}\lambda_{a}\,\|\,\mathcal{A}\,\|\,_{1}+\frac{\lambda_{e}}{2}\,\|\,\mathcal{X}-\mathcal{A}-\mathcal{C}\times_{2}\bm{U}_{2}\times_{3}\bm{U}_{3}\,\|\,_{F}^{2}$
(7)

###### Proposition 3.2

The optimization problem in Equation 7 admits a closed-form solution using the
soft-thresholding operator,

$\mathcal{A}^{*}=\mathrm{sgn}(\mathcal{X}-\mathcal{C}\times_{2}\bm{U}_{2}\times_{3}\bm{U}_{3})\odot\max(\bm{0},|\mathcal{X}-\mathcal{C}\times_{2}\bm{U}_{2}\times_{3}\bm{U}_{3}|-\frac{\lambda_{a}}{\lambda_{e}})$
(8)

The proof follows directly on using soft thresholding in solving the
orthogonal LASSO problem and is therefore omitted here.

####  3.4.3 Calculating the Core Tensor

The subproblem for optimizing $\mathcal{C}$ is

$\displaystyle\underset{\mathcal{C}}{\text{min}}\frac{\lambda_{z}}{2}\,\|\,\mathcal{C}-\mathcal{C}\times_{1}\bm{Z}\,\|\,_{F}^{2}+\frac{\lambda_{e}}{2}\,\|\,\mathcal{X}-\mathcal{A}-\mathcal{C}\times_{2}\bm{U}_{2}\times_{3}\bm{U}_{3}\,\|\,_{F}^{2}$
(9)

###### Proposition 3.3

Define $\bm{M}\triangleq\bm{X}_{(1)}-\bm{A}_{(1)}$ and
$\bm{P}\triangleq(\bm{U}_{3}\otimes\bm{U}_{2})^{\top}$. The optimization
problem in Equation 9 admits a closed-form solution

$\displaystyle\bm{C}_{(1)}^{*}=\lambda_{e}(\lambda_{z}(\bm{I}-\bm{Z})^{\top}(\bm{I}-\bm{Z})+\lambda_{e}\bm{I})^{-1}\bm{M}\bm{P}^{\top}$
(10)

The resulting matrix $\bm{C}_{(1)}^{*}$ is folded back to the original tensor
$\mathcal{C}$ after the calculation.

####  3.4.4 Calculating the Orthonormal Bases

When all other parameters are fixed, the subproblem to update $\bm{U}_{2}$ is

$\displaystyle\underset{\bm{U}_{2}}{\text{min}}$
$\displaystyle\frac{\lambda_{e}}{2}\,\|\,\mathcal{X}-\mathcal{A}-\mathcal{C}\times_{2}\bm{U}_{2}\times_{3}\bm{U}_{3}\,\|\,_{F}^{2}$
(11) $\displaystyle\bm{U}_{2}^{\top}\bm{U}_{2}=\bm{I}$

###### Proposition 3.4

Define
$\bm{Q}\triangleq(\bm{X}_{(2)}-\bm{A}_{(2)})(\bm{U}_{3}\otimes\bm{I})\bm{C}_{(2)}^{\top}$
and let $\bm{Q}=\hat{\bm{U}}\bm{D}\hat{\bm{V}}^{\top}$ be the singular value
decomposition of $\bm{Q}$. Thus, the optimization problem in Equation 11
admits a closed-form solution for the optimal values for the first orthogonal
basis $\bm{U}_{2}^{*}$

$\displaystyle\bm{U}_{2}^{*}=\hat{\bm{U}}\hat{\bm{V}}^{\top}$ (12)

Finally, the update of $\bm{U}_{3}^{*}$ is similar, therefore, omitted. In
conclusion, the update of $\bm{U}_{2}^{*}$ and $\bm{U}_{3}^{*}$ both yield
closed-form solutions.

##  4 Simulation Study & Results

In this section, we will evaluate the proposed methods via a simulation study.
We will first introduce the simulation setup in Section 4.1. Then, the
proposed method will be evaluated and compared with several clustering-based
benchmark methods in Section 4.3.

###  4.1 Simulation Setup

In this section, we formalize the simulation procedure we use to generate
clustered 3-order tensors with dimensionality $N\times I_{2}\times I_{3}$.
First, the procedure calls for the parameters in Table 2. After defining these
parameters and setting the random number generating seed, we follow the
following procedure to generate the data:

  * •

For each cluster $k$ and mode of tensor $l=2,3$, generate random factor
matrices $\bm{U}_{k,l}\in\mathbb{R}^{I_{l}\times P_{l}}$ from standard normal
distribution, then orthogonalize $\bm{U}_{k,l}$ via the QR decomposition.

  * •

For each cluster $k$, generate a random core tensor
$\mathcal{C}_{k}\in\mathbb{R}^{N\times P_{2}\times P_{3}}$ from standard
normal distributions.

  * •

For each cluster $k$, generate
$\mathcal{X}_{k}=\mathcal{C}_{k}\times_{2}\bm{U}_{k1}\times_{3}\bm{U}_{k2}+\epsilon$,
where noise is added on each element of $\mathcal{X}$ as
$\epsilon\sim\mathcal{N}(0,\sigma^{2})$,  with $\sigma^{2}=0.25$.

  * •

Elementwise add sparse anomaly
$\mathcal{X}\leftarrow\mathcal{X}+\psi\mathcal{A}_{p}$. Each entry of
$\mathcal{A}_{p}$ is generated from the Bernoulli distribution with
probability. Here $p$ is the anomaly ratio. In each entry where
$\mathcal{A}_{p}$ is $1$, indicating anomaly is present, we add the anomaly
intensity constant $\psi$ to the data.

Figure 3: Illustration depicting the generation process of simulated data Table 2: Required parameters for simulated clustered data generation and their notation. Notation |  Explanation
---|---
$K\in\mathbb{Z}^{+}$ |  Number of clusters
$N\in\mathbb{Z}^{+}$ |  Number of samples per cluster
$D_{l}\in\mathbb{Z}^{+}$ |  Ambient dimension of $l$th mode
$P_{l}\in\mathbb{Z}^{+}$ |  Intrinsic dimension of $l$th mode
$\sigma\in\mathbb{R}^{+}$ |  Standard deviation of random noise
$\mathcal{A}_{p}\in\mathbb{B}^{D_{1}\times\dots\times D_{L}}$ |  Sparse anomaly index tensor with $0<p<1$ of entries being $1$ and others $0$
$\psi\in\mathbb{R}$ |  Sparse anomaly intensity

Figure 3 provides a visual representation of the simulation data generation
process. This procedure allows us to control key elements of the data
generating process, such as the number of clusters, the cluster dimensions,
the prevalence of anomalies, and how gross the anomalies are. We will now
experiment with these factors and their effect on how well our model clusters
samples and/or detects anomalies.

###  4.2 Tuning Parameters Selection

The proposed LRTSD has three tuning parameters: $\lambda_{e}$ penalizes the
norm of the self-regression matrix; $\lambda_{z}$ penalizes the self-
regression error of the core tensor; $\lambda_{a}$ penalizes the norm of the
anomaly tensor. These three tuning parameters are in comparison to the
reconstruction error term, which has a coefficient of 1.

Related to the anomaly tuning parameter $\lambda_{a}$, one can select based on
the percentage of nonzero anomaly entries expected in the dataset. For
example, we can select the anomaly entries as a fixed percentage if we have
some rough ideas on the percentage of the anomaly entries. If such a
percentage is unknown, Otsu’s method can be used to automatically select the
tuning parameter $\lambda_{a}$, which aims to search for $\lambda_{a}$ that
minimizes the intra-class variance, defined as a weighted sum of variances of
the two classes: normal class and abnormal class. This procedure is discussed
in detail in (Yan et al. 2017). To select other tuning parameters, such as
$\lambda_{e}$ and $\lambda_{z}$, we propose using the normalized cut score.
For $\lambda_{e}$ and $\lambda_{z}$, $\lambda_{e}$ yields smaller diagonal
values in the self-regression matrix $\bm{Z}$; higher $\lambda_{z}$ yields a
more precise self-regression. To address the challenge of unknown clustering
accuracy, we propose using the normalized cut to evaluate the affinity matrix
produced by the algorithm (Zhang et al. 2023) and select the tuning parameters
$\lambda_{e}$ and $\lambda_{z}$ that produce the most reasonable affinity
matrix by minimizing the normalized cut score. This measurement is defined by
$NC=\sum_{c=1}^{C}\frac{W_{c}^{out}}{W_{c}^{in}+W_{c}^{out}}$, where $c$ is
the index of the cluster and $W_{c}^{in}$ and $W_{c}^{out}$ are the sum of the
weights of the graph within the cluster and the weights of the graph outside
the cluster. Smaller NC means better clustering.

Figure 4:  Clustering Accuracy of Different Tuning Parameter Combinations

Some examples of clustering accuracy for different combinations of tuning
parameters are shown in Figure 4. From Figure 4, we can see that (1)
$\lambda_{e}$, $\lambda_{z}$, $\lambda_{a}$ cannot be too small, since we will
almost lose a penalty term when any of the parameters are too small, which
affects the clustering accuracy. (2) As $\lambda_{e}$ and $\lambda_{z}$
increase, the clustering precision increases. When $\lambda_{e}$ continues to
increase, the clustering accuracy will decrease because the corresponding
terms will dominate optimization. (3) $\lambda_{a}$ cannot be too small or too
big, because we will not penalize the anomaly term when it is too small, and
it will be too hard to discover any anomaly when $\lambda_{a}$ is too big.

###  4.3 Performance Evaluation and Comparison

We generate a simulated three-dimensional tensor with three groups, thirty
samples per cluster, fifty as the ambient dimensionality of each mode, and
five as the intrinsic dimension on each mode. Therefore, the dimensionality of
the data tensor is $90\times 50\times 50$. This can be considered
representative of a spatio-temporal observation with one spatial dimension and
two temporal dimensions with the inherent clustered structure and redundant
information with respect to dimensionality. In this section, we compare the
clustering accuracy of our models and three benchmark models under different
anomaly ratios, and anomaly intensity, where anomaly ratios are defined as the
ratio of anomaly entry compared to the entire tensor dimensions and anomaly
intensity denotes $\psi$. The clustering accuracy is defined by the maximum
proportion of matching labels of the $90$ objects in the first dimension over
all permutations of the labels. The benchmark methods we compared are
introduced below:

  * •

K-means: For each of the $90$ objects, we reshape the two-dimensional
$50\times 50$ matrix to a vector of length $2500$. We make $3$ groups using
K-means based on the Euclidean distance.

  * •

Robust Tucker+K means: We first use the Robust Tucker decomposition introduced
in (Heng et al. 2023, Xue et al. 2017) to reduce the dimension to $15\times
15$ for each object. We then use K means to generate $3$ clusters based on the
vectorized core tensor.

  * •

Multi-mode Tensor Space Clustering (He and Atia 2022): This is a method based
on low-rank tensor representation (LTRR). We will use single-mode tensor space
clustering for two reasons. Firstly, all the other methods use a single mode,
and thus it makes for a fair comparison. Secondly, when we increased the
number of modes, the running time increased dramatically.

  * •

Robust Subspace Clustering (Peng et al. 2015): This method is based on the
thresholding ridge regression (TRR) and is compatible with the matrix data.
Therefore, we first reshaped the $90\times 50\times 50$ tensor into the
$90\times 2500$ matrix and made $3$ clusters.

When comparing the performance of LRTSD and the benchmark methods, we use a
grid search of the tuning parameters and report the best clustering precision
in different settings, namely $\lambda_{e}$, $\lambda_{c}$, $\lambda_{a}$ in
LRTSD, representation rank $r$ in LTRR, balance and thresholding parameter
$\lambda$ and $k$ in TRR. We report the mean accuracy of 20 replications of
data generation and clustering accuracy for each anomaly ratio & anomaly
intensity combination and each method. The clustering accuracy result is shown
in Table 3 with their standard deviation over the 20 replications in
parentheses. In each combination of anomaly intensity and anomaly ratio, the
best clustering accuracy is highlighted in bold. From the clustering accuracy
result, we can see that the K-means and Robust Tucker + K means methods yield
the two lowest accuracy proportion. This is because they are not designed
spedifically for clustering, and suffer from the curse of dimensionality when
vectorized. When the anomaly ratio and intensity are not large, LRTSD, TRR,
and LTRR have almost perfect clustering accuracy. As the intensity and ratio
of the anomalies increase, the LTRR and TRR start to show errors. The
threshold at which the two methods start to have incorrect clustering labels
is smaller than that of LRTSD. Meanwhile, although LRTSD achieves the highest
clustering accuracy, it does not require much more time compared to benchmark
methods. The average processing times for K-means, Robust Tucker + K means,
LTRR, TRR, and LRTSD are 0.06 second, 0.68 second, 0.08 second, 0.51 second,
0.71 second, respectively.

Table 3: Clustering accuracy for LRTSD and other three benchmarks Ratio | 0.05 | 0.1 | 0.125 | 0.15 | 0.175 | 0.2 | 0.225
---|---|---|---|---|---|---|---
Intensity | Method
4 | K mean | 0.49 | 0.45 | 0.45 | 0.47 | 0.47 | 0.45 | 0.45
(0.02) | (0.03) | (0.03) | (0.02) | (0.02) | (0.02) | (0.03)
Robust Tucker+K means | 0.5 | 0.5 | 0.49 | 0.45 | 0.48 | 0.47 | 0.45
(0.03) | (0.03) | (0.03) | (0.02) | (0.02) | (0.02) | (0.02)
LTRR |  1(0) |  1(0) |  1(0) |  1(0) |  1(0) |  1(0) |  1(0)
TRR | 1 | 1 | 1 | 0.99 | 0.98 | 0.97 | 0.97
(0) | (0) | (0) | (0.01) | (0.01) | (0.01) | (0.01)
LRTSD |  1(0) |  1(0) |  1(0) |  1(0) |  1(0) |  1(0) |  1(0)
5 | K mean | 0.47 | 0.44 | 0.46 | 0.44 | 0.46 | 0.44 | 0.44
(0.02) | (0.02) | (0.03) | (0.03) | (0.03) | (0.03) | (0.03)
Robust Tucker+K means | 0.47 | 0.47 | 0.48 | 0.45 | 0.43 | 0.43 | 0.43
(0.03) | (0.03) | (0.02) | (0.03) | (0.03) | (0.03) | (0.02)
LTRR | 1 | 1 | 1 | 1 | 0.99 | 0.98 | 0.97
(0) | (0) | (0) | (0) | (0.01) | (0.01) | (0.02)
TRR | 0.99 | 0.99 | 0.98 | 0.96 | 0.99 | 0.96 | 0.98
(0.01) | (0.01) | (0.01) | (0.02) | (0.01) | (0.02) | (0.01)
LRTSD |  1(0) |  1(0) |  1(0) |  1(0) |  1(0) |  1(0) |  1(0)
6 | K mean | 0.44 | 0.44 | 0.43 | 0.44 | 0.43 | 0.42 | 0.42
(0.02) | (0.02) | (0.03) | (0.03) | (0.03) | (0.04) | (0.03)
Robust Tucker+K means | 0.46 | 0.46 | 0.45 | 0.43 | 0.41 | 0.41 | 0.4
(0.03) | (0.02) | (0.04) | (0.02) | (0.02) | (0.03) | (0.02)
LTRR | 1 | 1 | 0.99 | 0.99 | 0.96 | 0.94 | 0.87
(0) | (0) | (0.01) | (0.01) | (0.02) | (0.02) | (0.05)
TRR | 0.99 | 0.96 | 0.97 | 0.96 | 0.9 | 0.87 | 0.81
(0.01) | (0.02) | (0.02) | (0.02) | (0.04) | (0.05) | (0.05)
LRTSD |  1(0) |  1(0) |  1(0) |  1(0) |  1(0) |  1(0) |  1(0)
7 | K mean | 0.44 | 0.41 | 0.4 | 0.41 | 0.4 | 0.39 | 0.41
(0.02) | (0.03) | (0.03) | (0.03) | (0.04) | (0.03) | (0.03)
Robust Tucker+K means | 0.45 | 0.43 | 0.43 | 0.4 | 0.39 | 0.4 | 0.39
(0.04) | (0.03) | (0.03) | (0.04) | (0.02) | (0.03) | (0.02)
LTRR | 1 | 0.99 | 0.98 | 0.96 | 0.83 | 0.74 | 0.64
(0) | (0.01) | (0.01) | (0.02) | (0.07) | (0.09) | (0.12)
TRR | 0.97 | 0.94 | 0.9 | 0.79 | 0.76 | 0.64 | 0.63
(0.01) | (0.02) | (0.04) | (0.1) | (0.09) | (0.11) | (0.08)
LRTSD |  1(0) |  1(0) |  1(0) |  1(0) |  1(0) |  1(0) |  1(0)
8 | K mean | 0.44 | 0.39 | 0.39 | 0.4 | 0.39 | 0.41 | 0.4
(0.02) | (0.02) | (0.03) | (0.02) | (0.03) | (0.03) | (0.03)
Robust Tucker+K means | 0.42 | 0.41 | 0.41 | 0.4 | 0.4 | 0.39 | 0.39
(0.03) | (0.04) | (0.03) | (0.03) | (0.02) | (0.02) | (0.03)
LTRR | 1 | 0.94 | 0.86 | 0.74 | 0.58 | 0.54 | 0.53
(0) | (0.01) | (0.04) | (0.07) | (0.09) | (0.08) | (0.08)
TRR | 0.98 | 0.86 | 0.82 | 0.78 | 0.74 | 0.51 | 0.47
(0.01) | (0.05) | (0.08) | (0.11) | (0.09) | (0.06) | (0.07)
LRTSD | 1 | 1 | 1 | 1 | 1 | 1 | 0.99
(0) | (0) | (0) | (0) | (0) | (0) | (0.01)

###  4.4 Sensitivity Analysis

More specifically, we aim to answer two major questions by comparing the
proposed methods with several benchmark methods in terms of anomaly detection
accuracy and clustering accuracy.

  * •

Question 1: If sparse and gross anomalies are present, does anomaly detection
lead to better clustering?

  * •

Question 2: Does the tensor decomposition and dimension reduction lead to
better clustering?

####  4.4.1 The Effect of Dimensionality Reduction on Clustering Accuracy

We generate a simulated three-dimensional tensor with three clusters, thirty
samples per cluster, and 125 as the data dimensionality on each mode, namely,
$90\times 125\times 125$, where the intrinsic dimensionality is $5$ in each
cluster. Therefore, the total intrinsic dimensionality of the tensor is
actually $5\times 3=15$. This can be considered representative of a
spatiotemporal observation with one spatial dimension and two temporal
dimensions with the inherent clustered structure and redundant information
with respect to dimensionality.

For each seed, we generate some data and apply dimensionality reduction by a
certain factor on each mode. The dimension reduction factor
$\rho_{d}=[(\prod_{l}I_{l}/P_{l})^{1/l}]$ is defined as the geometric mean of
the data dimensionality $I_{l}$ over the intrinsic dimension $P_{l}$ averaged
on different mode $l$ of the tensor. For one factor of the experiment, we do
not apply any dimensionality reduction (which is denoted by a factor of 1 on
the $x$-axis in Figure 5). In other words, we directly apply subspace
clustering on the high-dimensional tensor, optimizing the following objective:

$\underset{\bm{Z}}{\text{min}}\quad\lambda\,\|\,\bm{Z}\,\|\,^{2}_{F}+\,\|\,\bm{X}_{1}-\bm{X}_{1}\bm{Z}\,\|\,_{F}^{2}$

Figure 5 shows the result of difference from minimum accuracy. For each
replication, we obtain the clustering accuracy for different dimensionality
reduction factors. Then they are subtracted from the minimum clustering
accuracy within that replication and get the difference for each
dimensionality reduction factor. After $20$ replications, we get the result of
$20$ accuracy differences for each dimensionality reduction factor. Their mean
is plotted in Figure 5. Observing Figure 5, there is an optimal range of
dimensionality reduction ratio around $8$ to $10$, which roughly corresponds
to the true intrinsic dimension ratio, i.e., $125/15=8.33$. This clearly
suggests that clustering is more effective at the right proportion of
dimensionality reduction, supporting our claim that dimension reduction
improves the tensor clustering performance.

####  4.4.2 The Effect of Anomaly Detection on Clustering Performance

The parameters for generating the simulation data are set to be the same as
those in Section 4.3.

We factorize our experiment over the properties of the added sparse anomaly.
Specifically, we alter how prevalent the anomaly is and how large its
intensity is. At each factor of the experiment, we define the control group as
the absence of anomaly detection, in other words, removing the term
$\frac{\lambda_{a}}{2}\,\|\,\mathcal{A}\,\|\,_{1}$ from Equation 2. The
treatment group is Equation 2 as is; therefore, the algorithm tries to detect
the sparse anomaly. We do this over multiple seeds, where in each seed the
exact same observations are used for both the control and treatment groups. We
report the difference in clustering accuracy between the treatment group and
the control group to identify whether anomaly detection would increase
clustering accuracy or not. The clustering accuracy is defined as almost the
same as the well-known clustering accuracy, except that the permutation of the
labels is allowed such that it maximizes the agreement between the predicted
labels and the ground truth.

Figure 5: The evolution of clustering accuracy over multiple replications by
increasing dimensionality factor on each mode.  Figure 6: The average of
pairwise differences in clustering accuracy between when anomaly detections is
allowed and when it is not, across various anomaly conditions and ten data
generating seeds.

In Figure 6, we observe a clear relationship between the severity of the
anomalies and how much anomaly detection improves clustering performance. As
the intensity of the anomaly increases and/or the prevalence of anomalies
increases, the application of anomaly detection becomes even more significant
in terms of finding the right groupings between the samples of the data. This
can be considered quite an expected finding, but it is an important sanity
check of our optimization procedure, and it finalizes the discussion of why
simultaneous anomaly detection is crucial for the quality of cluster
discovery.

##  5 Case Study & Results

We apply our proposed methodology to the smart card data collected by MTR, the
company that operates Hong Kong’s subway system. Smart card data emerged as an
important ingredient for urban mobility analysis, as it can be collected
cheaply and contains detailed information about how people move within a city
(Pelletier et al. 2011). Using this information, decisions at the tactical or
strategic level can be made in a data-driven, thus more effective manner.

The data span an 82-day time horizon from January 1 2017 and record every tap-
in and tap-out event from 85 stations of the system. We resample these events
hourly for each day within the operating hours of the system. Thus, we end up
with the number of people entering a station every operating hour of every
station.

Such data can be represented as a three-dimensional tensor. The first mode of
the tensor represents the stations, while the second and third modes represent
the day and the time of day, respectively. Thus, we expect to obtain clusters
at the station level using the clustering objective and point anomalies at
hourly resolution.

Understanding commonalities between stations is important from an urban
mobility perspective. The stations that are similar to each other in terms of
passenger inflow may reveal what kind of features about those stations are the
most influential in people’s decision to take the subway on a specific day and
time of day.

Detecting anomalies is also important from a managerial point of view.
Generally, urban mobility shows strong predictability due to the regular
weekly schedules of city dwellers. However, rare events could occur, and the
subway management would rather take precautions in the face of such events.
Detection of past point anomalies may reveal the underlying reasons for one-
off events that lead to extreme passenger inflow, such as concerts,
celebrations, or extreme weather events. Armed with this knowledge, subway
management can be more accurate in predicting future rare events.

###  5.1 Station Clustering

Figure 7: Eigenvalue spectrum of the affinity matrix to locate discontinuous
segments

We trained our model with the smart card data and recovered three clusters.
The number of clusters was determined using the method introduced in Von
Luxburg (2007). We first plot the logarithm of the eigenvalue spectrum of the
Laplacian of the affinity matrix, then locate the discontinuous segments, as
shown in Figure 7.

Figure 8: Hourly inflow patterns per cluster on a randomly picked week, across
weekdays and weekend days. The color coding matches with the stations in
Figure 9.  Figure 9: Hong Kong subway stations on a map, clustered three-way.
Colors represent different clusters. Circle radius represents a typical
walking distance of 400 meters around the station.

Figure 8 and Figure 9, when observed jointly, reveal important insight into
how the stations are clustered. Weekdays play the primary role in determining
the cluster of a station. As shown in the Figure 8, the first cluster (Cluster
ID: 0, colored in green) has a strong inflow peak at the morning peak hour.
These stations are also mainly located in the residential areas of Hong Kong
(i.e., west and east end of Hong Kong Island, west side of New Territories,
east side of Kowloon), as shown in the Figure 9, Likewise, for the second
cluster (Cluster ID: 1, colored in orange), we observe a strong evening peak
inflow as opposed to the mornings, and they are mostly located in the business
region such as the Central area of Hong Kong Island and the opposite side in
Kowloon, as shown in Figure 9. As we could conclude so far, for the most part,
this can be explained by the land use, specifically, whether it is residential
or work and office space. The third cluster (Cluster ID: 2, colored in purple)
does not have a dominant peak between the morning and evening, and the rest of
the inflow is rather evenly distributed throughout the day. While not as
prominent as weekdays, weekends (especially Saturday) patterns also differ
from one cluster to another.

Interpretation of the third cluster requires some contextual investigation.
Hong Kong draws millions of visitors each year, especially from the bordering
Shenzhen City, Guangdong Province of mainland China. The third cluster is a
collection of stations where visitors enter Hong Kong, with the primary
intention being cross-border shopping. The absence of taxes and duties on
goods in Hong Kong makes it an appealing destination for Chinese visitors to
shop. Since shopping is not a strictly time-bound activity, we observe inflows
that are dispersed more homogeneously around the day, both on weekdays and
weekends.

We now bring in the land use information 111Open data via
https://www.pland.gov.hk/pland_en/info_serv/open_data/landu/index.html to the
picture. In Figure 10, the first and second clusters can strictly be separated
by the prevalence of residential areas or commercial and industrial areas,
respectively. The only station from the third cluster, Sha Tin station, is
surrounded by many shopping malls, with an iconic New Town Plaza which offers
around 2 million sq. ft. of exceptional more than 400 shopping, dining, and
lifestyle facilities. Thus, it is a popular shopping destination for
passengers visiting from mainland China.

Figure 10: Sha Tin and Tai Wai regions with land use information. The circle
color coding of clusters matches with the stations in Figure 9. The land use
color coding logic is as follows: Red: residential, Green: commercial, Orange:
industrial, Brown: institutional, governmental, and Grey: parks and
recreations.

We can analyze even more thoroughly with granular information. The key
takeaway is that the clustering of stations yields important insight into the
factors that determine passenger inflows.

---
(a) Tsim Sha Tsui

---
(b) Hong Kong

---
(c) Admiralty

---
(d) Hung Hom

Figure 11: Examples of rare events in five different stations, detected as
point anomalies over passenger inflows.

###  5.2 Detected Point Anomalies

We now focus on the point anomalies detected by the proposed model in the data
from the smart card. As we have discussed before, one-off events are an
important determinant of irregular patterns in passenger inflows.

Lunar New Year celebrations in Hong Kong, which took place between 28 and 31
of January in Victoria Harbor between Hong Kong Island and the Kowloon
Peninsula, yield a number of other rare events. One of the events is the Lunar
New Year Fireworks Display. In Figure 11, we observe the peaks related to that
event in the Tsim Sha Tsui station, Admiralty and Hong Kong, with locations
pinned in Figure 9. These stations are close to places that are considered the
best places to watch the fireworks display.

Another interesting rare event takes place at the Hung Hom station. In Figure
11, we observed an unusual peak at 8 pm for four consecutive days from
February 16 to February 19. This suggests the scheduled end of a four-day
event, perhaps a conference, given the proximity of the station to Hong Kong
Polytechnic University.

In general, detected anomalies are promising in terms of the amount of insight
they can provide for public transportation management.

##  6 Conclusion

In this paper, we focused on the problem of spatio-temporal modeling with
tensors. Specifically, we proposed a tensor decomposition-based model that is
capable of simultaneous clustering and point anomaly detection, all of which
are important challenges to be addressed in this kind of data. Our simulation
study shows that, in general, point anomaly detection increases clustering
accuracy when point anomalies exist in data and vice versa. As for the case
study, we trained our model on a collection of smart card transactions from
the subway system of Hong Kong. Our model was able to effectively recover
station clusters that are semantically consistent both temporally and
spatially. The model was also able to recover rare events in terms of point
anomalies, which offers valuable insights and assistance for managerial
analysis and operations.

##  7 Appendix

###  7.1 Proof of Proposition 3.1

###### Proof 7.1

We rearrange the objective functions as follows:

$\displaystyle\mathcal{L}_{\bm{Z}}$
$\displaystyle\triangleq\frac{1}{2}\,\|\,\bm{Z}\,\|\,_{F}^{2}+\frac{\lambda_{z}}{2}\,\|\,\bm{C}_{(1)}-\bm{Z}\bm{C}_{(1)}\,\|\,_{F}^{2}$
$\displaystyle=\frac{1}{2}tr(\bm{Z}^{\top}\bm{Z})+\frac{\lambda_{z}}{2}tr(\bm{C}_{(1)}^{\top}\bm{C}_{(1)}-2\bm{C}_{(1)}^{\top}\bm{Z}\bm{C}_{(1)}+\bm{C}_{(1)}^{\top}\bm{Z}^{\top}\bm{Z}\bm{C}_{(1)})$
$\displaystyle=\frac{1}{2}tr(\bm{Z}^{\top}\bm{Z})+\frac{\lambda_{z}}{2}(tr(\bm{C}_{(1)}^{\top}\bm{C}_{(1)})-2tr(\bm{C}_{(1)}^{\top}\bm{Z}\bm{C}_{(1)})+tr(\bm{C}_{(1)}^{\top}\bm{Z}^{\top}\bm{Z}\bm{C}_{(1)}))$
$\displaystyle=\frac{1}{2}tr(\bm{Z}^{\top}\bm{Z})+\frac{\lambda_{z}}{2}(tr(\bm{C}_{(1)}^{\top}\bm{C}_{(1)})-2tr(\bm{C}_{(1)}^{\top}\bm{Z}\bm{C}_{(1)})+tr(\bm{C}_{(1)}\bm{C}_{(1)}^{\top}\bm{Z}^{\top}\bm{Z}))$

Taking the derivative of $\mathcal{L_{\bm{Z}}}$ with respect to $\bm{Z}$ we
get:

$\displaystyle\frac{\partial\mathcal{L}_{\bm{Z}}}{\partial\bm{Z}}$
$\displaystyle=\bm{Z}-\lambda_{z}(\bm{C}_{(1)}\bm{C}_{(1)}^{\top})+\lambda_{z}\bm{C}_{(1)}\bm{C}_{(1)}^{\top}\bm{Z}$

Setting $\frac{\partial\mathcal{L}_{\bm{Z}}}{\partial\bm{Z}}=\bm{0}$, we find
the closed-form solution $\bm{Z}^{*}$ as in Equation (6).

###  7.2 Proof of Proposition 3.3

###### Proof 7.2

Rearranging the terms, we get the following objective function:

$\displaystyle\mathcal{L}_{\bm{C}_{(1)}}$
$\displaystyle\triangleq\frac{\lambda_{z}}{2}\,\|\,\bm{C}_{(1)}-\bm{Z}\bm{C}_{(1)}\,\|\,_{F}^{2}+\frac{\lambda_{e}}{2}\,\|\,\bm{M}-\bm{C}_{(1)}\bm{P}\,\|\,_{F}^{2}$
$\displaystyle\begin{aligned}
=\frac{\lambda_{z}}{2}tr(\bm{C}_{(1)}^{\top}(\bm{I}-\bm{Z})^{\top}(\bm{I}-\bm{Z})\bm{C}_{(1)})+\\\
\frac{\lambda_{e}}{2}tr(\bm{M}^{\top}\bm{M}-2\bm{M}^{\top}\bm{C}_{(1)}\bm{P}+\bm{P}^{\top}\bm{C}_{(1)}^{\top}\bm{C}_{(1)}\bm{P})\end{aligned}$
$\displaystyle\begin{aligned}
=\frac{\lambda_{z}}{2}tr(\bm{C}_{(1)}^{\top}(\bm{I}-\bm{Z})^{\top}(\bm{I}-\bm{Z})\bm{C}_{(1)})+\frac{\lambda_{e}}{2}tr(\bm{M}^{\top}\bm{M})\\\
-\lambda_{e}tr(\bm{P}\bm{M}^{\top}\bm{C}_{(1)})+\frac{\lambda_{e}}{2}tr(\bm{P}\bm{P}^{\top}\bm{C}_{(1)}^{\top}\bm{C}_{(1)})\end{aligned}$

Taking the derivative of $\mathcal{L}_{\bm{C}_{(1)}}$ with respect to
$\bm{C}_{(1)}$ we get:

$\frac{\partial\mathcal{L}_{\bm{C}_{(1)}}}{\partial\bm{C}_{(1)}}=\lambda_{z}(\bm{I}-\bm{Z})^{\top}(\bm{I}-\bm{Z})\bm{C}_{(1)}-\lambda_{e}\bm{M}\bm{P}^{\top}+\lambda_{e}\bm{P}\bm{P}^{\top}\bm{C}_{(1)}$

Here, notice that
$\bm{P}\bm{P}^{\top}=(\bm{U}_{3}^{\top}\otimes\bm{U}_{2}^{\top})(\bm{U}_{3}\otimes\bm{U}_{2})=(\bm{U}_{3}^{\top}\bm{U}_{3})\otimes(\bm{U}_{2}^{\top}\bm{U}_{2})=\bm{I}$,
because of the properties of the Kronecker product and $\bm{U}_{2}$ and
$\bm{U}_{3}$ being orthonormal matrices. This helps us simplify the closed-
form solution $\bm{C}_{(1)}^{*}$ which is found by setting
$\frac{\partial\mathcal{L}_{\bm{C}_{(1)}}}{\partial\bm{C}_{(1)}}=\bm{0}$.

###  7.3 Proof of Proposition 3.4

###### Proof 7.3

By unfolding the tensor into the 2nd mode, we have

$\displaystyle\underset{\bm{U}_{2}}{\arg\min}\frac{\lambda_{e}}{2}\|\mathcal{X}-\mathcal{A}-\mathcal{C}\times_{2}\bm{U}_{2}\times_{3}\bm{U}_{3}\|_{F}^{2}$
$\displaystyle=$
$\displaystyle\underset{\bm{U}_{2}}{\arg\min}\|\bm{X}_{(2)}-\bm{A}_{(2)}-\bm{U}_{2}\bm{C}_{(2)}(\bm{U}_{3}\otimes\bm{I})^{T}\|_{F}^{2}$
$\displaystyle=$
$\displaystyle\underset{\bm{U}_{2}}{\arg\min}\|\bm{B}_{1}-\bm{U}_{2}\bm{B}_{2}\|_{F}^{2},$

where $\bm{B}_{1}=\bm{X}_{(2)}-\bm{A}_{(2)}$ and
$\bm{B}_{2}=\bm{C}_{(2)}(\bm{U}_{3}\otimes\bm{I})^{T}$. In this form, this
problem has the same structure as the orthogonal Procrustes problem
(Schönemann 1966).

The solution is given by the SVD decomposition on
$\bm{Q}=\bm{B}_{1}\bm{B}_{2}^{T}=(\bm{X}_{(2)}-\bm{A}_{(2)})(\bm{U}_{3}\otimes\bm{I})\bm{C}_{(2)}^{\top}$
as $\bm{Q}=\hat{\bm{U}}\bm{D}\hat{\bm{V}}^{\top}$. Finally, the solution is
$\bm{U}_{2}^{*}=\hat{\bm{U}}\hat{\bm{V}}^{\top}$.

This work is partially funded by the Hong Kong Research Grants Council (RGC)
under the grants RGC GRF 16201718 and 16216119.

## References

  * Aggarwal (2015) Aggarwal CC (2015) Outlier analysis. _Data mining_ , 237–263 (Springer).
  * Bahadori et al. (2014) Bahadori MT, Yu QR, Liu Y (2014) Fast multivariate spatio-temporal analysis via low rank tensor learning. _Advances in neural information processing systems_ , 3491–3499.
  * Chandola et al. (2009) Chandola V, Banerjee A, Kumar V (2009) Anomaly detection: A survey. _ACM computing surveys (CSUR)_ 41(3):1–58.
  * Drakopoulos et al. (2019) Drakopoulos G, Spyrou E, Mylonas P (2019) Tensor clustering: A review. _2019 14th International Workshop on Semantic and Social Media Adaptation and Personalization (SMAP)_ , 1–6 (IEEE).
  * Du and Zhang (2014) Du B, Zhang L (2014) A discriminative metric learning based anomaly detection method. _IEEE Transactions on Geoscience and Remote Sensing_ 52(11):6844–6857.
  * Fanaee-T and Gama (2016) Fanaee-T H, Gama J (2016) Tensor-based anomaly detection: An interdisciplinary survey. _Knowledge-Based Systems_ 98:130–147.
  * Fu et al. (2014) Fu Y, Gao J, Tien D, Lin Z (2014) Tensor lrr based subspace clustering. _2014 International Joint Conference on Neural Networks (IJCNN)_ , 1877–1884 (IEEE).
  * Fu et al. (2016) Fu Y, Gao J, Tien D, Lin Z, Hong X (2016) Tensor lrr and sparse coding-based subspace clustering. _IEEE transactions on neural networks and learning systems_ 27(10):2120–2133.
  * Gahrooei et al. (2021) Gahrooei MR, Yan H, Paynabar K, Shi J (2021) Multiple tensor-on-tensor regression: An approach for modeling processes with heterogeneous sources of data. _Technometrics_ 63(2):147–159.
  * Gao et al. (2015) Gao H, Nie F, Li X, Huang H (2015) Multi-view subspace clustering. _Proceedings of the IEEE international conference on computer vision_ , 4238–4246.
  * Guo et al. (2019) Guo K, Zhang T, Xu X, Xing X (2019) Low-rank tensor thresholding ridge regression. _IEEE Access_ 7:153761–153772.
  * He and Atia (2022) He Y, Atia GK (2022) Multi-mode tensor space clustering based on low-tensor-rank representation. _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 36, 6893–6901.
  * Heng et al. (2023) Heng Q, Chi EC, Liu Y (2023) Robust low-rank tensor decomposition with the l2 criterion. _Technometrics_ 65(4):537–552.
  * Hitchcock (1927) Hitchcock FL (1927) The expression of a tensor or a polyadic as a sum of products. _Journal of Mathematics and Physics_ 6(1-4):164–189.
  * Hu and Work (2021) Hu Y, Work DB (2021) Robust tensor recovery with fiber outliers for traffic events. _ACM Trans. Knowl. Discov. Data_ 15(1), ISSN 1556-4681.
  * Huang et al. (2008) Huang H, Ding C, Luo D, Li T (2008) Simultaneous tensor subspace selection and clustering: the equivalence of high order svd and k-means clustering. _Proceedings of the 14th ACM SIGKDD international conference on Knowledge Discovery and Data mining_ , 327–335.
  * Huang et al. (2023) Huang S, Zhang H, Zeng H, Pižurica A (2023) From model-based optimization algorithms to deep learning models for clustering hyperspectral images. _Remote Sensing_ 15(11):2832.
  * Jegelka et al. (2009) Jegelka S, Sra S, Banerjee A (2009) Approximation algorithms for tensor clustering. _International Conference on Algorithmic Learning Theory_ , 368–383 (Springer).
  * Kolda and Bader (2009) Kolda TG, Bader BW (2009) Tensor decompositions and applications. _SIAM review_ 51(3):455–500.
  * Li (2021) Li Z (2021) Tensor topic models with graphs and applications on individualized travel patterns. _2021 IEEE 37th International Conference on Data Engineering (ICDE)_ , 2756–2761 (IEEE).
  * Li et al. (2020a) Li Z, Sergin ND, Yan H, Zhang C, Tsung F (2020a) Tensor completion for weakly-dependent data on graph for metro passenger flow prediction. _proceedings of the AAAI conference on artificial intelligence_ , volume 34, 4804–4810.
  * Li et al. (2022a) Li Z, Yan H, Tsung F, Zhang K (2022a) Profile decomposition based hybrid transfer learning for cold-start data anomaly detection. _ACM Transactions on Knowledge Discovery from Data (TKDD)_ .
  * Li et al. (2020b) Li Z, Yan H, Zhang C, Tsung F (2020b) Long-short term spatiotemporal tensor prediction for passenger flow profile. _IEEE Robotics and Automation Letters_ 5(4):5010–5017.
  * Li et al. (2022b) Li Z, Yan H, Zhang C, Tsung F (2022b) Individualized passenger travel pattern multi-clustering based on graph regularized tensor latent dirichlet allocation. _Data Mining and Knowledge Discovery_ 1–32.
  * Lin et al. (2018) Lin C, Zhu Q, Guo S, Jin Z, Lin YR, Cao N (2018) Anomaly detection in spatiotemporal data via regularized non-negative tensor analysis. _Data Mining and Knowledge Discovery_ 32(4):1056–1073.
  * Liu et al. (2013) Liu G, Lin Z, Yan S, Sun J, Yu Y, Ma Y (2013) Robust recovery of subspace structures by low-rank representation. _IEEE Trans. Pattern Anal. Mach. Intell._ 35(1):171–184.
  * Liu et al. (2019) Liu G, Zhang Z, Liu Q, Xiong H (2019) Robust subspace clustering with compressed data. _IEEE Transactions on Image Processing_ 28(10):5161–5170.
  * Lock (2018) Lock EF (2018) Tensor-on-tensor regression. _Journal of Computational and Graphical Statistics_ 27(3):638–647.
  * Lu et al. (2008a) Lu H, Plataniotis KN, Venetsanopoulos AN (2008a) Mpca: Multilinear principal component analysis of tensor objects. _IEEE transactions on Neural Networks_ 19(1):18–39.
  * Lu et al. (2008b) Lu H, Plataniotis KN, Venetsanopoulos AN (2008b) Uncorrelated multilinear principal component analysis through successive variance maximization. _Proceedings of the 25th international conference on Machine learning_ , 616–623.
  * Mao et al. (2022) Mao Z, Li Z, Li D, Bai L, Zhao R (2022) Jointly contrastive representation learning on road network and trajectory. _arXiv preprint arXiv:2209.06389_ .
  * Nomikos and MacGregor (1994) Nomikos P, MacGregor JF (1994) Monitoring batch processes using multiway principal component analysis. _AIChE Journal_ 40(8):1361–1375.
  * Pan and Yang (2009) Pan SJ, Yang Q (2009) A survey on transfer learning. _IEEE Transactions on knowledge and data engineering_ 22(10):1345–1359.
  * Parsons et al. (2004) Parsons L, Haque E, Liu H (2004) Subspace clustering for high dimensional data: a review. _Acm sigkdd explorations newsletter_ 6(1):90–105.
  * Pelletier et al. (2011) Pelletier MP, Trépanier M, Morency C (2011) Smart card data use in public transit: A literature review. _Transportation Research Part C: Emerging Technologies_ 19(4):557–568.
  * Peng et al. (2015) Peng X, Yi Z, Tang H (2015) Robust subspace clustering via thresholding ridge regression. _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 29.
  * Shen et al. (2022) Shen B, Xie W, Kong ZJ (2022) Smooth robust tensor completion for background/foreground separation with missing pixels: Novel algorithm with convergence guarantee. _Journal of Machine Learning Research_ 23(217):1–40.
  * Sofuoglu and Aviyente (2021) Sofuoglu SE, Aviyente S (2021) Low-rank on graphs plus temporally smooth sparse decomposition for anomaly detection in spatiotemporal data. _ICASSP 2021-2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , 5614–5618 (IEEE).
  * Sun et al. (2022) Sun B, Zhao Z, Liu D, Gao X, Yu T (2022) Tensor decomposition-inspired convolutional autoencoders for hyperspectral anomaly detection. _IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing_ 15:4990–5000.
  * Sun and Li (2019) Sun WW, Li L (2019) Dynamic tensor clustering. _Journal of the American Statistical Association_ 114(528):1894–1907.
  * Tucker (1966) Tucker LR (1966) Some mathematical notes on three-mode factor analysis. _Psychometrika_ 31(3):279–311.
  * Vidal (2011) Vidal R (2011) Subspace clustering. _IEEE Signal Process. Mag._ 28(2):52–68.
  * Vidal and Favaro (2014) Vidal R, Favaro P (2014) Low rank subspace clustering (LRSC). _Pattern Recognit. Lett._ 43:47–61.
  * Vidal et al. (2008) Vidal R, Tron R, Hartley R (2008) Multiframe motion segmentation with missing data using powerfactorization and gpca. _International Journal of Computer Vision_ 79(1):85–105.
  * von Luxburg (2007) von Luxburg U (2007) A tutorial on spectral clustering. _Stat. Comput._ 17(4):395–416.
  * Von Luxburg (2007) Von Luxburg U (2007) A tutorial on spectral clustering. _Statistics and computing_ 17:395–416.
  * Wang et al. (2023) Wang M, Zhen Y, Pan Y, Xu Z, Guo R, Zhao X (2023) Tensorized hypergraph neural networks. _arXiv preprint arXiv:2306.02560_ .
  * Xue et al. (2017) Xue N, Papamakarios G, Bahri M, Panagakis Y, Zafeiriou S (2017) Robust low-rank tensor modelling using tucker and cp decomposition. _2017 25th European Signal Processing Conference (EUSIPCO)_ , 1185–1189 (IEEE).
  * Yan et al. (2015) Yan H, Paynabar K, Shi J (2015) Image-Based process monitoring using Low-Rank tensor decomposition. _IEEE Trans. Autom. Sci. Eng._ 12(1):216–227.
  * Yan et al. (2017) Yan H, Paynabar K, Shi J (2017) Anomaly detection in images with smooth background via smooth-sparse decomposition. _Technometrics_ 59(1):102–114.
  * Yan et al. (2018) Yan H, Paynabar K, Shi J (2018) Real-Time monitoring of High-Dimensional functional data streams via Spatio-Temporal smooth sparse decomposition. _Technometrics_ 60(2):181–197.
  * Yang et al. (2015) Yang C, Robinson D, Vidal R (2015) Sparse subspace clustering with missing entries. _International Conference on Machine Learning_ , 2463–2472.
  * Yang et al. (2019) Yang S, Wu J, Xu Y, Yang T (2019) Revealing heterogeneous spatiotemporal traffic flow patterns of urban road network via tensor decomposition-based clustering approach. _Physica A: Statistical Mechanics and its Applications_ 526:120688, ISSN 0378-4371.
  * Zhang et al. (2015) Zhang C, Fu H, Liu S, Liu G, Cao X (2015) Low-rank tensor constrained multiview subspace clustering. _Proceedings of the IEEE international conference on computer vision_ , 1582–1590.
  * Zhang et al. (2020) Zhang C, Yan H, Lee S, Shi J (2020) Dynamic multivariate functional data modeling via sparse subspace learning. _Technometrics_ 1–14.
  * Zhang et al. (2023) Zhang C, Zheng B, Tsung F (2023) Multi-view metro station clustering based on passenger flows: a functional data-edged network community detection approach. _Data Mining and Knowledge Discovery_ 37(3):1154–1208.
  * Zhao et al. (2022a) Zhao X, Dai Q, Wu J, Peng H, Liu M, Bai X, Tan J, Wang S, Philip SY (2022a) Multi-view tensor graph neural networks through reinforced aggregation. _IEEE Transactions on Knowledge and Data Engineering_ 35(4):4077–4091.
  * Zhao et al. (2022b) Zhao X, Yan H, Hu Z, Du D (2022b) Deep spatio-temporal sparse decomposition for trend prediction and anomaly detection in cardiac electrical conduction. _IISE Transactions on Healthcare Systems Engineering_ 12(2):150–164.
  * Zhao et al. (2021) Zhao Y, Yan H, Holte S, Mei Y (2021) Rapid detection of hot-spots via tensor decomposition with applications to crime rate data. _Journal of Applied Statistics_ 1–27.
  * Zhao et al. (2020) Zhao Y, Yan H, Holte SE, Kerani RP, Mei Y (2020) Rapid detection of hot-spot by tensor decomposition on space and circular time with application to weekly gonorrhea data. _THE XIIITH INTERNATIONAL WORKSHOP ON INTELLIGENT STATISTICAL QUALITY CONTROL 2019_.