# Fast Hypergraph Regularized Nonnegative Tensor Ring Factorization Based on
Low-Rank Approximation

Xinhai Zhao  Yuyuan Yu  Guoxu Zhou<EMAIL_ADDRESS>Qibin Zhao
<EMAIL_ADDRESS>Weijun Sun  School of Automation, Guangdong University of
Technology, Guangzhou, China  The Center for Advanced Intelligence Project
(AIP), RIKEN , Tokyo, Japan  Guangdong Key Laboratory of IoT Information
Technology, Guangzhou, China

###### Abstract

For the high dimensional data representation, nonnegative tensor ring (NTR)
decomposition equipped with manifold learning has become a promising model to
exploit the multi-dimensional structure and extract the feature from tensor
data. However, the existing methods such as graph regularized tensor ring
decomposition (GNTR) only models the pair-wise similarities of objects. For
tensor data with complex manifold structure, the graph can not exactly
construct similarity relationships. In this paper, in order to effectively
utilize the higher-dimensional and complicated similarities among objects, we
introduce hypergraph to the framework of NTR to further enhance the feature
extraction, upon which a hypergraph regularized nonnegative tensor ring
decomposition (HGNTR) method is developed. To reduce the computational
complexity and suppress the noise, we apply the low-rank approximation trick
to accelerate HGNTR (called LraHGNTR). Our experimental results show that
compared with other state-of-the-art algorithms, the proposed HGNTR and
LraHGNTR can achieve higher performance in clustering tasks, in addition,
LraHGNTR can greatly reduce running time without decreasing accuracy.

###### keywords:

Clustering , Low-rank approximation , Hypergraph , Tensor ring decomposition

††journal: Journal of LaTeX Templates

##  1 Introduction

Nowadays, extracting meaningful features from high-dimensional data has always
been an important topic in pattern recognition and machine learning.
Nonnegative matrix factorization (NMF) has been received plenty of attention
[1]. However, with the development of data mining or acquisition and multi-
sensor technology, more and more high-dimensional and large data are
generated. They also emphasize that the limitation of the flat-view matrix
models [2]. Tensors are the multi-dimensional extension of matrices, which can
preserve the structural information of high-order data. In recent years,
tensor ring decomposition (TR) has been proposed and demonstrated great
potentials in image completion and clustering [3, 4]. The tensor ring
decomposition essentially expresses the information of each dimension in
tensors in the form of a three-order core tensor. Moreover, the TR complexity
is increasing linearly with the dimension increasing and it can alleviate the
curse of dimension. In other words, TR is highly suitable for higher-order
tensor analysis. To learn the localized parts of tensor objects, Yu _et al._
proposed the non-negative tensor ring decomposition (NTR) [5].

Nonnegative tensor models have the abilities to discover local-parts based
information and take full advantages of the multi-linear structure of data.
Many researchers recently have found that data can be thought of as a low-
dimensional non-linear manifold embedded in high-dimensional data space [6,
7]. Because manifold learning can learn the intrinsic manifold structure of
data, nonnegative tensor decomposition methods that combine with manifold
learning showed better performance in many tasks and hence gained great
popularity in recent years. Cai _et al._ explored the geometrical information
by constructing the nearest neighbours graph and proposed the graph
regularized nonnegative matrix factorization (GNMF) [8]. However, the GNMF can
only perform feature extraction from vector data. Yu _et al._ proposed the
graph regularized non-negative tensor ring decomposition (GNTR) to learn the
localized parts of tensor objects and capture the manifold information[5]. But
the learned graph from GNTR only models the pairs of similarity relationships
between two objects and the manifold information of data can not be
sufficiently learned especially in the high-dimension.

There are two main problems in the above-mentioned methods: (1) they do not
take the character of higher-order manifold similarity relations of tensor
data into consideration, for example, in graphs, we can easily obtain the
similarity relations between two samples, but we can not predict if there
exist more than two samples having similarity. Therefore, the operation that
complex relationships among objects are simply squeezed as pair-wise forms
inevitably results in manifold information loss. (2) most of the manifold
methods adopt the Gaussian kernel to measure the distance between two samples.
It has been proved that the normal two-order graph may be easily affected by
the radius parameter $\sigma$ of the Gaussian kernel [9]. These two major
issues can largely be overcome by using hypergraph. The hypergraph is not
limited to learning only one connection relationship of both data samples, but
also the multi-level connection relationship among the multiple data samples
[10]. A graph is a special case of hypergraph with each edge containing two
vertices [11]. Zeng _et al._ incorporated the manifold regularized hypergraph
with the standard NMF framework. They proposed the hypergraph regularized
nonnegative matrix factorization (HNMF) [12]. But it is also limited in its
ability to discover the multi-linear structure of data. Yin _et al._ pointed
that contrasting the normal graph method, the hypergraph has the better
ability to extract the local and high-order geometry structure information. So
they combined the nonnegative Canonical Polyadic (CP) decomposition with
hypergraph and proposed the HyperNTF [13]. However, this model is difficult to
be optimized. And it also has a slow convergence rate.

When large volumes of high-dimensional data become increasingly common, it is
time-consuming obviously to analyze these data due to their high storage
requirements and high computational complexity. Fortunately, the meaningful
features usually also have relatively low-rank properties [14, 15]. A much
smaller number of potential variables and components can contain most of the
meaningful features of data. Therefore, low-rank approximation methods can be
used to preprocess the obtained data. Zhou _et al._ aiming at the problems of
slow convergence and excessive computational complexity, they introduced the
low-rank approximation to NMF and nonnegative Tucker decomposition (NTD) [16,
17].

In this paper, we propose a novel non-negative tensor ring decomposition
method to learn the local-parted features and capture the high-dimensional
manifold structure information simultaneously. The combination can learn the
manifold information completely and is significantly beneficial to improve the
clustering performance. This proposed method is called hypergraph regularized
non-negative tensor ring decomposition (HGNTR). Considering that the
computation of HGNTR has a large complexity, we introduce the method of low-
rank approximations into the HGNTR, then propose the LraHGNTR to reduce
running time and suppress noise. The main contributions of this article
include the following aspects:

  * (1)

HGNTR integrates the novel method of hypergraph with nonnegative tensor ring
decomposition, can effectively learn the local-parted attributes of tensor
data, and discover more high-dimensional manifold information to enhance the
efficiency of feature extraction in the clustering tasks.

  * (2)

We develop an effective iterative algorithm based on multiplicative updating
rules (MUR) to optimize the HGNTR model. Meanwhile, we apply the low-rank
approximation trick to significantly reduce the computational complexity of
HGNTR, which is called LraHGNTR.

  * (3)

Our experimental results show that both proposed algorithms can achieve higher
performance metrics in clustering tasks, in addition, LraHGNTR can greatly
reduce running time and even filter the noise from the original data.

The rest of this article is as follows: In Part 2, we firstly review the
concepts and basic knowledge. In Part 3, we propose the HGNTR method and
develop the MUR to update the model, and then develop the LraHGNTR. In Part 4,
we take the different kinds of experiments to valid our method. The
experimental results on the database prove the effectiveness of our theory.
Part 5 is followed by the conclusion.

##  2 Notations and Preliminaries

We simply review the related notations and the definitions of nonnegative
tensor decomposition as follows, and the basic notations are briefly reviewed
in Table 1.

Table 1: The notations

Notations | Descriptions | Notations | Descriptions
---|---|---|---
$x$ | A scalar | x | A vector
$\mathbf{X}$ | A matrix | $\mathcal{X}$ | A tensor
$\|\cdot\|_{F}$ | Frobenius norm | $\mathbf{I}$ | Identity matrix
$\operatorname{tr}\\{\cdot\\}$ | Trace operation | $\otimes$ | Kronecker product
$\times_{n}$ | Mode-$n$ product | $\overline{\times}^{1}$ | Multi-linear product

###  2.1 Definition

Definition 1 (Mode-$n$ unfolding) Given a tensor
$\mathcal{X}\in\mathbb{R}^{I_{1}\times I_{2}\times\cdots\times I_{N}}$, its
standard mode-$n$ unfolding matrix can be expressed as
$\mathbf{X}_{(n)}\in\mathbb{R}^{I_{n}\times I_{1}\cdots I_{n-1}I_{n+1}\cdots
I_{N}}$, by fixing all the indices except $I_{n}$ [18]. If we transforms the
unfolding matrix to the tensor, it is the inversion of the above operation and
can be denoted as follow:

$\mathcal{X}\leftarrow folding\left({\mathbf{X}_{(n)}}\right).$

The another mode-$n$ unfolding matrix of tensor $\mathcal{X}$ is denoted by
$\mathbf{X}_{[n]}\in\mathbb{R}^{I_{n}\times I_{n+1}\cdots I_{N}I_{1}\cdots
I_{n-1}}$, , by fixing all the indices except $I_{n}$, which is often used in
TR operations [19]. The distinction between two kinds of unfolding is the
sequence of the $N-1$ indices.

Definition 2 (Mode-$n$ product) The mode-$n$ product of a tensor
$\mathcal{X}\in\mathbb{R}^{I_{1}\times I_{2}\times\cdots\times I_{N}}$ with a
matrix $\mathbf{U}\in\mathbb{R}^{J_{n}\times I_{n}}$ is denoted by
$\mathcal{X}\times_{n}\mathbf{U}\in\mathbb{R}^{I_{1}\times\cdots\times
I_{n-1}\times J_{n}\times I_{n+1}\times\cdots\times I_{N}}$. See the [18] for
the details. Element-wisely, we have the following formula:

$\left(\mathcal{X}\times_{n}\mathbf{U}\right)_{i_{1}\cdots
i_{n-1}ji_{n+1}\cdots i_{N}}=\sum_{i_{n}=1}^{I_{n}}x_{i_{1}i_{2}\ldots
i_{N}}u_{ji_{n}}.$ (1)

Definition 3 (Multi-linear product) Given the tensor
$\mathcal{G}^{(n)}\in\mathbb{R}^{R_{n}\times I_{n}\times R_{n+1}}$ and
$\mathcal{G}^{(n+1)}\in\mathbb{R}^{R_{n+1}\times I_{n+1}\times R_{n+2}}$,
since the tensors $\mathcal{G}^{(n)}$ and $\mathcal{G}^{(n+1)}$ have an
equivalent mode size $R_{n+1}$, they can be integrated into a single tensor by
applying multi-linear product. And the single tensor can be denoted as
$\mathcal{G}^{(n,n+1)}\in\mathbb{R}^{R_{n}\times I_{n}I_{n+1}\times R_{n+2}}$
[20]. For multi-linear product, we have the following formula element-wisely:

$\mathcal{G}^{\left(n,n+1\right)}_{r_{n},\left(i_{n}-1\right)I_{n+1}+i_{n+1},r_{n+2}}=\sum_{r_{n+1}=1}^{R_{n+1}}\mathcal{G}^{\left(n\right)}_{r_{n},i_{n},r_{n+1}}\mathcal{G}^{(n+1)}_{r_{n+1},i_{n+1},r_{n+2}},$
(2)

where $i_{n}=1,2,\cdots,I_{n}$ and $i_{n+1}=1,2,\cdots,I_{n+1}$. For
simplicity, it can be rewritten as:

$\mathcal{G}^{\left(n,n+1\right)}=\mathcal{G}^{(n)}\overline{\times}^{1}\mathcal{G}^{(n+1)}.$
(3)

###  2.2 Nonnegative Tucker Decomposition

The nonnegative Tucker decomposition (NTD) make the tensor into the
multiplication of nonnegative core tensor and factor matrices in tensor’s each
mode [21]. It can be written as follow:

$\displaystyle\mathcal{A}=\mathcal{G}\times_{1}\mathbf{A}^{(1)}\times_{2}\mathbf{A}^{(2)}\cdots\times_{n}\mathbf{A}^{(N)},$
(4)

where $\mathcal{G}\in\mathbb{R}_{+}^{R_{1}\times R_{2}\times\cdots\times
R_{N}}$ denotes the nonnegative core tensor,
$\mathbf{A}^{(n)}\in\mathbb{R}_{+}^{I_{n}\times R_{n}},n=1,2,\cdots,N$ denote
the nonnegative factor matrices, $\mathbf{A}^{(n)}\geq 0$, and
$\mathcal{G}\geq 0$. The $\bm{R}=\left[R_{1},R_{2},\cdots,R_{N}\right]$ is the
nonnegative rank of NTD. However, the existence of the core tensor also
increases the computation complexity of the model and limits the ability to
represent higher-dimensional tensors.

###  2.3 Graph Regularized Nonnegative Tensor Ring Decomposition

Tensor ring (TR) decomposition is a more fundamental and general decomposition
model than the Canonical Polyadic (CP) decomposition [22] and the Tucker
decomposition (TD) [19]. The research results also proved the model’s
effectiveness in various applications, e.g, tensor-based image completion [23,
24], hyperspectral image denoising [20] , neural network compression [25, 26,
3], etc. Through the TR model, a high-order tensor
$\mathcal{X}\in\mathbb{R}^{I_{1}\times I_{2}\times\cdots\times I_{N}}$ can be
expressed as the circular contract products over a series of 3rd-order core
tensors. Each element of the tensor $\mathcal{X}$ can be written the as
follow:

$\mathcal{X}\left(i_{1},i_{2},\ldots,i_{N}\right)=\operatorname{Tr}\left\\{\prod_{n=1}^{N}\mathbf{G}_{n}(i_{n})\right\\},$
(5)

where $\mathbf{G}_{n}(i_{n})\in\mathbb{R}^{R_{n}\times R_{n+1}}$ is the
$i_{n}$-th lateral slice matrix of the core tensor
$\mathcal{G}^{(n)}\in\mathbb{R}^{R_{n}\times I_{n}\times R_{n+1}}$.
$\bm{R}=\left[R_{1},R_{2},\cdots,R_{N}\right]$ is called TR-rank. The TR
restricts the rank of the border core tensors to be equal, such as
$R_{1}=R_{N+1}$, not must be value-1, in other words, both border core tensors
are 3$rd$-order tensors. Comparing with the Tensor Train decomposition (TT)
[27] where the rank of edge core tensors $R_{1}=R_{N+1}$ must be restricted to
value-1 and it can only keep relatively limited connection and interaction
with between edge core tensors, the border core tensors of TR
$\mathcal{G}^{(1)}$ and $\mathcal{G}^{(N)}$ can fully connect or interact each
other directly enough. So TR can construct into a ring-like structure.
Meanwhile, it is well known that the trace operations of matrices have cyclic
invariance, such as
$\operatorname{Tr}\left(ABC\right)=\operatorname{Tr}\left(CAB\right)$, and
this property also holds for tensor ring decomposition. So the TR has the
circular dimensional permutation invariance. These special attributes can make
the TR model have some important numeric properties in the clustering task.

To further sufficiently learn local-parted information, Yu _et al._ proposed
the Nonnegative Tensor Ring decomposition (NTR) [5]. The NTR is based on the
TR model, it can be written as follow:

$\displaystyle\min$
$\displaystyle\frac{1}{2}\left\|\mathcal{X}-\operatorname{NTR}\left(\mathcal{G}^{(1)},\mathcal{G}^{(2)},\cdots,\mathcal{G}^{(N)}\right)\right\|^{2}_{F}$
(6) $\displaystyle\text{s.t.}\ \mathcal{G}^{(n)}\geq 0,\quad n=1,2,\cdots,N,$

where $\mathcal{G}^{(n)}\in\mathbb{R}_{+}^{R_{n}\times I_{n}\times R_{n+1}}$
denotes the $n$-th nonnegative core tensor.
$\operatorname{NTR}\left(\mathcal{G}^{(1)},\mathcal{G}^{(2)},\cdots,\mathcal{G}^{(N)}\right)$
can be thought of as the reconstruction of the core tensors. The Frobenius
norm is adopted to measure the similarity of them. According to the theorem
3.5 in [19], when we use alternating least squares (ALS) method, for the TR
decomposition of $\mathcal{X}$, its matrix-form can be written as
$\mathbf{X}_{[n]}=\mathbf{G}_{(2)}^{(n)}\left(\mathbf{G}_{\left[2\right]}^{\neq
n}\right)^{\top}$. The central idea of ALS is to optimize a core tensor when
the other core tensors are fixed, which has been extensively implemented in
the CP decomposition and TD. Thus, the NTR problem (6) can be rewritten as the
sub-problem framework:

$\displaystyle\min$
$\displaystyle\mathcal{F}_{NTR}^{(n)}=\frac{1}{2}\left\|\mathbf{X}_{[n]}-\mathbf{G}_{(2)}^{(n)}\left(\mathbf{G}_{\left[2\right]}^{\neq
n}\right)^{\top}\right\|_{F}^{2}$ (7) $\displaystyle\text{s.t.
}\mathbf{G}_{(2)}^{(n)}\geq 0,\quad n=1,2,\cdots,N,$

where $\mathbf{G}_{(2)}^{(n)}\in\mathbb{R}_{+}^{I_{n}\times R_{n}R_{n+1}}$
denotes the classical mode-$2$ unfolding matrix of $\mathcal{G}^{(n)}$.
$\bm{R}=\left[R_{1},R_{2},\cdots,R_{N}\right]$ denotes nonnegative rank of NTR
and $\mathbf{G}^{\neq n}_{\left[2\right]}\in\mathbb{R}_{+}^{I_{n+1}\cdots
I_{N}I_{1}\cdots I_{n-1}\times R_{n}R_{n+1}}$ denotes the another type of
mode-2 unfolding matrix of sub-chain tensor $\mathcal{G}^{\neq
n}\in\mathbb{R}_{+}^{R_{n+1}\times I_{n+1}\cdots I_{N}I_{1}\cdots I_{n}\times
R_{n}}$. The sub-chain tensor $\mathcal{G}^{\neq n}$ denotes the multi-linear
product of all core tensors except the $n$-th core tensor and can be written
as follow:

$\mathcal{G}^{\neq
n}=\mathcal{G}^{(n+1)}\overline{\times}^{1}\cdots\overline{\times}^{1}\mathcal{G}^{(N)}\overline{\times}^{1}\mathcal{G}^{(1)}\overline{\times}^{1}\cdots\overline{\times}^{1}\mathcal{G}^{(n-1)}.$

Detailed instructions about $\mathcal{G}^{\neq n}$ can be found in the
Definition 3.3 [19]. To enable NTR to observe the manifold geometrical
structure of high-dimensional data, Yu _et al._ propose a graph regularized
nonnegative tensor ring decomposition (GNTR) [5]. However, the learned graph
form GNTR only models the single relationship between two objects, and
therefore does not accurately indicate the high-dimensional manifold
information.

###  2.4 Hypergraph

The idea of hypergraph manifold learning is inspired by the theory of simple
two-dimensional graph [13, 12]. In a normal graph, an edge is only connected
with two vertices and the weight of the edge merely indicates the relationship
between two connected vertices. However, in reality, the relationship of
objects is more complex, there is a great possibility that there is a close
connection between three vertices or more than. In this situation, the normal
graph can not deal with this complex data well. The main distinction between
hypergraph and normal two-order graph is that the hypergraph uses a subset of
the vertices as an edge and each edge of a hypergraph usually connects more
than two vertices. So the method of hypergraph manifold learning can promote
the clustering performance [28, 29, 30, 31, 32]. $Fig.\ref{fig:hypergraph}$ is
an example that shows the difference between hypergraph and normal graph.

(a) (b)

Figure 1: The difference between the hypergraph and normal graph. (a) A normal
graph. (b) A hypergraph consists of the nodes and lines that denotes their
relationships. (c) The incidence matrix $\mathbf{H}$ of the hypergraph.

We assume that $\mathbf{G}=\left(V,E,\mathbf{W}\right)$ is a weighted
hypergraph. $V$ is a finite set of vertices in the hypergraph. $E$ is the set
of hyperedges. And each hyperedge $e$ is a subset of $V$, then $\cup_{e\in
E}=V$. $\mathbf{W}$ is a diagonal matrix that represents the weights of
hyperedge in which its elements are $w$.

The relationship of nodes and hyperedges in the hypergraph can be expressed by
an indicator matrix $\mathbf{H}$ with the size of
$\left|V\right|\times\left|E\right|$. And the incidence matrix $\mathbf{H}$ of
hypergraph $\mathbf{G}$ is defined as follows element-wisely [13]:

$h\left(v,e\right)=\left\\{\begin{matrix}1&if\>\>\>v\in e\\\ 0&if\>\>\>v\notin
e.\\\ \end{matrix}\right.$ (8)

The affinity matrix $\mathbf{A}$ of hypergraph is calculated as follow:

$\mathbf{A}_{ij}=\mathbf{exp}\left(-\frac{\left\|v_{i}-v_{j}\right\|^{2}}{\sigma^{2}}\right).$
(9)

where $\sigma$ is the average distance, $v_{i}$ and$v_{j}$ denotes two
vertices in the hypergraph. It denotes the probabilistic hypergraph [33]. It
is important to note that since hypergraph can establish multiple connections
among multiple vertices, they can learn a more stable and compact
representation than normal graph. This special characteristic causes the
hypergraph to be insensitive to hyper-parameters $\sigma$. And the weight
matrix $\mathbf{W}$ of hypergraph is calculated as follow:

$\mathbf{W}_{i}=\sum_{v_{j}\in e_{i}}A\left(i,j\right).$ (10)

The degree of hyperedge $e$ with a symbol of $d_{E}\left(e\right)$, measures
the amount of vertices or nodes $v\in V$, is defined as
$d_{E}\left(e\right)=\begin{matrix}\sum_{v\in
V}\end{matrix}h\left(v,e\right)$. The degree of each vertex $v\in V$ is the
sum of weights for the hyperedges, is denoted as
$d_{V}\left(v\right)=\begin{matrix}\sum_{e\in
E}\end{matrix}w\left(e\right)h\left(v,e\right)$. Here we denote
$\mathbf{D}_{V}$ as a diagonal matrix where its elements is
$d_{V}\left(v\right)$. And $\mathbf{D}_{E}$ can be denoted by a diagonal
matrix where its elements is $d_{E}\left(e\right)$.

After obtaining a hypergraph that encodes the relationships among the samples,
the data analysis tasks can be fulfilled by grouping the vertices in a way
that measures the degree of similarity between samples, also known as
hypergraph cutting missions from a graph learning perspective [11]. Then we
can have the following hypergraph structure loss function:

$\displaystyle\Omega=\frac{1}{2}\sum_{e\in E}\sum_{v_{i},v_{j}\in
e_{i}}\frac{w(e)h(v_{i},e)h(v_{j},e)}{d_{E}(e)}\|\frac{\mathbf{F}(v_{i})}{\sqrt{d(v_{i})}}-\frac{\mathbf{F}(v_{j})}{\sqrt{d(v_{j})}}\|_{F}^{2},$
(11)

where $\mathbf{F}=\left[f_{1},\cdots,f_{m}\right]$ denotes the $n\times m$
matrix. And $f_{i}$ is a $n$-dimensional vector where each element indicates
the probability that a vertex is a part of sub-graph. The Eq.(11) is designed
to learn a set of hypergraph cuts that retain as many hyperedges as possible
when partitioning. Notably, the sub-graphs are clusters of data after
grouping. Thus, the complex data geometrical information can be preserved in
the space of learned hypergraph cuts with much lower dimensionality than one
of the sample spaces [34]. After a few steps of calculation [11], the Eq.(11)
can be converted into a matrix form, as follows:

$\displaystyle\Omega=$
$\displaystyle\operatorname{tr}\left\\{\mathbf{F}^{T}(\mathbf{D}_{v}-\mathbf{H}\mathbf{W}\mathbf{D}_{E}^{-1}\mathbf{H}^{T})\mathbf{F}\right\\}$
(12) $\displaystyle=$
$\displaystyle\operatorname{tr}\left\\{\mathbf{F}^{T}\mathbf{L}_{hyper}\mathbf{F}\right\\},$

where $\mathbf{L}_{hyper}=\mathbf{D}_{v}-\mathbf{S}$, and
$\mathbf{S}=\mathbf{H}\mathbf{W}\mathbf{D}_{E}^{-1}\mathbf{H}^{T}$.

##  3 Hypergraph Regularized Nonnegative Tensor Ring Factorization

In this part, in order to capture the multi-linear structural information and
high-dimensional geometric manifold information of the data simultaneously, we
innovatively add a hypergraph regularization term to the NTR model and propose
the HGNTR model based on NTR decomposition. This model can express the
original tensor data as a series of non-negative tensor circular contract
products and completely capture the high-dimensional and complex manifold
structural information. And we develop an effective algorithm to optimize the
HGNTR model. At the same time, considering that the meaningful part of the
data usually has a certain low-rank property, to reduce the computational
complexity significantly, then we propose a low-rank approximation method
(LraHGNTR).

###  3.1 Objective Function of HGNTR

Here we propose the hypergraph regularized NTR model, called HGNTR, to
factorize the nonnegative tensors. And our goal is to use the hypergraph to
learn the manifold information of the last dimension. It is also be said that
to learn the manifold structure of the last dimensional unfolding matrix if we
write the HGNTR in the sub-problem framework. The reason is that in tensor
data of our experimental database, the last dimension of data usually denotes
the samples and defaults the number of objects, and the hypergraph regularized
term is naturally imposed on the $N$-th core tensor. Given a tensor
$\mathcal{X}\in\mathbb{R}^{I_{1}\times I_{2}\times\cdots\times I_{N}}$, the
objective function of HGNTR is as follows:

$\displaystyle\min_{\mathbf{G}_{(2)}^{(n)}}\mathcal{F}_{HGNTR}^{(n)}=$
$\displaystyle\frac{1}{2}\left\|\mathbf{X}_{[n]}-\mathbf{G}_{(2)}^{(n)}\left(\mathbf{G}_{\left[2\right]}^{\neq
n}\right)^{\top}\right\|_{F}^{2}+\frac{\beta}{2}\operatorname{tr}\left(\left(\mathbf{G}_{(2)}^{(N)}\right)^{\top}\mathbf{L}_{hyper}\mathbf{G}_{(2)}^{(N)}\right)$
(13) $\displaystyle\text{s.t. }\mathbf{G}^{(n)}_{(2)}\geq 0,\quad
n=1,2,\cdots,N,$

in which the $\mathbf{L}_{hyper}\in\mathbb{R}^{I_{N}\times I_{N}}$ is the
hypergraph Laplacian matrices; $\lambda\geq 0$ is a trade-off parameter to
control the punishment of regularization term. $Fig.\ref{fig:tensor_orl}$ and
$Fig.\ref{fig:tensor_gt}$ explain the principle and mechanism of HGNTR
decomposition on ORL and GT database. The brief introduction of the databases
is presented in Section 4.2. The last core tensor can be thought of as the
feature tensor, and the other $N-1$ core tensors can be thought of as the
basic tensors. By applying hypergraph regularization to the last core tensor,
the higher-order similarity relations are preserved in the feature tensors
that are decomposed by HGNTR.

Figure 2:  Illustration picture of HGNTR on three-order database. HGNTR
decomposes the tensor into three core tensor, and expresses the information of
each dimension as corresponding core tensor. Figure 3:  Illustration picture
of HGNTR on four-order database. HGNTR decomposes the tensor into four core
tensor, and expresses the information of each dimension as corresponding core
tensor.

###  3.2 Optimization of HGNTR

To solve the learned objective function in Eq.(13), we adopt the method of
Multiplicative Updating Rules (MUR). Noticing that the equation
$\left\|\mathbf{A}\right\|_{F}^{2}=\operatorname{Tr}(\mathbf{A}\mathbf{A}^{\top})$,
we define $\mathbf{a}_{ik}$ as the Langrange multiplier for constraints
$\mathbf{G}^{(n)}_{(2)}\geq 0$. And $\mathbf{A}=\left[\mathbf{a}_{ik}\right]$.
The Langrange function $\mathcal{L}_{HGNTR}^{(n)}$ is as follow:

$\displaystyle\mathcal{L}_{HGNTR}^{(n)}$
$\displaystyle=\frac{1}{2}\operatorname{tr}\left(\mathbf{X}_{[n]}\left(\mathbf{X}_{[n]}\right)^{\top}-2\mathbf{G}_{(2)}^{(n)}\left(\mathbf{G}_{\left[2\right]}^{\neq
n}\right)^{\top}\left(\mathbf{X}_{[n]}\right)^{\top}+\mathbf{G}_{(2)}^{(n)}\left(\mathbf{G}_{\left[2\right]}^{\neq
n}\right)^{\top}\mathbf{G}_{\left[2\right]}^{\neq
n}\left(\mathbf{G}_{\left(2\right)}^{\left(n\right)}\right)^{\top}\right)$
(14)
$\displaystyle+\frac{\beta}{2}\operatorname{tr}\left(\left(\mathbf{G}_{(2)}^{(n)}\right)^{\top}\mathbf{L}_{hyper}\mathbf{G}_{(2)}^{(n)}\right)+\operatorname{Tr}\left(\mathbf{A}\left(\mathbf{G}_{(2)}^{(n)}\right)^{T}\right).$

In the case of $n=N$, the derivatives of $\mathcal{L}_{HGNTR}^{(n)}$ with
respect to $\mathbf{G}^{(n)}_{(2)}$ is as follow :

$\displaystyle\frac{\partial\mathcal{L}_{HGNTR}^{(n)}}{\partial\mathbf{G}_{(2)}^{(n)}}=-\mathbf{X}_{[n]}\mathbf{G}_{\left[2\right]}^{\neq
n}+\mathbf{G}_{(2)}^{(n)}\left(\mathbf{G}_{\left[2\right]}^{\neq
n}\right)^{\top}\mathbf{G}_{\left[2\right]}^{\neq
n}+\beta\left(\mathbf{G}_{(2)}^{(n)}\right)^{\top}\mathbf{L}_{hyper}+\mathbf{A}.$
(15)

By using Karush-Kuhn-Tucker (KKT) condition, element-wisely
$\mathbf{a}_{ik}\left(\mathbf{G}_{(2)}^{(n)}\right)_{ik}=0$, we can get the
follow update rules for each element in case of $n=N$:

$\displaystyle\left(\mathbf{G}_{(2)}^{(n)}\right)_{ik}\leftarrow\left(\mathbf{G}_{(2)}^{(n)}\right)_{ik}\frac{\left(\mathbf{X}_{[n]}\mathbf{G}_{\left[2\right]}^{\neq
n}+\beta\left(\mathbf{G}_{\left(2\right)}^{(n)}\right)^{\top}\mathbf{S}\right)_{ik}}{\left(\mathbf{G}_{(2)}^{(n)}\left(\mathbf{G}_{\left[2\right]}^{\neq
n}\right)^{\top}\mathbf{G}_{\left[2\right]}^{\neq
n}+\beta\left(\mathbf{G}_{\left(2\right)}^{(n)}\right)^{\top}\mathbf{D}_{v}\right)_{ik}}.$
(16)

While $n=1,2,\cdots,N-1$, we just need to ignore the hypergraph regularized
term. It is rewritten as follow:

$\displaystyle\left(\mathbf{G}_{(2)}^{(n)}\right)_{ik}\leftarrow\left(\mathbf{G}_{(2)}^{(n)}\right)_{ik}\frac{\left(\mathbf{X}_{[n]}\mathbf{G}_{\left[2\right]}^{\neq
n}\right)_{ik}}{\left(\mathbf{G}_{(2)}^{(n)}\left(\mathbf{G}_{\left[2\right]}^{\neq
n}\right)^{\top}\mathbf{G}_{\left[2\right]}^{\neq n}\right)_{ik}}.$ (17)

The whole algorithm procedure is given in Algorithm 1. Convergence can be
proved by constructing auxiliary functions. The method is the same as the
convergence proof in [8]. The valid experiment is designed to show the
convergence property of HGNTR and LraHGNTR in section 4.7.

Algorithm 1 HGNTR based on MU method

0: Tensor $\mathcal{X}\in\mathbb{R}^{I_{1}\times I_{2}\times\cdots\times
I_{N}}$, nonnegative multi-way rank $r_{n}$ for $n=1,2,\cdots,N$, maximum
number of iterations of $t_{max}$, balance parameter $\beta$, the number of
nearest neighbours $k$.

0: Core tensors $\mathcal{G}_{n}$ for $n=1,2,\cdots,N$.

1: Initialize $\mathcal{G}_{n}\in\mathbb{R}^{R_{n}\times I_{n}\times R_{n+1}}$
for $n=1,2,\cdots,N$ as random tensors from the uniform distribution between 0
and 1.

2: Calculate hypergraph Laplacian matrix $\mathbf{L}_{hyper}$ by
$\mathbf{L}_{hyper}=\mathbf{D}_{V}-\mathbf{S}$.

3: repeat

4: for $n$ = 1 to $N$ do

5: Obtain mode-2 unfolding matrix of sub-chain $\mathbf{G}_{[2]}^{\neq n}$.

6: for $t$ = 1 to $t_{max}$ do

7: if $n=N$ then

8: Updating $\mathbf{G}_{(2)}^{(n)}$ and by Eq.(16).

9: else

10: Updating $\mathbf{G}_{(2)}^{(n)}$ and by Eq.(17) .

11: end if

12: end for

13: $\mathcal{G}_{n}\leftarrow$
folding$\left({\mathbf{G}_{(2)}^{(n)}}^{t_{max}}\right)$.

14: end for

15: until convergence

###  3.3 Fast HGNTR Algorithm Based on Low-rank Approximation

In Eq.(16) and Eq.(17), $\mathbf{X}_{[n]}\mathbf{G}_{\left[2\right]}^{\neq n}$
results in a large amount of computational complexity, especially in the
condition of large-scale data. If $\mathbf{X}_{[n]}$ can be replaced by a
smaller term, it is possible to reduce the running time by a significant
amount. Considering the well ability of Tucker decomposition (TD) [18] to
approximate tensor data, and the fact that it also can decompose high-
dimensional tensor into a series of low-rank factors, we choose the Tucker
decomposition (TD) for the low-rank approximation.

The low-rank approximation has been successfully used in the field of matrix
decomposition [35, 36, 37]. In the field of tensor, the fast tensor
decomposition methods are proposed in [38] , which are based on fiber
sampling. Sequence extraction based on tSVD is also proposed in [39], which
can achieve a balance between accuracy and efficiency. In addition, the theory
of low-rank approximation has been applied to the NMF and NTF [16, 17], and
the experiment results have demonstrated the effectiveness in reducing
complexity and denoising. These methods can be used for dimensionality
reduction of tensor data.

In this paper, we approximately change the original tensor $\mathcal{X}$ into
$\tilde{\mathcal{X}}$ by the TD. Then the low-rank approximate term is used in
iterative processes of HGNTR. We call this fast HGNTR algorithm the low-rank
HGNTR (LraHGNTR).Details about denoising can be found in the section 4.5. The
LraHGNTR can be explained as follows in two steps:

The low-rank approximation (LRA) step: obtain the LRA result of $\mathcal{X}$
by normal Tucker decomposition (TD). It is denoted as follow:

$\displaystyle\mathcal{X}\approx\mathcal{\tilde{X}}$
$\displaystyle=\left\llbracket\mathcal{C};\mathbf{U}^{(1)},\mathbf{U}^{(2)},\cdots,\mathbf{U}^{(N)}\right\rrbracket$
(18)
$\displaystyle=\mathcal{C}\times_{1}\mathbf{U}^{(1)}\times_{2}\mathbf{U}^{(2)}\cdots\times_{N}\mathbf{U}^{(N)}.$

the
$\mathcal{C}_{+}\in\mathbb{R}_{+}^{\tilde{R_{1}}\times\tilde{R_{2}}\times\cdots\times\tilde{R_{N}}}$
denotes Tucker core tensor, and
$\mathbf{U}^{(n)}\in\mathbb{R}_{+}^{{I_{n}\times\tilde{R_{n}}}},n=1,2,\cdots,N$
denote factor matrices. Here,
$\left(\mathbf{U}^{n}\right)^{\top}\mathbf{U}^{n}=\mathbf{I}$. And
$\bm{\tilde{R}}=\left[\tilde{R}_{1},\tilde{R}_{2},\cdots,\tilde{R}_{N}\right]$
is used to control the degree and error of approximation. The effect of rank
$\bm{\tilde{R}}$ on of the approximation accuracy will be described in detail
in Section 4.4.

The NTR step: we need to perform NTR to minimize
$\frac{1}{2}\left\|\tilde{\mathcal{X}}-\operatorname{NTR}\left(\mathcal{G}^{(1)},\mathcal{G}^{(2)},\cdots\mathcal{G}^{(N)}\right)\right\|_{F}^{2}$.
Now, we show how the derivatives computing of $\mathcal{L}_{HGNTR}^{(n)}$ on
$\mathbf{G}^{(n)}_{(2)}$ can be accelerated after $\mathcal{X}$ is replaced
with its low-rank approximation $\tilde{\mathcal{X}}$. If the precondition
exists, which tensor $\mathcal{X}$ has not only the structure of Tucker, but
also the structure of NTR, we can have the following formula:

$\displaystyle\tilde{\mathcal{X}}$
$\displaystyle=~{}\mathcal{C}\times_{1}\mathbf{U}^{(1)}\times_{2}\mathbf{U}^{(2)}\cdots\times_{N}\mathbf{U}^{(N)}$
(19)
$\displaystyle\approx\operatorname{NTR}\left(\mathcal{G}^{(1)},\mathcal{G}^{(2)},\cdots\mathcal{G}^{(N)}\right),$

where $\mathcal{G}^{(n)}\in\mathbb{R}^{R_{n}\times I_{n}\times R_{n+1}}$.
Considering $\left(\mathbf{U}^{n}\right)^{\top}\mathbf{U}^{n}=\mathbf{I}$ ,we
can get the following:

$\displaystyle\mathcal{C}\approx\operatorname{NTR}\left(\mathcal{G}^{(1)},\mathcal{G}^{(2)},\cdots\mathcal{G}^{(N)}\right)\times_{1}\left(\mathbf{U}^{(1)}\right)^{T}\times_{2}\left(\mathbf{U}^{(2)}\right)^{T}\times_{3}\cdots\times_{N}\left(\mathbf{U}^{(N)}\right)^{T}.$
(20)

By storing $\mathcal{C}$, $\mathcal{G}$ and $\mathbf{U}$ to reconstruct the
approximation $\tilde{\mathcal{X}}$, we can consume much less computational
complexity and storage especially when $\tilde{R_{n}}\ll I_{n}$. By the
Theorem 4.2 in the [19], Eq.(20) can be rewritten as follow:

$\displaystyle\mathcal{C}\times_{n}\mathbf{U}^{(n)}$ (21)
$\displaystyle\approx$
$\displaystyle\operatorname{NTR}\left(\mathcal{G}^{(1)}\times_{2}\left(\mathbf{U}^{(1)}\right)^{T},\cdots,\mathcal{G}^{(n-1)}\times_{2}\left(\mathbf{U}^{(n-1)}\right)^{T},\mathcal{G}^{n},\mathcal{G}^{(n+1)}\times_{2}\left(\mathbf{U}^{(n+1)}\right)^{T},\cdots\mathcal{G}^{(N)}\times_{2}\left(\mathbf{U}^{(N)}\right)^{T}\right)$
$\displaystyle=$
$\displaystyle\operatorname{NTR}\left(\mathcal{Z}^{(1)},\cdots,\mathcal{Z}^{(n-1)},\mathcal{G}^{(n)},\mathcal{Z}^{(n+1)},\cdots,\mathcal{Z}^{N}\right),$

where
$\mathcal{Z}^{(n)}=\mathcal{G}^{(n)}\times_{2}\left(\mathbf{U}^{(n)}\right)^{T}\in\mathbb{R}^{R_{n}\times\tilde{R}_{n}\times
R_{n+1}},~{}n=1,\cdots,n-1,n+1,\cdots,N$. The term
$\mathbf{X}_{[n]}\mathbf{G}_{\left[2\right]}^{\neq n}$ can be written as
follow:

$\displaystyle~{}\mathbf{X}_{[n]}\mathbf{G}_{[2]}^{\neq n}$ (22)
$\displaystyle\approx$
$\displaystyle~{}\mathbf{U}^{(n)}\mathbf{C}_{(n)}\left(\mathbf{U}^{\neq
n}\right)^{T}\mathbf{G}^{\neq n}_{[2]}$ $\displaystyle=$
$\displaystyle~{}\mathbf{U}^{(n)}\mathbf{C}_{(n)}\mathbf{Z}^{\neq n}_{[2]}.$

where
$\mathbf{C}_{(n)}\in\mathbb{R}^{\tilde{R}_{n}\times\tilde{R}_{1}\cdots\tilde{R}_{n-1}\tilde{R}_{n+1}\cdots\tilde{R}_{N}}$
is the mode-n unfolding matrix of $\mathcal{C}$ and $\mathbf{Z}^{\neq
n}_{\left[2\right]}\in\mathbb{R}^{\tilde{R}_{n+1}\cdots\tilde{R}_{N}\tilde{R}_{1}\cdots\tilde{R}_{n-1}\times
R_{n}R_{n+1}}$ denotes the mode-2 unfolding matrix of sub-chain tensor
$\mathcal{Z}^{\neq n}$. $\mathcal{Z}^{\neq n}$ and $\mathbf{U}^{\neq n}$ is
defined as follow:

$\mathcal{Z}^{\neq
n}=\mathcal{Z}^{(n+1)}\overline{\times}^{1}\cdots\overline{\times}^{1}\mathcal{Z}^{(N)}\overline{\times}^{1}\mathcal{Z}^{(1)}\overline{\times}^{1}\cdots\overline{\times}^{1}\mathcal{Z}^{(n-1)},$
$\mathbf{U}^{\neq
n}=\left[\mathbf{U}^{(N)}\otimes\cdots\otimes\mathbf{U}^{(n+1)}\otimes\mathbf{U}^{(n-1)}\otimes\cdots\otimes\mathbf{U}^{(1)}\right]^{T}.$

Finally, we need to change the $\mathbf{X}_{[n]}\mathbf{G}_{[2]}^{\neq n}$ by
Eq.(22). And we can get the following update rules for each element in case of
$n=N$:

$\displaystyle\left(\mathbf{G}_{(2)}^{(n)}\right)_{ik}\leftarrow\left(\mathbf{G}_{(2)}^{(n)}\right)_{ik}\frac{\left(\mathbf{U}^{(n)}\mathbf{C}_{(n)}\mathbf{Z}_{[2]}^{\neq
n}+\beta\left(\mathbf{G}_{\left(2\right)}^{(n)}\right)^{\top}\mathbf{S}\right)_{ik}}{\left(\mathbf{G}_{(2)}^{(n)}\left(\mathbf{G}_{\left[2\right]}^{\neq
n}\right)^{\top}\mathbf{G}_{\left[2\right]}^{\neq
n}+\beta\left(\mathbf{G}_{\left(2\right)}^{(n)}\right)^{\top}\mathbf{D}_{v}\right)_{ik}}.$
(23)

While $n=1,2,\cdots,N-1$, the update rules are rewritten as follow:

$\displaystyle\left(\mathbf{G}_{(2)}^{(n)}\right)_{ik}\leftarrow\left(\mathbf{G}_{(2)}^{(n)}\right)_{ik}\frac{\left(\mathbf{U}^{(n)}\mathbf{C}_{(n)}\mathbf{Z}_{[2]}^{\neq
n}\right)_{ik}}{\left(\mathbf{G}_{(2)}^{(n)}\left(\mathbf{G}_{\left[2\right]}^{\neq
n}\right)^{\top}\mathbf{G}_{\left[2\right]}^{\neq n}\right)_{ik}}.$ (24)

The whole LraHGNTR algorithm is given in Algorithm 2.

Algorithm 2 LraHGNTR based on MU method

0: Tensor $\mathcal{X}\in\mathbb{R}^{I_{1}\times I_{2}\times\cdots\times
I_{N}}$, nonnegative multi-way rank $r_{n}$ for $n=1,2,\cdots,N$, multi-linear
rank $\tilde{\mathbf{R}}_{n}$, maximum number of iterations of $t_{max}$,
balance parameter $\beta$, the number of nearest neighbours $k$.

0: Core tensors $\mathcal{G}_{n}$ for $n=1,2,\cdots,N$.

1: Initialize $\mathcal{G}_{n}\in\mathbb{R}^{R_{n}\times I_{n}\times R_{n+1}}$
for $n=1,2,\cdots,N$ as random tensors from the uniform distribution between 0
and 1.

2: Calculate hypergraph Laplacian matrix $\mathbf{L}_{hyper}$ by
$\mathbf{L}_{hyper}=\mathbf{D}_{V}-\mathbf{S}$.

3: repeat

4: for $n$ = 1 to $N$ do

5: Obtain the sub-chain tensor $\mathcal{Z}^{\neq n}$.

6: for $t$ = 1 to $t_{max}$ do

7: if $n=N$ then

8: Updating $\mathbf{G}_{(2)}^{(n)}$ by Eq.(23).

9: else

10: Updating $\mathbf{G}_{(2)}^{(n)}$ by Eq.(24).

11: end if

12: end for

13: $\mathcal{G}_{n}\leftarrow$
folding$\left({\mathbf{G}_{(2)}^{(n)}}^{t_{max}}\right)$.

14: end for

15: until convergence.

###  3.4 Computational Complexity Analysis

The HGNTR is efficiently updated by the method of Multiplicative Updating
Rule. To simplify, we assume that $I=I_{n}$ and $R=R_{n}R_{n+1}\,$,
$n=1,2,\cdots,N$. For the each iteration of the HGNTR algorithm, the tensor
$\mathcal{G}^{\neq n}$ is calculated and its computational cost is
$\mathcal{O}_{\mathcal{G}^{\neq
n}}=(I^{2}+I^{3}+\cdots+I^{N-1})R^{\frac{3}{2}}$. It is important to note that
we do not need to compute $\left(\mathbf{G}_{\left[2\right]}^{\neq
n}\right)^{\top}\mathbf{G}_{\left[2\right]}^{\neq n}$ explicitly, instead, we
just use the self-contraction operation for each small core tensors
$\mathcal{G}^{(n)},n=1,2,\cdots,n=N-1$ in sequential order, then compute their
product. The similar approaches has been used for CP and Tucker decomposition.
By this operation, the computational complexity of
$(\mathbf{G}_{\left[2\right]}^{\neq
n})^{\top}\mathbf{G}_{\left[2\right]}^{\neq n}$ is changing from
$(I^{N-1})R^{2}$ to $(N-1)IR^{2}+(N-2)R^{3}$. We can conclude that as the
dimension increases, the computational complexity increases linearly, making
this operation more useful to represent higher-order tensor data such as four-
order or more.

The total computational complexity of HGNTR in case of $n=1,2,\cdots,N-1$ is
$\mathcal{O}_{\mathcal{G}^{\neq
n}}+\mathcal{O}((I^{N}R+IR+IR^{2}+\left(N-1\right)IR^{2}+\left(N-2\right)R^{3}))$.
In the case of $n=N$, the HGNTR only adds the cost of calculating the graph
regularization term. As a result of the $\mathbf{L}_{hyper}$ is usually
sparse, the computational complexity of the hypergraph term is small. Then the
computational cost of HGNTR is $\mathcal{O}_{\mathcal{G}^{\neq
n}}+\mathcal{O}(I^{N}R+IR+IR^{2}+\left(N-1\right)IR^{2}+\left(N-2\right)R^{3}+2I^{2}R)$.

For the LraHGNTR, we assume that $I=I_{n}$ and
$\tilde{R}=\tilde{R}_{n}\tilde{R}_{n+1},n=1,2,\cdots,N$. The computing
complexity of the mode-2 product of core tensors and factor matrices in
Eq.(21) is $\mathcal{O}_{M}=\left(N-1\right)\tilde{R}RI$. The total
computational complexity of LraHGNTR in case of $n=1,2,\cdots,N-1$ is
$\mathcal{O}_{M}+\mathcal{O}(I(\tilde{R}^{N-1}(R+\tilde{R}))+IR^{2}+IR+(N-1)IR^{2}+(N-2)R^{3})$.
In the case of $n=N$, the total computational cost of LraHGNTR is
$\mathcal{O}_{M}+\mathcal{O}(I(\tilde{R}^{N-1}(R+\tilde{R}))+IR^{2}+2I^{2}R+IR+(N-1)IR^{2}+(N-2)R^{3})$.

The Table 2 shows the total floating-point counts of HGNTR and LraHGNTR in
detail when updating. From the above complexity analysis, we can conclude
that, when the data size is large, and $R=\tilde{R}_{n}\ll I$, the
computational complexity is greatly reduced and the speed of LraHGNTR is
significantly ahead.

Algorithm |  HGNTR |  LraHGNTR
---|---|---
Addition |  $NI^{N}R+NIR^{2}+$ |  $NI(\tilde{R}^{N-1}(R+\tilde{R}))+2IR$
$N(N-1)IR^{2}+2IR$ |  $NIR^{2}+N(N-1)IR^{2}$
$2I^{2}R+N(N-2)R^{3}$ |  $+N(N-2)R^{3}+2I^{2}R$
Multiplication |  $NI^{N}R+NIR^{2}+NIR+$ |  $NI(\tilde{R}^{N-1}(R+\tilde{R}))+NIR+$
$+N(N-1)IR^{2}$ |  $NIR^{2}+N(N-1)IR^{2}$
$+2I^{2}R+N(N-2)R^{3}$ |  $+N(N-2)R^{3}+2I^{2}R$
Division |  $NIR$ |  $NIR$
Overall |  $\mathcal{O}\left(NI^{N}R\right)$ |  $\mathcal{O}\left(NI(\tilde{R}^{N-1}(R+\tilde{R}))\right)$
Table 2: Computational operation floating-point counts of HGNTR and LraHGNTR
for each iteration when updating.

##  4 Experiment

In this section, we firstly introduce performance metrics such as accuracy
(ACC), normalized mutual information (NMI), and purity (PUR). Then we simply
review the databases that are used for our experiment. We select six proposed
algorithms as the comparison in the experiments to prove our model
effectiveness. We do clustering experiments to count the performance metrics
of our algorithm on different datasets. After that, we test the accuracy of
the noisy-added database. Finally, we will do the parameter sensitivity
experiment and convergence experiment to test the robustness and stability of
HGNTR and LraHGNTR. And we make the basic visualization to verify our
excellent performance of clustering again.

###  4.1 Performance Metrics

To evaluate the effectiveness of each algorithm, we adopt three metrics such
as ACC, NMI, and PUR. The definition of ACC is as follows:

$\operatorname{ACC}\left(y_{i},\hat{y}_{i}\right)=\frac{1}{n}\sum_{i=1}^{n}\delta\left(y_{i},\operatorname{map}\left(\hat{y}_{i}\right)\right),$

where $n$ is the total amount of objects of database. $y_{i}$ and
$\hat{y}_{i}$ respectively denotes the cluster label of the object and the
true label of the object. $\operatorname{map}\left(\cdot\right)$ denotes the
permutation mapping function, which is responsible for mapping each cluster
label $y_{i}$ to the equivalent label of the database. $\delta\left(\right)$
is the delta function. The AC is used to measure how similar the true label
and the experiment-got label . For example, if the object label $y_{i}$ and
the real label $\hat{y}_{i}$ are equal, then
$\left(y_{i},\operatorname{map}\left(\hat{y}_{i}\right)\right)=1$, if not,
then $\left(y_{i},\operatorname{map}\left(\hat{y}_{i}\right)\right)=0$.

By combining with the information theory, the agreement between two cluster
parts can be measured with mutual information (MI). The MI is widely used in
clustering applications. MI between the two sets of cluster labels
$C^{\prime}$ and true labels $C$ is defined by

$\operatorname{MI}\left(C,C^{\prime}\right)=\sum_{c_{i}\in C,c_{i}^{\prime}\in
C^{\prime}}p\left(c_{i},c_{i}^{\prime}\right)\cdot\log_{2}\frac{p\left(c_{i},c_{i}^{\prime}\right)}{p\left(c_{i}\right)\cdot
p\left(c_{i}^{\prime}\right)},$

in which $p\left(c_{i}\right)$ and $p\left(c_{i}^{\prime}\right)$ denote the
probability that the object that is selected arbitrarily from the databases
belongs to of category $c_{i}$ and category $c_{i}^{\prime}$. And
$p\left(c_{i},c_{i}^{\prime}\right)$ is defined as the joint probability that
object belongs to category $c_{i}$ and category $c_{i}^{\prime}$ at the same
time. To force the score to have an upper bound, we use the NMI as one of
evaluation measures and the definition of NMI is denoted as follows:

$\operatorname{NMI}\left(C,C^{\prime}\right)=\frac{\operatorname{MI}\left(C,C^{\prime}\right)}{\max\left(H(C),H\left(C^{\prime}\right)\right)},$

where $H\left(C\right)$ and $H\left(C^{\prime}\right)$ are defined as the
entropies of the true label set $C$ and the cluster label set $C^{\prime}$. It
is obvious to know the score ranges of
$\operatorname{NMI}\left(C,C^{\prime}\right)$ is from $0$ to $1$. If the two
label sets are the same, the $\operatorname{NMI}\left(C,C^{\prime}\right)=1$,
otherwise $\operatorname{NMI}\left(C,C^{\prime}\right)=0$ .

The PUR measures the degree that how many each cluster contains data points
are from single category. Given a certain category $C_{i}$ which size is
$n_{i}$, the definition of PUR is as follows [40]:

$\operatorname{PUR}(C_{i})=\frac{1}{n_{i}}\max_{h}\left(n_{i}^{h}\right),$

where $\max_{h}\left(n_{i}^{h}\right)$ is the number that represents the
number of categories from category $C_{i}$ belongs to category $h$.

###  4.2 Databases

In the experiment, we use two databases to prove the effectiveness of our
method. A brief introduction of the two databases is as follows:

  * 1.

ORL Database: The ORL face database contains 400 grey-scale images of 40
different individuals. For the image in each category, these images are
collected in a different light, facial expressions, and facial details. All
the images are taken against a dark, uniform background, and the front face.
In this experiment, each image is adjusted to the size of $32\times 27$. We
construct a 3rd-order tensor $\mathcal{X}\in\mathbb{R}^{32\times 27\times
400}$.

  * 2.

GT Database: The Georgia Tech face database contains images of 50 people. All
people in the database are represented by 15 color images with a cluttered
background that is taken at resolution $640\times 480$ pixels. The average
size of the faces in these images is 150x150 pixels. The pictures show frontal
or tilted faces with different facial expressions, lighting conditions, and
scales. We downsampled every image to $40\times 30\times 3$, finally we get
the 4rd-order tensor $\mathcal{X}\in\mathbb{R}^{40\times 30\times 3\times
750}$.

  * 3.

Face95 Database: The Face95 dataset consists of color facial images,
represented by 72 individuals, each with 20 exclusive images. Hence, there are
1440 face images in the dataset. The dataset contains frontal faces with mild
expressions. We downsampled the database to $50\times 40\times 3$, finally we
get the 4rd-order tensor $\mathcal{X}\in\mathbb{R}^{50\times 40\times 3\times
1000}$.

  * 4.

Face96 Database: The facial images in the Face96 dataset are full color and
consist of 152 individuals, each with 20 images. Therefore, there are 3040
images in the dataset. All faces in this dataset have complicated backgrounds
comparing to ORL and Face95. Furthermore, the images were taken under various
lighting situations and the faces contain different ranges of expressions. We
downsampled every image to $49\times 49\times 3$, finally we get the 4rd-order
tensor $\mathcal{X}\in\mathbb{R}^{49\times 49\times 3\times 1200}$.

###  4.3 Comparison Algorithms

In order to ensure the fairness of the experiment, our algorithms all use
Multiplicative Updating (MU) rules to update the model. To verify the
effectiveness of our proposed algorithms, we compare our algorithms with the
following state-of-the-art algorithms in our tasks.

  * 1.

K-means: One of the most famous unsupervised clustering algorithms.

  * 2.

NMF: Nonnegative matrix factorization aims to learn the parts-based basis of
data objects [1].

  * 3.

GNMF: Graph regularized nonnegative matrix factorization (GNMF), the NMF
considering the manifold geometric structure in the data [8].

  * 4.

HNMF: Hypergraph regularized nonnegative matrix factorization (HNMF), the NMF
considering the high-dimensional manifold geometric structure in the data
[12].

  * 5.

NTR-MU: Nonnegative tensor ring decomposition based on MU optimization rule to
learn the parts-based basis of tensor data objects [5].

  * 6.

GNTR-MU: Graph regularized nonnegative tensor ring decomposition based on MU
optimization rule, the NTR considering the manifold geometric structure in the
data [5].

###  4.4 Clustering Tasks for Real-world Database

In this part of the experiment, we test the comprehensive clustering
performance of extracted low-dimensional features of our proposed algorithms
HGNTR and LraHGNTR on databases, including their subsets. The method of
getting subsets is that we randomly select the partial $n$ samples from each
database and merge $n$ samples as a subset of each database. Considering the
number of categories, we divide the ORL into
$k_{ORL}\in\\{10,15,20,25,30,35,40\\}$ categories. We divide the GT into
$k_{GT}\in\\{10,20,30,40,50\\}$ categories. We divide the Face95 into
$k_{Face95}\in\\{10,20,30,40,50\\}$ categories. We divide the Face96 into
$k_{Face96}\in\\{20,30,40,50,60\\}$ categories. Meanwhile, while running the
matrix methods such as NMF, GNMF, HNMF, we need to transform the tensor into
the matrix by the Matlab reshape operation. All the experiments are performed
on a personal computer with an i5-9500T at 2.20GHz, 16GB memory, Windows 10,
and Matlab 2020b. For the ORL, the maximum number of iterations is set to 10.
For the other databases, maximum number of iterations is set to 5. We have
counted the ACC, NMI, PUR, and their STDs obtained from the experimental
results in Table 3, Table 4, Table 5, and Table 6. Through the experimental
results, we can get the following conclusions:

  * -

For the matrix method, the effectiveness of HNMF is usually better than that
of GNMF and NMF, indicating that the manifold information obtained by our
hypergraph regularization can better reflect the complex manifold geometry
information of the data. And it also can help improve the recognizability of
low-dimensional character. Meanwhile, we have also observed that with the
growth of sample categories, in other words, the amount of data, HNMF shows
limited ability to obtain hypergraph information. It proves the necessity that
we need to develop the tensor method of manifold learning that is more
suitable for high-dimensional data.

  * -

GNTR is usually better than NTR, shows that the learned manifold geometry can
effectively enhance the identifiability of low-dimensional features. HGNTR is
usually better than HNMF, indicating that our algorithm HGNTR can learn low-
dimensional features better by learning spatial structure information in the
form of tensors.

  * -

HGNTR is usually better than general other algorithms, which illustrates the
model in which tensor learning algorithm and hypergraph are combined together,
can learn the complex structure information, comparing with the state-of-the-
art methods. Hypergraph is an effective method, it can find more than two
connection relationships between data points. It can improve the accuracy of
clustering.

  * -

When the data volume (number of categories) of the three datasets is small, we
find that sometimes the performance metrics of tensor methods are a little
weaker than matrix methods. The reason is that tensor methods are more
effective in the condition of larger amounts of data. It also reflects the
advantages of the tensor approaches.

We also count the running time of each database under different categories in
$Fig.\ref{fig:runtime1}$. Through the experimental results, we can get the
following conclusions:

  * -

Noting that the question of running time of HGNTR and the computational
complexity, we consider using a low-rank Tucker approximation method to
replace the greatly complex term in the iterative process. Experiments have
shown that the performance of HGNTR and LraHGNTR are similar, but LraHGNTR can
greatly reduce the running time. It verifies that the low-rank Tucker
approximation can effectively reduce the computational cost while retaining
the accuracy .This is a key point of our algorithm leading other algorithms.

  * -

Through the comparison of subsets, we also found that as the category of
subsets increases, in other words, as the number of data increases, compared
with other algorithms, our algorithm HGNTR and LraHGNTR outperforms more
obviously. Especially for LraHGNTR, as the amount of data increases, the
performance of our algorithm LraHGNTR drops significantly in the running time.

By the experiment, we can conclude that the effect of rank $\bm{\tilde{R}}$
for LraHGNTR, if we choose a lower rank, our running time will be greatly
reduced, but the accuracy will be lost. If we choose a little larger rank, our
accuracy will increase. So if we choose a suitable rank, we can retain the
accuracy, and get the ideal running time.

| $k_{ORL}$ | K-means | NMF | GNMF | HNMF | NTR-MU | GNTR-MU | HGNTR-MU | lraHGNTR-MU
---|---|---|---|---|---|---|---|---|---
ACC | 10 | 75.80$\pm$3.35 | 77.70$\pm$3.26 | 74.60$\pm$2.87 | 75.30$\pm$3.40 | 71.70$\pm$2.75 | 78.00$\pm$2.70 | 78.30$\pm$3.49 |  78.10$\pm$3.14
15 | 69.93$\pm$2.97 | 75.06$\pm$5.57 | 73.80$\pm$4.48 | 76.33$\pm$4.01 | 68.53$\pm$4.80 | 75.40$\pm$2.40 | 77.06$\pm$3.37 |  76.66$\pm$2.96
20 | 69.60$\pm$3.31 | 72.55$\pm$2.99 | 75.15$\pm$2.69 | 78.85$\pm$3.85 | 70.00$\pm$2.69 | 78.00$\pm$4.18 |  80.05$\pm$2.87 | 80.10$\pm$2.99
25 | 68.64$\pm$4.24 | 72.80$\pm$4.58 | 76.12$\pm$2.34 | 75.90$\pm$3.05 | 68.60$\pm$3.08 | 77.64$\pm$2.66 |  80.68$\pm$3.69 | 81.20$\pm$2.95
30 | 73.36$\pm$2.66 | 69.36$\pm$2.66 | 73.80$\pm$2.41 | 72.96$\pm$3.39 | 66.16$\pm$4.91 | 73.60$\pm$2.66 |  74.00$\pm$2.51 | 76.06$\pm$2.21
35 | 66.54$\pm$4.69 | 69.08$\pm$2.43 | 72.85$\pm$2.30 | 72.32$\pm$2.15 | 63.71$\pm$2.97 | 73.65$\pm$2.54 | 74.08$\pm$1.60 |  73.91$\pm$2.87
40 | 66.95$\pm$3.54 | 67.93$\pm$2.01 | 73.28$\pm$2.77 | 71.47$\pm$2.73 | 68.55$\pm$3.43 | 74.70$\pm$2.85 |  75.65$\pm$2.09 | 76.55$\pm$2.10
NMI | 10 | 85.54$\pm$1.06 | 84.19$\pm$2.24 | 84.71$\pm$1.32 | 84.78$\pm$0.82 | 80.74$\pm$3.43 | 84.07$\pm$1.31 |  85.40$\pm$1.06 | 85.22$\pm$0.78
15 | 84.13$\pm$1.47 | 85.30$\pm$3.55 | 85.59$\pm$1.99 |  88.11$\pm$1.65 | 80.79$\pm$1.76 | 86.22$\pm$1.30 | 87.94$\pm$1.52 | 87.93$\pm$1.15
20 | 82.99$\pm$2.08 | 84.03$\pm$2.21 | 85.60$\pm$1.86 | 88.12$\pm$2.46 | 83.05$\pm$2.04 | 87.36$\pm$1.62 |  89.15$\pm$0.95 | 89.18$\pm$1.72
25 | 84.62$\pm$1.71 | 85.40$\pm$2.75 | 87.38$\pm$1.67 | 87.81$\pm$1.05 | 82.79$\pm$1.78 | 87.35$\pm$0.80 |  89.98$\pm$1.20 | 90.06$\pm$0.96
30 | 84.13$\pm$1.77 | 84.46$\pm$1.55 | 87.20$\pm$0.76 | 87.02$\pm$1.35 | 82.78$\pm$2.37 | 87.36$\pm$0.99 |  88.20$\pm$1.25 | 88.76$\pm$0.84
35 | 83.87$\pm$2.17 | 85.13$\pm$1.18 | 86.55$\pm$1.42 | 87.13$\pm$1.21 | 80.89$\pm$1.61 | 87.02$\pm$1.15 |  88.06$\pm$1.01 | 88.11$\pm$0.91
40 | 84.59$\pm$1.43 | 84.08$\pm$0.68 | 87.78$\pm$1.07 | 87.44$\pm$0.98 | 84.78$\pm$1.65 | 88.51$\pm$0.77 |  89.81$\pm$1.07 | 89.91$\pm$0.73
PUR | 10 | 80.10$\pm$2.60 | 80.70$\pm$2.45 | 78.90$\pm$1.59 | 79.90$\pm$2.96 | 75.30$\pm$3.26 | 81.30$\pm$1.41 |  81.60$\pm$2.01 | 81.70$\pm$2.05
15 | 76.06$\pm$2.34 | 79.33$\pm$4.71 | 79.33$\pm$3.23 | 81.40$\pm$2.72 | 72.20$\pm$2.86 | 80.53$\pm$1.62 |  81.66$\pm$2.49 | 81.73$\pm$1.87
20 | 73.30$\pm$3.09 | 76.10$\pm$3.06 | 78.25$\pm$2.72 | 81.55$\pm$3.70 | 73.85$\pm$2.82 | 80.75$\pm$3.29 |  82.65$\pm$1.90 | 82.50$\pm$2.83
25 | 73.64$\pm$3.36 | 76.40$\pm$3.77 | 79.20$\pm$1.98 | 79.04$\pm$1.96 | 72.12$\pm$2.83 | 80.12$\pm$2.00 |  82.88$\pm$2.69 | 83.28$\pm$1.70
30 | 71.60$\pm$3.30 | 73.36$\pm$2.66 | 76.96$\pm$2.07 | 76.43$\pm$2.72 | 74.23$\pm$3.45 | 77.20$\pm$1.91 |  77.46$\pm$2.05 | 78.96$\pm$1.95
35 | 70.97$\pm$3.46 | 72.94$\pm$2.28 | 75.83$\pm$2.32 | 76.17$\pm$1.64 | 67.45$\pm$2.92 | 76.54$\pm$1.96 |  76.74$\pm$2.12 | 77.45$\pm$2.07
40 | 71.33$\pm$2.95 | 71.53$\pm$0.76 | 76.83$\pm$2.03 | 75.47$\pm$1.78 | 72.95$\pm$2.87 | 78.30$\pm$1.92 |  79.45$\pm$1.52 | 79.65$\pm$2.00
Table 3: Performance metrics comparisons on ORL database between K-means, NMF, GNMF, HNMF, NTR-MU, GNTR-MU, HGNTR-MU and LraHGNTR-MU methods. We divide the database into different subset. Each subset contains $k_{ORL}$ categories.  | $k_{GT}$ | K-means | NMF | GNMF | HNMF | NTR-MU | GNTR-MU | HGNTR-MU | lraHGNTR-MU
---|---|---|---|---|---|---|---|---|---
ACC | 10.00 | 62.66$\pm$2.49 | 66.40$\pm$3.07 | 62.13$\pm$3.31 | 67.33$\pm$2.44 | 62.66$\pm$3.36 |  74.80$\pm$4.09 | 74.00$\pm$2.49 | 75.33$\pm$2.90
20.00 | 53.20$\pm$4.97 | 54.80$\pm$1.82 | 46.53$\pm$4.18 | 54.33$\pm$0.62 | 57.33$\pm$4.23 | 61.26$\pm$2.03 |  61.53$\pm$0.93 | 62.66$\pm$3.30
30.00 | 46.48$\pm$2.67 | 48.27$\pm$1.70 | 43.91$\pm$2.24 | 49.06$\pm$0.92 | 48.53$\pm$6.14 | 53.78$\pm$2.29 |  58.22$\pm$2.48 | 59.02$\pm$1.81
40.00 | 46.47$\pm$0.85 | 47.26$\pm$2.57 | 45.76$\pm$0.38 | 48.50$\pm$1.56 | 53.90$\pm$8.34 | 52.70$\pm$1.03 |  55.37$\pm$2.49 | 55.80$\pm$1.55
50.00 | 45.71$\pm$0.76 | 45.72$\pm$1.16 | 43.33$\pm$1.80 | 49.49$\pm$1.20 | 43.47$\pm$6.21 | 50.80$\pm$2.86 |  52.00$\pm$2.37 | 53.25$\pm$2.33
NMI | 10.00 | 72.76$\pm$2.74 | 72.97$\pm$2.30 | 68.06$\pm$2.74 | 69.00$\pm$3.09 | 66.45$\pm$1.77 | 74.25$\pm$3.79 |  74.69$\pm$2.70 | 75.55$\pm$1.64
20.00 | 66.99$\pm$3.01 | 65.97$\pm$1.67 | 58.30$\pm$3.96 | 66.81$\pm$1.17 | 67.42$\pm$3.68 | 71.89$\pm$1.97 |  73.01$\pm$0.87 | 73.04$\pm$2.10
30.00 | 62.21$\pm$0.92 | 63.50$\pm$2.17 | 59.33$\pm$1.63 | 63.57$\pm$0.80 | 64.45$\pm$5.05 | 67.70$\pm$0.62 |  70.85$\pm$1.32 | 71.51$\pm$0.91
40.00 | 63.26$\pm$0.83 | 63.54$\pm$1.39 | 62.98$\pm$0.84 | 65.06$\pm$0.93 | 68.69$\pm$5.55 | 68.71$\pm$0.41 | 71.24$\pm$0.76 |  70.82$\pm$0.65
50.00 | 64.37$\pm$0.58 | 63.64$\pm$0.48 | 62.76$\pm$1.58 | 66.12$\pm$0.86 | 63.30$\pm$5.52 | 69.25$\pm$2.31 |  70.64$\pm$2.23 | 71.42$\pm$2.30
PUR | 10.00 | 66.26$\pm$3.60 | 69.60$\pm$2.64 | 64.93$\pm$2.38 | 67.46$\pm$2.68 | 64.66$\pm$3.62 |  74.80$\pm$4.09 | 74.53$\pm$2.23 | 76.13$\pm$1.52
20.00 | 56.87$\pm$3.73 | 57.47$\pm$1.72 | 49.33$\pm$3.74 | 57.73$\pm$1.48 | 60.67$\pm$3.66 | 64.60$\pm$2.01 |  64.86$\pm$0.98 | 65.06$\pm$3.18
30.00 | 49.38$\pm$2.07 | 50.89$\pm$2.07 | 46.67$\pm$1.84 | 52.35$\pm$1.06 | 52.08$\pm$6.55 | 57.37$\pm$2.30 |  61.11$\pm$1.72 | 61.86$\pm$1.09
40.00 | 49.20$\pm$0.43 | 50.00$\pm$2.57 | 48.03$\pm$0.55 | 51.33$\pm$1.75 | 56.73$\pm$7.16 | 55.93$\pm$0.75 |  58.36$\pm$2.07 | 59.10$\pm$1.13
50.00 | 48.24$\pm$0.92 | 47.78$\pm$1.26 | 46.48$\pm$2.02 | 51.65$\pm$1.23 | 46.43$\pm$6.79 | 54.49$\pm$2.98 |  55.49$\pm$2.38 | 56.16$\pm$2.98
Table 4: Performance metrics comparisons on GT database between K-means, NMF, GNMF, HNMF, NTR-MU, GNTR-MU, HGNTR-MU and LraHGNTR-MU methods. We divide the database into different subset. Each subset contains $k_{GT}$ categories.  | $k_{Face95}$ | K-means | NMF | GNMF | HNMF | NTR-MU | GNTR-MU | HGNTR-MU | lraHGNTR-MU
---|---|---|---|---|---|---|---|---|---
ACC | 10 | 76.50$\pm$6.76 | 79.20$\pm$4.20 | 84.90$\pm$1.51 | 82.40$\pm$5.09 | 73.40$\pm$4.76 | 85.20$\pm$2.41 |  85.60$\pm$2.10 | 86.10$\pm$2.56
20 | 59.05$\pm$1.67 | 58.10$\pm$2.82 | 65.50$\pm$1.74 | 60.40$\pm$3.92 | 56.45$\pm$1.68 | 64.05$\pm$4.66 |  64.80$\pm$4.42 | 66.00$\pm$0.53
30 | 50.96$\pm$3.47 | 49.53$\pm$1.70 | 54.73$\pm$1.69 | 55.60$\pm$1.42 | 51.66$\pm$1.53 | 55.06$\pm$1.78 | 56.50$\pm$1.87 |  56.30$\pm$2.62
40 | 47.90$\pm$3.10 | 45.70$\pm$1.90 | 50.00$\pm$1.40 | 49.25$\pm$2.06 | 46.02$\pm$1.61 | 50.17$\pm$1.80 |  50.58$\pm$1.96 | 51.43$\pm$1.83
50 | 44.82$\pm$1.99 | 44.18$\pm$1.38 | 49.24$\pm$1.20 | 46.20$\pm$2.65 | 45.74$\pm$1.54 | 49.68$\pm$1.57 | 51.02$\pm$2.11 |  50.48$\pm$1.89
NMI | 10 | 79.67$\pm$2.88 | 78.58$\pm$3.59 | 83.55$\pm$1.62 | 84.07$\pm$1.74 | 74.78$\pm$1.69 | 83.48$\pm$1.89 |  84.20$\pm$2.37 | 85.04$\pm$2.45
20 | 70.64$\pm$1.34 | 69.52$\pm$0.94 | 73.71$\pm$1.00 | 71.66$\pm$1.78 | 67.95$\pm$0.93 | 74.85$\pm$2.77 | 74.00$\pm$1.88 |  74.39$\pm$0.53
30 | 68.53$\pm$2.26 | 67.18$\pm$0.79 | 70.48$\pm$0.78 | 71.12$\pm$0.59 | 68.15$\pm$0.82 | 70.25$\pm$0.98 |  71.35$\pm$0.60 | 71.49$\pm$1.85
40 | 67.39$\pm$1.46 | 65.60$\pm$0.98 | 68.57$\pm$0.51 | 68.65$\pm$0.67 | 66.64$\pm$1.14 | 68.56$\pm$1.10 |  69.39$\pm$0.82 | 69.45$\pm$0.89
50 | 67.53$\pm$0.76 | 66.20$\pm$0.57 | 69.06$\pm$0.32 | 67.88$\pm$1.37 | 67.24$\pm$1.31 | 69.83$\pm$0.80 |  70.26$\pm$0.95 | 70.65$\pm$0.66
PUR | 10 | 78.60$\pm$4.90 | 79.40$\pm$3.87 | 84.90$\pm$1.51 | 84.10$\pm$4.15 | 74.70$\pm$3.03 | 85.20$\pm$2.41 |  85.70$\pm$2.07 | 86.20$\pm$2.51
20 | 61.90$\pm$1.82 | 60.80$\pm$1.67 |  67.15$\pm$1.65 | 63.10$\pm$2.58 | 59.35$\pm$1.64 | 67.10$\pm$4.10 | 66.80$\pm$2.95 | 67.75$\pm$0.50
30 | 54.13$\pm$3.00 | 52.53$\pm$1.32 | 57.70$\pm$1.50 | 59.23$\pm$1.53 | 54.53$\pm$1.74 | 57.53$\pm$1.27 |  59.17$\pm$1.10 | 58.93$\pm$2.03
40 | 50.92$\pm$2.41 | 49.00$\pm$2.07 | 53.27$\pm$1.44 | 52.47$\pm$1.54 | 49.70$\pm$1.23 | 53.02$\pm$1.30 |  53.75$\pm$2.16 | 54.78$\pm$1.60
50 | 48.50$\pm$1.57 | 47.60$\pm$0.83 | 52.32$\pm$0.61 | 50.22$\pm$2.23 | 48.86$\pm$1.68 | 52.66$\pm$1.32 | 53.86$\pm$1.54 |  53.44$\pm$1.64
Table 5: Performance metrics comparisons on Face95 database between K-means, NMF, GNMF, HNMF, NTR-MU, GNTR-MU, HGNTR-MU and LraHGNTR-MU methods. We divide the database into different subset. Each subset contains $k_{Face95}$ categories.  | $k_{Face96}$ | K-means | NMF | GNMF | HNMF | NTR-MU | GNTR-MU | HGNTR-MU | lraHGNTR-MU
---|---|---|---|---|---|---|---|---|---
ACC | 20 | 76.80$\pm$6.96 | 77.00$\pm$3.82 | 77.65$\pm$5.70 | 79.05$\pm$2.04 | 66.35$\pm$3.89 | 75.90$\pm$4.13 |  80.05$\pm$4.73 | 81.00$\pm$6.14
30 | 74.40$\pm$2.48 |  76.83$\pm$2.69 | 72.76$\pm$3.87 | 72.70$\pm$3.93 | 71.63$\pm$2.45 | 71.66$\pm$2.72 | 76.67$\pm$4.79 | 78.87$\pm$3.87
40 | 69.85$\pm$4.97 | 74.60$\pm$2.35 | 66.97$\pm$3.66 | 69.17$\pm$4.10 | 66.27$\pm$2.21 | 71.57$\pm$2.71 | 72.25$\pm$2.24 |  72.47$\pm$4.03
50 | 68.46$\pm$1.46 | 68.76$\pm$1.33 | 69.02$\pm$3.10 | 71.38$\pm$2.49 | 66.64$\pm$1.98 | 71.32$\pm$4.15 |  74.20$\pm$3.00 | 74.32$\pm$2.06
60 | 67.50$\pm$2.42 | 70.06$\pm$2.06 | 70.78$\pm$2.78 | 69.56$\pm$1.45 | 64.76$\pm$3.37 | 71.95$\pm$3.03 |  72.43$\pm$4.10 | 72.50$\pm$2.35
NMI | 20 | 85.55$\pm$2.77 | 85.85$\pm$1.56 | 89.25$\pm$1.97 | 89.96$\pm$0.70 | 78.44$\pm$2.90 | 88.37$\pm$1.74 | 89.98$\pm$2.52 |  89.40$\pm$2.94
30 | 86.75$\pm$1.73 | 88.23$\pm$1.11 | 88.81$\pm$1.47 | 89.58$\pm$1.80 | 85.57$\pm$1.88 | 89.87$\pm$1.06 | 91.20$\pm$1.10 |  91.08$\pm$1.31
40 | 86.77$\pm$1.40 | 88.02$\pm$0.56 | 87.52$\pm$0.52 | 88.05$\pm$1.50 | 84.37$\pm$0.95 | 89.37$\pm$0.57 | 89.57$\pm$0.81 |  89.55$\pm$0.85
50 | 86.52$\pm$1.01 | 87.20$\pm$0.40 | 88.04$\pm$0.95 | 88.10$\pm$1.08 | 85.08$\pm$0.72 | 89.62=1.08 | 90.37$\pm$0.85 |  89.96$\pm$0.53
60 | 86.88$\pm$1.03 | 88.57$\pm$0.47 | 88.55$\pm$0.81 | 88.29$\pm$0.70 | 85.10$\pm$1.87 | 90.29$\pm$0.83 |  90.86$\pm$0.52 | 91.02$\pm$0.43
PUR | 20 | 79.60$\pm$5.19 | 80.80$\pm$2.60 | 82.50$\pm$4.34 | 83.40$\pm$1.77 | 70.30$\pm$4.00 | 81.20$\pm$2.55 | 83.10$\pm$4.28 | 84.10$\pm$4.97
30 | 77.73$\pm$2.03 | 80.36$\pm$1.76 | 78.36$\pm$2.83 | 78.70$\pm$3.18 | 75.10$\pm$1.73 | 77.83=1.92 | 81.30$\pm$3.47 | 82.93$\pm$2.36
40 | 74.55$\pm$3.39 | 78.57$\pm$1.38 | 73.32$\pm$2.32 | 75.32$\pm$3.02 | 70.35$\pm$1.73 | 77.20$\pm$1.79 |  77.50$\pm$2.03 | 76.75$\pm$2.60
50 | 72.70$\pm$1.23 | 73.48$\pm$1.01 | 74.26$\pm$2.15 | 76.04$\pm$1.81 | 69.98$\pm$1.94 | 76.78$\pm$3.18 | 78.48$\pm$2.41 |  78.34$\pm$1.54
60 | 72.15$\pm$1.86 | 75.13$\pm$1.35 | 75.30$\pm$1.89 | 74.41$\pm$0.62 | 68.81$\pm$2.96 | 76.65$\pm$2.40 | 77.73$\pm$2.02 |  77.00$\pm$0.82
Table 6: Performance metrics comparisons on Face96 database between K-means,
NMF, GNMF, HNMF, NTR-MU, GNTR-MU, HGNTR-MU and LraHGNTR-MU methods. We divide
the database into different subset. Each subset contains $k_{Face96}$
categories.

(a) (b)

Figure 4: The running time of various algorithms subset in different number of
categories. (a) ORL. (b) GT.

(a) (b)

Figure 5: The running time of various algorithms subset in different number of
categories. (a) ORL. (b) GT.

###  4.5 Clustering Tasks for Real-world Database with Gaussian Noise

In this part, we count the performance of our algorithm on the noisy dataset.
In the case of noiseless, we can easily find that LraHGNTR is similar to HGNTR
in accuracy, but LraHGNTR is faster. We compare various algorithm performances
by adding the Gaussian noise on ORL and GT database. The degree of noise is
measured by the signal-to-noise ratio (SNR). Before the experiment, as a
result of the negative value may appear in datasets after adding very heavy
Gaussian noise, we make the negative value in each database are truncated at
zero to ensure the negativity for each algorithm except LraHGNTR. The reason
is that LraHGNTR can directly deal with non-negative elements. The
$Fig.\ref{noise_1}$, $Fig.\ref{noise_2}$, $Fig.\ref{noise_3}$ and
$Fig.\ref{noise_4}$ show the comparsion with noise. Through experiment
results, we can get the following conclusions:

  * -

For the datasets with adding Gaussian noise, our algorithms can still achieve
good results compared with other algorithms, and the running time of LraHGNTR
is also reduced. It proves the effectiveness of our algorithms again.

  * -

In the case of mild noise, our approximate algorithm LraHGNTR can achieve
almost the same performance as the HGNTR. And in the case of high noise, such
as SNR = $10\,dB$, the LraHGNTR is obviously better than HGNTR, which
demonstrates that the low-rank approximation method is beneficial to filter
noise.

  * -

Low-rank approximation methods can eliminate the noise effectively. The reason
is that in most datasets with noise, noise usually does not have the low rank
property, so replacing the original tensor $\mathcal{X}$ with an approximate
tensor $\tilde{\mathcal{X}}$ can ignore the noise and thus avoid the influence
of noise.

(a) (b) (c)

Figure 6: Performance metrics $\left(\%\right)$ comparisons on between
K-means, NMF, GNMF, HNMF, NTR-MU, GNTR-MU, HGNTR-MU and LraHGNTR-MU methods
for noise-added ORL and GT database. (a) ACC on ORL. (b) ACC on GT. (c) NMI on
ORL.

(a) (b) (c)

Figure 7: Performance metrics $\left(\%\right)$ comparisons on between
K-means, NMF, GNMF, HNMF, NTR-MU, GNTR-MU, HGNTR-MU and LraHGNTR-MU methods
for noisy-added ORL and GT database. (a) ACC on GT. (b) ACC on ORL. (c) PUR on
GT.

(a) (b) (c)

Figure 8: Performance metrics $\left(\%\right)$ comparisons on between
K-means, NMF, GNMF, HNMF, NTR-MU, GNTR-MU, HGNTR-MU and LraHGNTR-MU methods
for noisy-added ORL and GT database. (a) NMI on Face95. (b) PUR on Face96. (c)
PUR on Face95.

(a) (b) (c)

Figure 9: Performance metrics $\left(\%\right)$ comparisons on between
K-means, NMF, GNMF, HNMF, NTR-MU, GNTR-MU, HGNTR-MU and LraHGNTR-MU methods
for noisy-added ORL and GT database. (a) NMI on Face96. (b) PUR on Face95. (c)
PUR on Face96.

(a) (b)

Figure 10: Running time comparisons on between K-means, NMF, GNMF, HNMF, NTR-
MU, GNTR-MU, HGNTR-MU and LraHGNTR-MU methods for noise-added ORL and GT
database. (a) ORL. (b) GT.

(a) (b)

Figure 11: Running time comparisons on between K-means, NMF, GNMF, HNMF, NTR-
MU, GNTR-MU, HGNTR-MU and LraHGNTR-MU methods for noise-added Face95 and
Face96 database. (a) Face95. (b) Face96.

###  4.6 Parameter Sensitivity

In order to investigate the effect for our model on the situation of parameter
changing, we need to carry out parameter sensitivity test. Considering that
our algorithms has three parameters that can be adjusted, which are the number
of inner iterations $t$, the number of nearest neighbours $k$ and the graph
regularization parameter $\beta$, we use the method of controlling variables.
We divide the experiment into three cases: $\left(1\right)$ fix $k=5$ and
$\beta=0.1$, and choose $t\in\\{20,40,60,80\\}$; $\left(2\right)$ fix $t=20$
and $\beta=0.1$, and choose $k$ to vary from $3$ to $6$; $\left(3\right)$ fix
$t=20$ and $k=5$, and choose $\beta\in\\{0.1,0.2,0.3,0.4,0.5,0.6\\}$.

$Fig.\ref{beta change sensitivity orl}$ and $Fig.\ref{t change sensitivity
orl}$ present the clustering performance of the HGNTR and LraHGNTR across
different parameters on the two public databases. It can be observed that the
performances of HGNTR and LraHGNTR algorithms across different parameters are
quite stable. The ACC, NMI, and PUR change little when the number of inner
iterations $t$ rises from $20$ to $80$, so $t$ can be selected around $20$ to
balance the calculation cost and fitting error. It can be observed that the
ACC, NMI, and PUR are stable when $k$ in the range of $3-6$. It can be noted
that the hypergraph regularization parameter $\beta$ also has a relatively
weak effect on the performance of our both algorithms.

(a) (b) (c)

Figure 12: The algorithm performance curve on ORL database in different
parameters. (a) The graph regularization parameter $\beta$. (b) The number of
inner iterations $t$. (c) The the number of nearest neighbours $k$.

(a) (b) (c)

Figure 13: The algorithm performance curve on GT database in different
parameters. (a) The graph regularization parameter $\beta$. (b) The number of
inner iterations $t$. (c) The number of nearest neighbours $k$.

###  4.7 Convergence Analysis

We construct a tensor $\mathcal{X}\in\mathbb{R}^{40\times 40\times 40\times
40}$ by random selection obeying uniform distribution between 0 and 1, and add
independent Gaussian noisy with SNR = $10\,dB$ and SNR = $15\,dB$ on this
tensor. Then we decomposition it by HGNTR and LraHGNTR. We count the
convergence evolution of the error respectively in $Fig.\ref{convergence
curve}$.

Our algorithms can reach a stable value after one hundred iterations. It can
be proved that both algorithms are convergent. Particularly, we can see the
LraHGNTR has a faster convergence rate in the first half part of the iteration
process. It can verify that the convergence of the low-rank approximate method
is faster. And our low-rank approximation method always reaches the stable
convergence value firstly. Regretfully, in the condition of noisy-added tensor
data, when the convergence is about to complete, our convergence accuracy of
LraHGNTR is a little weaker than the HGNTR, there is a slight difference, this
may is related to the heavy Gaussian noise added on the certain synthetic
manual tensor data, it has completely random properties in the statistic. And
$fig.\ref{convergence_relative}$ shows the relative performance in terms of
running time and final fitting error. A similar conclusion has been shown in
[16].

(a) (b) (c)

Figure 14: The convergence curve on tensor $\mathcal{X}$ and noisy-added
$\mathcal{X}$.The horizontal axis is the number of iterations. The vertical
axis is the fitting error. (a) Original tensor. (b) Noisy-added tensor with
SNR = $10\,dB$. (c) Noisy-added tensor with SNR = $15\,dB$.  Figure 15: The
comparsion of performance in terms of runtime and the final fitting accuracy
on different algorithms. $f^{*}$ and $f^{*}_{max}$ denote the noiseless
fitting error and max noiseless fitting error. $t$ and $t_{max}$ denote the
each running time and max running time among them. $f$ and $f_{max}$ denote
the noisy-added fitting error and max noise-added fitting error.

###  4.8 Basis Visualization

In order to compare the ability that extracts the parts-based features from
database objects by each algorithm, we visual the basis matrices that are
extracted from the ORL and GT database in $Fig.\ref{fig:visual_orl}$ and
$Fig.\ref{fig:visual_gt}$. For the high-dimensional databases, such as ORL and
GT, the matrix methods show very limited capabilities, so that they produce
ghost faces, while the tensor methods can produce the more sparse
representation. The nonnegativity of the tensor decomposition produces about
two main advantages: (1) it can represent the certain potential structure of
the data and be physically interpretable. (2) it can reduce the dimensionality
of the data and memory storage costs. By comparing all methods, our algorithms
HGNTR and LraHGNTR have a higher degree of sparsity. And at the same time,
they retain more physically meaningful features.

(a) (b) (c) (d)

(e) (f) (g)

Figure 16: Visual images of base-matrix based on different algorithms for ORL
database. (a) NMF. (b) GNMF. (c) HNMF. (d) NTR-MU. (e) GNTR-MU. (f) HGNTR-MU.
(g) LraHGNTR-MU.

.

(a) (b) (c) (d)

(e) (f) (g)

Figure 17: Visual images of base-matrix based on different algorithms for GT
database. (a) NMF. (b) GNMF. (c) HNMF. (d) NTR-MU. (e) GNTR-MU. (f) HGNTR-MU.
(g) LraHGNTR-MU.

.

##  5 Conclusion

In this paper, we propose the method of HGNTR. Each dimension information of
the third-order tensor is expressed as the corresponding three-order core
tensor. At the same time, compared with ordinary manifold graph learning,
hypergraph manifold learning can explore more high-order manifold structure
information. An efficient optimization method based on multiplicative updating
rules (MUR) is developed. Considering the computational complexity of the
algorithm, we use the low-rank approximation method to reduce the running
time. The proposed HGNTR and LraHGNTR can achieve better performance than
state-of-the-art (SOTA) algorithms.

##  6 Acknowledgement

The authors would like to thank…

## References

  * [1] D. D. Lee, H. S. Seung, Learning the parts of objects by non-negative matrix factorization, Nature 401 (6755) (1999) 788–791.
  * [2] A. Cichocki, D. Mandic, L. De Lathauwer, G. Zhou, Q. Zhao, C. Caiafa, H. A. Phan, Tensor decompositions for signal processing applications: From two-way to multiway component analysis, IEEE signal processing magazine 32 (2) (2015) 145–163.
  * [3] Q. Zhao, M. Sugiyama, L. Yuan, A. Cichocki, Learning efficient tensor representations with ring-structured networks, in: ICASSP 2019-2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), IEEE, 2019, pp. 8608–8612.
  * [4] L. Yuan, C. Li, D. Mandic, J. Cao, Q. Zhao, Tensor ring decomposition with rank minimization on latent space: An efficient approach for tensor completion, in: Proceedings of the AAAI Conference on Artificial Intelligence, Vol. 33, 2019, pp. 9151–9158.
  * [5] Y. Yu, G. Zhou, N. Zheng, S. Xie, Q. Zhao, Graph regularized nonnegative tensor ring decomposition for multiway representation learning, arXiv preprint arXiv:2010.05657 (2020).
  * [6] M. Belkin, P. Niyogi, V. Sindhwani, Manifold regularization: A geometric framework for learning from labeled and unlabeled examples., Journal of machine learning research 7 (11) (2006).
  * [7] M. Belkin, P. Niyogi, Laplacian eigenmaps and spectral techniques for embedding and clustering., in: Nips, Vol. 14, 2001, pp. 585–591.
  * [8] D. Cai, X. He, J. Han, T. S. Huang, Graph regularized nonnegative matrix factorization for data representation, IEEE transactions on pattern analysis and machine intelligence 33 (8) (2010) 1548–1560.
  * [9] J. Yu, D. Tao, M. Wang, Adaptive hypergraph learning and its application in image classification, IEEE Transactions on Image Processing 21 (7) (2012) 3262–3272.
  * [10] A. Bretto, Hypergraph theory, An introduction. Mathematical Engineering. Cham: Springer (2013).
  * [11] D. Zhou, J. Huang, B. Schölkopf, Learning with hypergraphs: Clustering, classification, and embedding, Advances in neural information processing systems 19 (2006) 1601–1608.
  * [12] K. Zeng, J. Yu, C. Li, J. You, T. Jin, Image clustering by hyper-graph regularized non-negative matrix factorization, Neurocomputing 138 (2014) 209–217.
  * [13] W. Yin, Z. Ma, Q. Liu, Hyperntf: A hypergraph regularized nonnegative tensor factorization for dimensionality reduction, arXiv preprint arXiv:2101.06827 (2021).
  * [14] J. B. Tenenbaum, V. De Silva, J. C. Langford, A global geometric framework for nonlinear dimensionality reduction, science 290 (5500) (2000) 2319–2323.
  * [15] P.-G. Martinsson, V. Rokhlin, M. Tygert, A randomized algorithm for the decomposition of matrices, Applied and Computational Harmonic Analysis 30 (1) (2011) 47–68.
  * [16] G. Zhou, A. Cichocki, S. Xie, Fast nonnegative matrix/tensor factorization based on low-rank approximation, IEEE Transactions on Signal Processing 60 (6) (2012) 2928–2940.
  * [17] G. Zhou, A. Cichocki, Q. Zhao, S. Xie, Efficient nonnegative tucker decompositions: Algorithms and uniqueness, IEEE Transactions on Image Processing 24 (12) (2015) 4990–5003.
  * [18] T. G. Kolda, B. W. Bader, Tensor decompositions and applications, SIAM review 51 (3) (2009) 455–500.
  * [19] Q. Zhao, G. Zhou, S. Xie, L. Zhang, A. Cichocki, Tensor ring decomposition, arXiv preprint arXiv:1606.05535 (2016).
  * [20] Y. Chen, W. He, N. Yokoya, T.-Z. Huang, X.-L. Zhao, Nonlocal tensor-ring decomposition for hyperspectral image denoising, IEEE Transactions on Geoscience and Remote Sensing 58 (2) (2019) 1348–1362.
  * [21] Y.-D. Kim, S. Choi, Nonnegative tucker decomposition, in: 2007 IEEE Conference on Computer Vision and Pattern Recognition, IEEE, 2007, pp. 1–8.
  * [22] T. Hazan, S. Polak, A. Shashua, Sparse image coding using a 3d non-negative tensor factorization, in: Tenth IEEE International Conference on Computer Vision (ICCV’05) Volume 1, Vol. 1, IEEE, 2005, pp. 50–57.
  * [23] W. Wang, V. Aggarwal, S. Aeron, Efficient low rank tensor ring completion, in: Proceedings of the IEEE International Conference on Computer Vision, 2017, pp. 5697–5705.
  * [24] L. Yuan, C. Li, J. Cao, Q. Zhao, Randomized tensor ring decomposition and its application to large-scale data reconstruction, in: ICASSP 2019-2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), IEEE, 2019, pp. 2127–2131.
  * [25] J. A. Bengua, P. N. Ho, H. D. Tuan, M. N. Do, Matrix product state for higher-order tensor compression and classification, IEEE Transactions on Signal Processing 65 (15) (2017) 4019–4030.
  * [26] Y. Pan, J. Xu, M. Wang, J. Ye, F. Wang, K. Bai, Z. Xu, Compressing recurrent neural networks with tensor ring for action recognition, in: Proceedings of the AAAI Conference on Artificial Intelligence, Vol. 33, 2019, pp. 4683–4690.
  * [27] I. V. Oseledets, Tensor-train decomposition, SIAM Journal on Scientific Computing 33 (5) (2011) 2295–2317.
  * [28] C. Hong, J. Yu, J. Li, X. Chen, Multi-view hypergraph learning by patch alignment framework, Neurocomputing 118 (2013) 79–86.
  * [29] C. Wang, J. Yu, D. Tao, High-level attributes modeling for indoor scenes classification, Neurocomputing 121 (2013) 337–343.
  * [30] S. Agarwal, K. Branson, S. Belongie, Higher order learning with graphs, in: Proceedings of the 23rd international conference on Machine learning, 2006, pp. 17–24.
  * [31] Y. Huang, Q. Liu, F. Lv, Y. Gong, D. N. Metaxas, Unsupervised image categorization by hypergraph partition, IEEE Transactions on Pattern Analysis and Machine Intelligence 33 (6) (2011) 1266–1273.
  * [32] Z. Tian, T. Hwang, R. Kuang, A hypergraph-based learning algorithm for classifying gene expression and arraycgh data with prior knowledge, Bioinformatics 25 (21) (2009) 2831–2838.
  * [33] Y. Huang, Q. Liu, S. Zhang, D. N. Metaxas, Image retrieval via probabilistic hypergraph ranking, in: 2010 IEEE computer society conference on computer vision and pattern recognition, IEEE, 2010, pp. 3376–3383.
  * [34] S. Huang, H. Wang, Y. Ge, L. Huangfu, X. Zhang, D. Yang, Improved hypergraph regularized nonnegative matrix factorization with sparse representation, Pattern Recognition Letters 102 (2018) 8–14.
  * [35] P. Drineas, R. Kannan, M. W. Mahoney, Fast monte carlo algorithms for matrices ii: Computing a low-rank approximation to a matrix, SIAM Journal on computing 36 (1) (2006) 158–183.
  * [36] N. Halko, P.-G. Martinsson, J. A. Tropp, Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions, SIAM review 53 (2) (2011) 217–288.
  * [37] M. W. Mahoney, Randomized algorithms for matrices and data, arXiv preprint arXiv:1104.5557 (2011).
  * [38] C. F. Caiafa, A. Cichocki, Generalizing the column–row matrix decomposition to multi-way arrays, Linear Algebra and its Applications 433 (3) (2010) 557–573.
  * [39] N. Vannieuwenhoven, R. Vandebril, K. Meerbergen, On the truncated multilinear singular value decomposition, Numerical Analysis and Applied Mathematics Section (2011).
  * [40] A. Huang, et al., Similarity measures for text document clustering, in: Proceedings of the sixth new zealand computer science research student conference (NZCSRSC2008), Christchurch, New Zealand, Vol. 4, 2008, pp. 9–56.


printacmref=false

The authors have contributed equally to this work. In particular, the UWB
processing was provided by CC and MS, and the stereo-camera processing was
provided by FD, CH. The authors collected the data together at University of
Toronto. [1] [1] [1]

# STAR-loc: Dataset for STereo And Range-based localization

Frederike Dümbgen  University of Toronto Toronto  ,  Mohammed A. Shalaby
McGill University Montreal  ,  Connor Holmes  University of Toronto Toronto  ,
Charles C. Cossette  McGill University Montreal  ,  James R. Forbes  McGill
University Montreal  ,  Jerome Le Ny  Polytechnique Montreal Montreal  and
Timothy D. Barfoot  University of Toronto Toronto

###### Abstract.

This document contains a detailed description of the STAR-loc dataset. For a
quick starting guide please refer to the associated Github repository
(https://github.com/utiasASRL/starloc). The dataset consists of stereo camera
data (rectified/raw images and inertial measurement unit measurements) and
ultra-wideband (UWB) data (range measurements) collected on a sensor rig in a
Vicon motion capture arena. The UWB anchors and visual landmarks (Apriltags)
are of known position, so the dataset can be used for both localization and
Simultaneous Localization and Mapping (SLAM).

##  1\. Summary

This dataset contains multiple trajectories of a custom sensor rig inside a
motion capture (mocap) arena. Attached to the sensor rig are a stereo camera,
providing image streams and inertial measurement unit (IMU) data, and one or
two ultra-wideband (UWB) tags, providing distance measurements to up to 8
fixed and known UWB anchors. Also fixed in the mocap arena are 55 Apriltag
apriltag  landmarks of known position, which are detected online by the stereo
camera. A depiction of the experimental setup can be found in Figure 1.

A compact version of the dataset can be found on Github and is the recommended
starting point for using the dataset. It contains all pre-processed csv files
of the data, as well as some convenience functions for reading the data. If
required, the full dataset can be found on Google Drive. In addition to the
aforementioned csv files, the full dataset also includes the raw bag files and
many analysis plots, making it easier to choose which part of the dataset to
use for the given application. Both the Github repository and drive are
structured as follows:

Figure 1. The STAR-loc dataset contains data collected on different
trajectories of a sensor rig (bottom right) in a Vicon mocap arena, providing
range and calibrated stereo measurements.

  * •

data/<dataset-name>/: folder containing all data of the run of name <dataset-
name> (see Section 3 for an overview of the different runs).

    * –

csv files of processed data (apriltag.csv, imu.csv, uwb.csv, gt.csv. See
data/README.md for information on the data fields.

    * –

calib.json: file with calibration information of this run.

    * –

(Drive only) bag files of this run (<dataset-name>_0.db3 and metadata.yaml)

    * –

(Drive only) some analysis plots of the data, see Appendix A for information
on the plots.

    * –

(Drive only) the video of the left camera with overlaid Apriltag detections.

  * •

mocap/: folder containg environment information.

    * –

uwb_markers_<version>.csv,mocap/environment_<version>.csv: csv files of the
UWB anchor and Apriltag landmark locations, respectively, used in Setup of
<version>. See Section 3 for descriptions of the different setup versions.

    * –

(Drive only) plots of the different setup versions.

    * –

(Drive only) vsk and tag config files used in calibration procedure.

##  2\. System description

A photo of the sensor rig used for data collection is shown in Figure 2.
Attached to the rig are the UWB tag(s) and the stereo camera. A laptop
111Lenovo P16 Thinkpad with 16GB RAM and RTXA4500 NVIDIA GPU (not shown) is
used for driving all sensors and collecting the data through the help of
dedicated Robot Operating System (ROS) nodes.

Figure 2 also depicts the different sensor frames including their transforms.
The frames are:

  * •

$\mathcal{F}_{r}$: robot frame (tracked by mocap system)

  * •

$\mathcal{F}_{c}$: camera frame

  * •

$\mathcal{F}_{t_{i}}$: $i$-th tag frame, with $i=1\ldots 3$

Figure 2. Left: the used sensor rig with one UWB tag and stereo camera
attached. Right: the corresponding CAD model including the definitions of the
sensor frames, showing the three possible UWB tag locations
$\mathcal{F}_{t_{1}}$, $\mathcal{F}_{t_{2}}$ and $\mathcal{F}_{t_{3}}$ and the
stereo camera frame $\mathcal{F}_{c}$.

###  2.1. UWB tag and anchors

The main UWB system used are custom-made boards provided by the Mobile
Robotics and Autonomous Systems Laboratory (MRASL), where each board is fitted
with DWM1000 UWB transceivers222Available at
https://www.qorvo.com/products/p/DWM1000.. The board is shown in Figure 6.
These boards are equipped with a STM32F405RG microcontroller333Available at
https://www.st.com/en/microcontrollers-microprocessors/stm32f405rg.html. to
interface with the DWM1000 modules, and the boards are then connected to the
onboard computer using USB. Details about the ranging and communication
protocols can be found in uwb-hardware . We also provide a small dataset with
the off-the-shelf development kit MDEK1001444Available at
https://www.qorvo.com/products/p/MDEK1001. (setup s5). We use custom data
collection scripts to record distance measurements.

Figure 3. Left: One of eight fixed UWB anchors. Right: subset of Apriltag
landmarks.

###  2.2. Visual landmarks

We use Apriltags apriltag , tag version Tag36h11, as known landmarks. The tags
are printed on US letter paper, glued to foam board, and spread throughout the
room at fixed locations of different heights and orientations, as shown in
Figure 6. We use the open-source libraries apriltag_ros and apriltag_msg for
Apriltag detection.

###  2.3. Stereo camera

We use the ZED2i stereo camera from stereolabs with a
$12\text{\,}\mathrm{c}\mathrm{m}$ baseline and set to $1920\times 1080$
resolution. For data collection, we employ the ROS-based camera driver
provided by stereolabs, using the image_transport_plugins ROS packages for
image compression.

###  2.4. Motion capture

We use a VICON mocap system for ground truth monitoring. For data collection,
we use the ros2-vicon-receiver package by OPT4SMART to publish the native data
stream from VICON as ROS messages.

##  3\. Dataset overview

Table 1. Overview of datasets setup | s1 | s2 | s3 | s4 | s5 |
---|---|---|---|---|---|---
landmarks | v1 | v1 | v2 | v2 | v3 |
trajectory |  |  |  |  |  |  description
loop-2d | x | x | x | x | x |  multiple loops, keeping sensor rig level
loop-2d-fast | x | x | x |  |  |  multiple loops, walking fast, keeping sensor rig level
loop-3d | x | x | x |  | x |  multiple loops, pointing sensor rig at different directions
zigzag |  | x | x | x |  |  a roughly piecewise linear trajectory, keeping sensor rig level
eight |  | x | x |  |  |  draw large eights in the air, pointing at an area dense in Apriltags
apriltag |  |  | x |  |  |  one loop with complete Apriltag detections
ell |  |  | x |  |  |  multiple ell-shaped trajectories with a lot of height variation
grid |  |  | x |  |  |  static measurements at 8 grid points, including one full circle at last point.
loop-3d-z |  |  | x |  |  |  loop with a lot of height variation

###  3.1. Runs and setups

An overview of all collected datasets is shown in Table 1. The datasets are
distinguished by their setup, landmarks and trajectory version. These aspects
are described in more details in the paragraphs below.

#### Trajectory versions

loop-2d

loop-2d-fast

loop-3d

zigzag

eight

Figure 4. Overview of the different standard trajectory types recorded with
many different setups. The shown data corresponds to setup s3.

apriltag

ell

grid

loop-3d

Figure 5. Overview of the trajectory types recorded with the setup s3.

The first five listed trajectory types (also denoted “standard” in what
follows) are repeated for different kinds of tag setups to facilitate their
comparison. The four remaining datasets are recorded for setup s3 only, and
serve as datasets to highlight specific phenomena. Example trajectories for
the 5 standard dataset types are shown in Figure 4.

An overview of the remaining datasets is shown in Figure 5.

#### Setup versions

The datasets include 5 different setups of the sensor rig:

  * •

s1: We use locations 1 and 2 (see Figure 2) for the UWB tags 1 and 2,
respectively, on the sensor rig. The rig is held just above the person’s
height using a single PVC tube.

  * •

s2: hand-held up, 2 tags: The sensor rig is mounted on a second PVC tube to
increase the proportion of line-of-sight (LOS) regions: the rig is now a half
metre above the person’s height.

  * •

s3: hand-held up, 1 tag: We use only location 3 (see Figure 2) for the UWB tag
1.

  * •

s4: mounted on Jackal robot, 1 tag at location 3.

  * •

s5: hand-held up, 1 Decawave tag at location 3.

#### Landmarks versions

The datasets include three different landmark layouts. The first two differ
only minimally in geometry; Version v2 is an improved setup of v1 after
analyzing noise patterns of both the mocap system and the UWB anchors. Version
v3 is the layout used for Decawave data and consists of only four anchors. The
different layouts are shown Figure 6.

v1

v2

v3

Figure 6. The three different UWB anchor layouts and two different Apriltag
landmark layouts used.

###  3.2. Collected data

In this Section, we report some statistics taken over all the collected
datasets. An individual analysis per dataset can be found in Appendix A. The
names of the datasets are of the format:

$\texttt{trajectory-version\\_setup-version}.$

#### Number of messages and rates

The collected data rates and number of messages received, for each dataset,
are shown in Table 2.

Table 2. Overview of all collected datasets | median rate [Hz] | duration [s]
---|---|---
data type | Apriltag | UWB | Vicon | IMU |
loop-2d_s1 | 14.9 | 266.1 | 48.0 | 193.8 | 65.6
loop-2d-fast_s1 | 14.9 | 260.7 | 48.0 | 194.1 | 34.2
loop-3d_s1 | 14.9 | 249.3 | 48.0 | 193.6 | 87.7
loop-2d_s2 | 14.9 | 251.3 | 48.0 | 193.3 | 86.9
loop-2d-fast_s2 | 14.9 | 246.4 | 48.1 | 192.8 | 35.3
loop-3d_s2 | 14.9 | 244.6 | 48.0 | 192.6 | 80.3
zigzag_s2 | 14.9 | 256.4 | 48.0 | 193.0 | 78.0
eight_s2 | 15.0 | 246.6 | 48.1 | 192.3 | 49.2
loop-2d_s3 | 14.9 | 248.6 | 48.0 | 192.9 | 105.3
loop-2d-fast_s3 | 14.9 | 247.8 | 48.0 | 193.1 | 55.5
loop-3d_s3 | 14.9 | 249.4 | 48.0 | 193.9 | 99.8
zigzag_s3 | 14.9 | 249.4 | 48.0 | 193.9 | 201.6
eight_s3 | 14.9 | 247.4 | 48.0 | 192.9 | 98.5
loop-2d_s4 | 11.3 | 247.4 | 48.3 | 196.2 | 194.1
zigzag_s4 | 14.3 | 249.1 | 48.0 | 198.0 | 218.2
grid_s3 | 14.9 | 248.6 | 48.0 | 192.7 | 112.9
ell_s3 | 14.9 | 248.7 | 48.0 | 192.6 | 111.1
loop-3d-z_s3 | 14.9 | 249.1 | 48.0 | 193.2 | 150.3
apriltag_s3 | 14.9 | 247.8 | 48.0 | 193.1 | 116.4
loop-2d_s5 | 14.9 | 40.0 | 48.0 | 193.2 | 82.0
loop-2d-v2_s5 | 14.9 | 40.0 | 48.0 | 193.6 | 101.0
loop-3d_s5 | 14.9 | 40.0 | 48.0 | 193.1 | 67.0
median rates [Hz] | 14.9 | 248.6 | 48.0 | 193.2 |

####  UWB quality

The quality of UWB measurements per dataset, measured by the (signed)
difference between the ground truth and measured distances, is shown in Figure
7. The measurements shown are before calibration. For performance after
calibration, please refer to the individual dataset plots provided in the
dataset on Google Drive.

Figure 7. UWB performance over all datasets.

#### Stereo quality

The quality of stereo measurements per dataset is measured by the difference
between the projected pixel values using the ground truth landmark locations
and camera pose, and the measured pixel values, as shown in Figure 8. The
measurements shown are before calibration.

Figure 8. Stereo camera performance over all datasets.

##  4\. Data processing

###  4.1. Stereo camera calibration

Stereo camera calibration involves two main steps:

  1. (1)

Data association: Involves associating Apriltag identifiers to the ground
truth Vicon locations.

  2. (2)

Extrinsic calibration: Involves solving for the transformation between the
Vicon frame attached to the camera rig and the left Zed2i camera frame (at
focal point).

Note that the rectified intrinsic parameters provided for the Zed2i camera
were assumed to be accurate and were not tuned in this calibration procedure.

Data association was performed using a camera sequence with Apriltaglabels to
manually label the set of Vicon ground truth landmark locations.

For extrinsic calibration, we provide three different methodologies: fixed
calibration (appendix _cal), global calibration (_cal_global) and individual
calibration (_cal_individual) For fixed calibration, the CAD model of the
sensor rig was used to establish a transformation between Vicon and camera
frames. For global calibration, we use the dataset apriltag_s3 to find the
transformation that minimizes the reprojection (also called photometric) error
using Gauss-Newton optimization state-estimation . For individual calibration,
we use the same procedure, but calculate the best transformation for each
dataset individually.

To reduce the effect of outliers, we perform the calibration three times,
removing outliers after each iteration. Outliers are pixel measurements with
errors higher than a fixed threshold (we used 100 for v2 and v3 and 300 for
v1).

###  4.2. UWB calibration

The ranging protocol used in these experiments is the double-sided two way
ranging (DS-TWR) protocol shown in Figure 9. Each of the 6 timestamps
$t_{i},i\in\\{1,\ldots,6\\}$ is recorded as part of the dataset, as well as
the received signal power at the timestamps $t_{2}$ and $t_{4}$.

Figure 9. The DS-TWR protocol used in the experiments to generate range
measurements between any two UWB tags.

A time-of-flight measurement can be generated from the recorded timestamps as

$t_{f}=\frac{1}{2}\left(\Delta t_{41}-\frac{\Delta t_{64}}{\Delta
t_{53}}\Delta t_{32}\right),$

where, as shown in Figure 9, $\Delta t_{ji}=t_{j}-t_{i}$. Nonetheless, the
range measurements are biased as shown in the blue distribution in Figure 11.
This error stems from multiple factors such as timestamping inaccuracies,
multipath propagation, and obstacles.

Figure 10. Histogram representing the distribution of the bias of the range
measurements for the zigzag s2 experiment. The distribution for the raw and
the calibrated measurements are shown.

Figure 11. The relation between the bias and the standard deviation of the
bias for the measurements as a function of the lifted received signal power.
These are obtained from a separate calibration experiment.

Table 3. Calibrated antenna delays for all tags using a separate calibration experiment. Tag ID | Antenna Delay [ns]
---|---
1 | -0.5672
2 | -0.2139
3 | -0.2499
4 | -0.5514
5 | -0.3756
6 | -0.3292
7 | -0.1349
9 | -0.2914
10 | -0.3000
11 | -0.3953
12 | -0.3069

To obtain more accurate range measurements, the calibration procedure
presented in shalaby2023  is implemented. The antenna delay values shown in
Table 3 and the power-to-bias relation shown in Figure 11 are obtained by
performing two separate experiments with 6 UWB tags fitted on 3 flying
quadcopters (2 per quadcopter). A relation between the uncertainty or standard
deviation of measurements as a function of the received signal power is also
learned and is shown in the bottom plot of Figure 11. More information on the
experiment and how these relations are obtained can be found in shalaby2023 .

The timestamps are then corrected using the antenna delays and the received-
signal power and the range measurements are recomputed, giving much less
biased measurements as shown in Figure 11. Both the raw and calibrated
measurements are provided when available, and they can be found under the
columns range and range_calib, respectively. The uncertainty of each
measurement is also provided under the std column.

##  5\. Acknowledgments

We would like to thank Joey Zou for his help in preparing the Decawave data
collection pipeline and rig adjustments. We also thank Connor Jong for his
help during experiments and Abhishek Goudar for advice on UWB data collection.

## References

  * (1) Mohammed Ayman Shalaby, Charles Champagne Cossette, James Richard Forbes, and Jerome Le Ny, _“Calibration and Uncertainty Characterization for Ultra-Wideband Two-Way-Ranging Measurements,”_ arXiv:2210.05888, 2022.
  * (2) Edwin Olson, _“AprilTag: A robust and flexible visual fiducial system,”_ in 2011 IEEE International Conference on Robotics and Automation, 2011, pp. 3400-3407, doi: 10.1109/ICRA.2011.5979561.
  * (3) Mohammed Ayman Shalaby, Charles Champagne Cossette, Jerome Le Ny, and James Richard Forbes._“Multi-Robot Relative Pose Estimation and IMU Preintegration Using Passive UWB Transceivers,”_ arXiv:2304.03837, 2023.
  * (4) Barfoot, Timothy D. State Estimation for Robotics. Cambridge University Press, 2017.

##  Appendix A Individual dataset plots

Each dataset folder in data/<dataset_name> contains a number of plots for fast
data inspection. The different plots are described in more detail below.

Figure 12. Ground truth data, from left to right: projection of trajectory to
x-y plan (gt-xy.png), translation over time (gt-t-xyz.png) and rotation over
time, using axis-angle representation (gt-t-angles).

Figure 13. Apriltag data, from left to right: mask of Apriltag detections over
time (apriltag_mask.png), black corresponding to at least one detection in
either the left or the right stereo camera frame, and number of synchronized
Apriltag detections over time (apriltag_number.png).

Figure 14. Stereo camera pixel errors, using no calibration (top,
apriltag_cal.pdf), one global calibration (middle, apriltag_cal_global.pdf)
and individual calibrations (bottom, apriltag_cal_individual.pdf). We also
provide plots cropped at $\pm$ 300 pixels (_crop).

Figure 15. UWB data, top: UWB bias before and after calibration
(uwb_bias_histogram.png), bottom: raw and calibrated measurements compared to
ground truth (uwb_measurements_tag_1.png).

  *[mocap]: motion capture
  *[IMU]: inertial measurement unit
  *[UWB]: ultra-wideband
  *[ROS]: Robot Operating System
  *[CAD]: Computer Aided Design
  *[LOS]: line-of-sight