t-HGSP: Hypergraph Signal Processing Using t-Product Tensor Decompositions

Graph signal processing (GSP) techniques are powerful tools that model complex relationships within large datasets, being now used in a myriad of applications in different areas including data science, communication networks, epidemiology, and sociology. Simple graphs can only model pairwise relationships among data which prevents their application in modeling networks with higher-order relationships. For this reason, some efforts have been made to generalize well-known graph signal processing techniques to more complex graphs such as hypergraphs, which allow capturing higher-order relationships among data. In this article, we provide a new hypergraph signal processing framework (t-HGSP) based on a novel tensor-tensor product algebra that has emerged as a powerful tool for preserving the intrinsic structures of tensors. The proposed framework allows the generalization of traditional GSP techniques while keeping the dimensionality characteristic of the complex systems represented by hypergraphs. To this end, the core elements of the t-HGSP framework are introduced, including the shifting operators and the hypergraph signal. The hypergraph Fourier space is also defined, followed by the concept of bandlimited signals and sampling. In our experiments, we demonstrate the benefits of our approach in applications such as clustering and denoising.

. In a co-authoship network, two different hypergraphs (b) H 1 (c) H 2 are mapped to the same simple graph in (a) by the clique expansion. Hyperedges are color-coded by publication, e.g. the red hyperedge (e 3 ) indicates that Carl, Dan, and Ed coauthored a publication. Since e 4 ⊂ e 3 , e 4 is ignored in the clique expansion representation. and signal representations have thus been extended for signals on graphs under the umbrella of GSP [5], [6], [7].
At the heels of GSP, a new and exciting set of tools for processing data on complex structures is emerging, coined Hypergraph Signal Processing (HGSP) [8], [9], [10], [11], [12]. The motivation lies first in the ability of hypergraphs to capture high-order interactions between more than two nodes. In the case of a co-authorship network, for instance, a hypergraph can better represent the relationships between the different groups of authors. While edges in conventional graphs ( Fig. 1(a)) can only model pairwise relationships between nodes, limiting GSP to single-way analysis, hyperedges connect more than two nodes as illustrated in Fig. 1(b)-(c). Hypergraph signals are thus defined as those associated to the vertices of a hypergraph, whose polyadic interactions are modeled by hyperedges. The literature on data processing in high-order networks is sparse compared to GSP [12], even though higher-order relationships are widespread in many complex systems such as multiple neurons firing at the same time [13], biochemical reactions having more than two proteins [1], people interacting in small groups of more than two people [14], and multilateral relationships among several related points (e.g., on a surface) in 3D point cloud [15].
A second motivation for advancing hypergraph signal processing is to overcome the limitation of standard graphs to model only single-layer connections. Many complex networks in nature exhibit multiple-layer interactions between two nodes, which can be modeled by multilayer graphs with each type of interaction defining a layer [16], [17], [18]. Mathematically, the matrix algebra of graphs cannot effectively model multi-layer relationships as all entries in a matrix are considered equivalent [9]; however, networks with multi-layer relationships have been successfully studied by being mapped to hypergraphs [19]. Among the applications of multilayer graphs are multi-sensor point clouds, imaging spectroscopy, and multilayer social networks. Point cloud processing, in particular, is an area of significant growth as it arises with new sensor technologies such as LiDAR [20], tomographic synthetic aperture radar [21], and time-of-flight spectral cameras [22], [23], [24]. Point clouds are ubiquitous in computer vision and are essential in the exploration of planetary bodies [25], [26].
The advantages of HGSP, however, are currently hindered by several technical challenges including: 1. Choosing an appropriate algebraic descriptor: In GSP, computation is enabled by encoding the graph structure in an adjacency matrix or its associated Laplacian. Hypergraphs, however, have many different algebraic descriptors for the same higher-order interaction structure, and choosing an appropriate algebraic descriptor is a key element when defining an HGSP framework. In fact, hypergraphs can be described by both matrices [27] and tensors [28]. GSP operations can easily be applied to hypergraph matrix-based descriptors. Nevertheless, these fail to capture the high-order structure of the data. Higher dimensional tensors, on the other hand, allow capturing deeper insights from high-order data.
2. Defining a hypergraph signal shifting operation: In GSP, the algebraic descriptors of the graph are called graph shift operators which are a generalization of the classical time delay in DSP. The shifting operation replaces a signal value at each node of a graph with a linear combination of the signal values at the neighbors of that node. When defining a new tensor-based HGSP framework, modeling the hypergraph signal interactions over hyperedges that connect more than two nodes is not as straightforward as in the graph case.
3. Defining a loss-free hypergraph Fourier transform: In GSP, the graph Fourier transform is defined based on the eigendecomposition of the graph shift operator. In HGSP when considering tensor algebraic descriptors, the hypergraph Fourier transform is not uniquely generalized since the concept of eigendecomposition of a matrix is not uniquely defined for high-order tensors. Previous efforts relied on an approximated, yet computationally intensive decomposition, that leads to subspace mappings of the signal, hence, providing a lossy Fourier transform.
This article tackles these limitations by proposing a new HGSP framework based on t-product factorizations. The main contributions of this article can be enumerated as follows: (i) formulate the core elements of the new HGSP framework, which includes the tensor-based algebraic representation of hypergraphs, hypergraph signals, and the hypergraph shifting operation; (ii) define a loss-free hypergraph Fourier transform using t-product factorizations; (iii) provide an analysis of frequency and the concept of bandlimited signals by defining a measure of the total variation of signals on hypergraphs; (iv) introduce the fundamentals for filtering, sampling, and recovery of hypergraph signals.

A. Graph Signal Processing Overview
An undirected simple graph on N data points is given by G = (V (G), E(G)), where V (G) is the set of nodes defined on the N data points and E(G) is the set of edges describing the pairwise interactions between these nodes. These interactions can also be represented by the symmetric adjacency matrix A ∈ R N ×N whose entries are A(u, v) = 1 if there is an edge connecting the nodes u and v and 0 otherwise. The combinatorial graph Laplacian associated with G and computed as L = D − A is a symmetric positive semidefinite matrix where the diagonal matrix D stores the degrees of the nodes in the graph and its entries are given by D(u, u) = v∈V (G) A(u, v). A real valued signal, x ∈ R N , on the graph, G, is defined as a function x : The component x(v) represents the value of the signal on the node v ∈ V (G). From G = (V (G), E(G)), GSP defines a shift operator, F, as a local operation that replaces the signal's value at each node with a linear combination of the signal's values from neighboring nodes according to Fx. Several graph shift operators such as the adjacency or the Laplacian can be considered in order to exploit different properties of the graph and to define a graph signal processing framework [6], [29].
The notion of linear filtering the graph signal, x, by the filter, Q, is achieved by multiplication y = Qx, resulting in a new graph signal y. While Q can be any arbitrary matrix, it is considered a linear shift invariant (LSI) operator if it satisfies the condition that QFx = FQx. Now since F is real and symmetric, it is diagonalizable by means of the eigenvalue decomposition as F = VΛV T where V is a unitary matrix whose column vectors are eigenvectors and Λ is a diagonal matrix of corresponding eigenvalues [30], [31]. Then, given a filtering weight function q : R → R, a linear shift invariant filter can be written as is the frequency response of the filter Q. By spanning the set R N such that any graph signal, x, can be written as a linear combination of column vectors from V, the Graph Fourier Transform (GFT) of a signal x on G can then be defined asx = V T x and its inverse GFT as x = Vx with the eigenvalues interpreted as frequencies [32]. As N -length vectors, the columns of V are, themselves, graph signals whose corresponding eigenvalues are measures of their total variations. The bandwidth of x is defined in terms of the nonzero components ofx. can be represented using matricial or tensorial algebraic descriptors. In the matrix-based approach, a hypergraph is represented by the incidence matrix B ∈ R N ×E such that B ue = 1 if the vertex u is an element in the hyperedge e and 0 otherwise. Even though the matrix-based representation of the hypergraph is simple, it fails to capture the high-dimensionality structure of the hypergraph. For this reason, in this work, we use the tensorial representation of hypergraphs to propose a new HGSP framework.

B. From Graphs to Hypergraphs
The adjacency tensor of a hypergraph H = (V (H), E(H)) with N nodes and M = m.c.e(H) is traditionally represented by an M th-order N -dimensional tensor A ∈ R N M defined as Observe that 1 enumerates all the possible combinations of c i positive integers {k 1 , k 2 , . . . , k c i } whose summation satisfies Note that a hypergraph with M = 2 degrades to a simple graph with a binary adjacency matrix A ∈ R N ×N ; hence, HGSP is a generalization of GSP.
The degree of a vertex v k ∈ V (H) is the number of hyper-

C. Related Work
The theory of signal processing on higher-order networks is largely unexplored compared to that on graphs. The most prominent abstractions for such polyadic data are simplicial complexes [33] and hypergraphs [11], [34]. Simplicial complexes are used as an alternative representation of hypergraphs, with the particular drawback of the inclusion property, characteristic of simplicial complexes, that is undesirable when representing interactions that are exclusive to multiple nodes but do not imply the interaction between all the subsets of nodes [11], [12]. For instance, in a co-authorship network ( Fig. 1), having a paper of three or more authors does not imply that these authors have written papers in pairs. While tensor-based hypergraph representations can differentiate these two cases, simple graphs, simplicial complexes, and even matrix-based hypergraph representations, which represent hypergraphs as graphs, cannot [12]. As shown in Fig. 1, a common matrix-based hypergraph representations, the clique expansion, replaces every hyperedge with a clique subgraph, clearly failing to capture lower-dimensional relationships and not providing an injective mapping for a hypergraph.
Recently, a tensor-based hypergraph signal processing framework was introduced [9], which proceeds analogously to GSP. In [9], much like the Fourier Transform in GSP is defined via the eigen-decomposition of the graph Laplacian or adjacency matrix, the Fourier transform in HGSP is defined via the canonical polyadic (CP) tensor decomposition of the hypergraph adjacency or Laplacian tensor. The canonical polyadic (CP) decomposition expresses a tensor as the finite sum of rank-one tensors ( Fig. 2(a)). For instance, the CP decomposition of a third-order tensor X ∈ R N 1 ×N 2 ×N 3 is represented , and c r ∈ R N 3 , and • represents the tensor outer product operation. Particularly, the symmetric orthogonal-CP decomposition [35], being a special case of the CP decomposition, is used in [9] to decompose a super-symmetric tensor F ∈ R N ×N ×N representing either the hypergraph adjacency or Laplacian tensor as where the f r 's form an orthogonal set. The integer R is the smallest number of rank-one tensors required to express F and is referred to as the CP rank of F [36]. The CP rank has no relaxation which limits its applications. While the CP-tensor decomposition approach to HGSP was shown effective in some image processing applications [9], [10], it has some drawbacks. The first drawback comes from the fact that the hypergraph spectrum space is obtained from the orthogonal-CP decomposition of the hypergraph adjacency tensor A, and, as shown in Fig. 3(a), this tensor does not have an exact CP orthogonal decomposition. In fact, most tensors generally do not have an exact CP orthogonal decomposition [37]. Current efforts have focused on finding symmetric orthogonal approximations [35]; however, errors add uncertainty, affecting the accuracy of the hypergraph Fourier transform which is undesirable, particularly when filtering on the Fourier domain. Second, finding the CP-rank R is an NP-hard problem which invariably poses high computational complexity [38]. Third, generally, the rank R is smaller than the number of nodes in the hypergraph N [9]. As a consequence, hypergraph signals in R N are invariably mapped to a subspace of R N which does not guarantee perfect recovery on the inverse hypergraph Fourier transform for any type of signal as depicted in Fig. 3(a).
An alternative setting for hypergraph signal processing that can overcome all of the above-mentioned limitations is therefore needed. Choosing the tensor decomposition that better parallels the characteristics of the eigendecomposition in traditional GSP frameworks becomes a key aspect when defining an HGSP framework. To this end, we adopt an alternative setting based on the more recently proposed tensor-tensor multiplication (t-product), in which the familiar tools of linear algebra are extended to better understand tensors [39]. The t-product of two small tensors a ∈ R 1×1×N 3 and b ∈ R 1×1×N 3 , known as tubal scalars, is the tensor c ∈ R 1×1×N 3 computed as c = a b, where denotes the circular convolution of two vectors. Thus, by considering tensors as matrices of tubal scalars, tensor factorizations built from t-products are analog to matrix factorizations  [41] and its (top-left) 3 rd -order adjacency tensor representation A. At the bottom of (a) and (b) is the signal after applying the forward and inverse hypergraph Fourier transform and the corresponding absolute error using: the symmetric orthogonal CP-based decomposition, and the proposed t-product factorization, respectively. At the top of (a) and (b) are the absolute errors of the tensor decompositions. The advantages of the proposed approach are readily seen.
such as the SVD, QR, and eigendecompositions [39]. The tensor singular value decomposition, based on this t-product, was coined t-SVD. The significance of these decompositions is that they allow for the extension of familiar matrix analysis to the multilinear setting while preserving the intrinsic structure of tensors and avoiding the loss of information. The algebraic framework for higher-order tensors is thus analogous to that of matrices such as the transpose, identity, orthogonality, and diagonalization. The t-SVD of a third-order tensor X ∈ R N 1 ×N 2 ×N 3 is for instance X = U * S * V T where * denotes the t-product [39], U ∈ R N 1 ×N 1 ×N 3 and V ∈ R N 2 ×N 2 ×N 3 are orthogonal tensors, and S ∈ R N 1 ×N 2 ×N 3 is an f-diagonal tensor where the frontal slices are diagonal matrices as shown in Fig. 2(b). Although the t-SVD was first demonstrated on 3rd-order tensors, it is easily extended to n-th (n > 3) order tensors [40]. The proposed research exploits this novel set of t-product factorizations to develop a new HGSP framework, dubbed as t-HGSP, that aims at being more stable and, more importantly, loss-free compared to previous tensor-based HGSP frameworks as demonstrated in Fig. 3(b). To formulate the t-product and its algebra as the backbone of our t-HGSP framework the following definitions need to be introduced.

D. t-Product Operations
Denoting vectors by bold lowercase letters (e.g. a), matrices as uppercase letters (e.g. A), and tensors as calligraphic letters (e.g. A), we first define the (i, j)-th tube scalar of the 3rd-order tensor A as a ij which is illustrated in Fig. 4(b). The j-th lateral slice of the tensor, shown in Fig. 4(c), is denoted as A j ≡ A(: , j, :) ∈ R N 1 ×1×N 3 which is a vector of tubal scalars. The k-th frontal slice depicted in Fig. 4(d) is denoted as A (k) ≡ A(:, : , k) ∈ R N 1 ×N 2 ×1 which is a matrix.
Definition 1 (t-product [40], [42]): As mentioned before, the t-product of two tubal scalars a ∈ R 1×1×N 3 and b ∈ R 1×1×N 3 is also a tubal scalar c ∈ R 1×1×N 3 computed as c = a b, where denotes the circular convolution of two vectors. Then, the t-product of two vectors of tubal In general, the t-product of two 3rd-order tensors A ∈ (2) . . .
where the operator bcirc(A) converts the set of frontal slices of the tensor A into a block circulant matrix and unfold(B) stacks vertically the set of frontal slices of B into a N 2 N 3 × N 4 matrix. The operator fold() reverses this process, fold(unfold(A)) = A. Circulant matrices are diagonalized by the discrete Fourier transform; hence, the t-product can be computed efficiently in the Fourier domain as explained in [42]. Using MATLAB notation, letÂ := fft(A, [], 3) denote the tensor obtained by applying the fast Fourier transform (FFT) along each tubal element of A. For the remainder of this article, the hat notation refers to a tensor that has gone through this operation. Thus, the t-product of A ∈ R N 1 ×N 2 ×N 3 and B ∈ R N 2 ×N 4 ×N 3 can also be computed by matrix multiplication of each pair of frontal slices in the Fourier Domain as The t-product can be easily extended to high-order tensors in a recursive manner. The t-product of Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply. (2) . . .
where A (l) and B (l) are order-(p − 1) tensors formed from holding the p-th index at l, respectively for l = 1, 2, . . . , N p [40]. Thus, each successive t-product operation involves tensors of one order less and at the base label of recursion there is a t-product of 3rd-order tensors. Given that the extension to high-order tensors is obtained by recursion, for simplicity and without loss of generality, we present our framework for the base case which is given by hypergraphs with M = 3 whose adjacency tensor is a 3rd-order tensor A ∈ R N ×N ×N .
Definition 2 (Transpose and Symmetric tensors [42]): The transpose of a 3rd-order tensor A ∈ R N 1 ×N 2 ×N 3 , denoted as A T , is the tensor obtained by transposing each of the frontal slices and then reversing the order of the transposed frontal slices 2 through N 3 . For a higher-order tensor . . , N p and then reversing the order of the A (l) s from l = 2 to l = N p as The that is a ≤ b if and only if each of the components of a in the Fourier domainâ (k) are pair-wise less than or equal than the components of b in the Fourier domainb (k) . Definition 4 (Norm [43]): The norm of a vector of tubal scalars X ∈ R N 1 ×1×N 3 is a 1 × 1 × N 3 tubal scalar and it can be computed as 3). Definition 5 (Inverse of a Tensor [42]): where I N 1 N 1 N 3 is the identity tensor whose first frontal slice is the N 1 × N 1 identity matrix, and the other frontal slices are all zeros. Definition 6 (t-eigendecomposition [42]): where Λ ∈ R N 1 ×N 1 ×N 3 is an f-diagonal tensor, and the following holds where λ j ∈ R 1×1×N 3 and V j ∈ R N 1 ×1×N 3 corresponds to an eigen-tuple and its corresponding eigen-matrix.
where V is an orthogonal tensor and the set of eigen matrices { V 1 , V 2 , . . . , V N 1 } form an orthonormal set (Definition [13] in the Supplemental Material). This decomposition can be efficiently computed in the Fourier domain by means of matrix eigendecomposition as shown in Algorithm 1 (for 3rd-order tensors).

III. HYPERGRAPH SHIFTING OPERATION
In order to make these t-product operations the backbone of our t-HGSP framework, we define the hypergraph shifting operation as where X is the hypergraph signal, Y is the one-time shifted/filtered signal and F is the shifting operator. As in GSP, F is any operator that captures the relational dependencies between nodes, including the adjacency tensor A and the Laplacian L. Motivated by the connection between the t-product and the Fourier transform [39] and the fact that, in this article, we are only considering undirected hypergraphs, we also define symmetric adjacency and Laplacian tensor descriptors. Note that, for M > 2, the tensors A and L introduced before are not symmetric according to Definition 2. Therefore, we define a symmetrization operator sym(A) that generates a symmetric version A s ∈ R N ×N ×(2N +1) of A ∈ R N ×N ×N , by adding a matrix of zeros 0 N ×N as the first frontal slice, dividing by 2, and reflecting the frontal slices of A along the third dimension as Now, if we let N s = 2N + 1, for a higher-order tensor obtained by recursively appending a (M − 1)thorder tensor of zeros O at the front, dividing by 2, and reflecting the (M − 1)th-order tensors A (l) along the last dimension as where the (M − 1)-order tensors are obtain from holding the M -th index at l, respectively, for l = 1, 2, . . . , N. When applied to the degree tensor and the Laplacian tensor, we obtain D s and L s , respectively. Notice that this operation is reversible. The symmetric hypergraph shift is then given by In Section IV-A, we include further discussions on the effects of the symmetrization operation given by the connection of the t-product and the Fourier transform. In GSP, the graph signal is defined as an N-length vector x = [x 1 , . . . , x N ] where each signal element is related to one node in the hypergraph. In the proposed framework, given that the shifting operation is defined by the t-product and the shifting operator is a tensor of dimension N × N × N s × · · · × N s , the hypergraph signal X s and its one-time shifted signal Y s should both be tensors of size N × 1 × N s × · · · × N s to have consistent operations. Thus, we now relate a tubal scalar ( be a set of one dimensional signals on the hypergraph. One can consider assigning to each tubal scalar (x s ) i,1 ∈ R 1×1×N s ×···×N s , the set of signals corresponding to that node, i.e. {x i,1 , . . . , x i,L }, 1 ≤ i ≤ N . Then a total number of L ≤ N (M −2) s signals can be jointly processed on the hypergraph. This setting can be particularly useful when processing signals that are correlated along the tubes, as it is the case of hyperspectral information. We demonstrate this in our experiments by considering the sparse representation of a hyperspectral point cloud.
Definition 9 (Hypergraph Signal from One Dimensional Signal): For a hypergraph with N nodes and m.c.e(H) = M , in [9], an alternative form of hypergraph signal is defined as an (M − 1)th-order N -dimensional tensor X computed as the outer product of an original signal in the hypergraph This definition of hypergraph signal can be easily adapted to the proposed t-product framework by expanding on a new second dimension, X = expand(X) where X is an M th-order tensor with dimensions N × 1 × N × · · · × N and by computing its symmetric version as X s = sym( X) such that X s ∈ R N ×1×N s ×···×N s . Notice that as in the CP-based HGSP framework [9], the hypergraph signal X s is just another representation of an original one dimensional signal x that aims at reflecting its properties in different dimensions. For instance, for a hypergraph with M = 3, the hypergraph signal highlights the properties of the 2-D signal components x i x j .
As an example, let the adjacency tensor be the shifting operator, i.e. F s = A s ∈ R N ×N ×N s , and consider the 3-uniform hypergraph in Fig. 5. The shifted signal in node v 7 is then computed as where (a s ) 7,2 and (a s ) 7,3 are tubal scalars of the symmetrized adjacency tensor A s that represent the hyperedge e 2 = {v 2 , v 3 , v 7 } and (a s ) 7,5 and (a s ) 7,6 represent the hyperedge e 3 = {v 5 , v 6 , v 7 }, which are the only two hyperedges that contain the node v 7 . Note that the entries of the adjacency tensor are very sparse, the tubal scalar (a s ) 7,2 , for instance, has only two elements different from zero: (a s ) 7,2 which corresponds to the entry a 723 ∈ A and its reflection (a s ) (13) 7,2 which is product of the symmetrization operation. Given that (a s ) (4) 7,2 = (a s ) (13) 7,2 = a 723 /2, we will use, for simplicity, the entries of the adjacency tensor a i,j,k ∈ A in this example. Then, for any hypergraph signal, the l-th frontal slice of the i-th tubal scalar (y s ) i,1 can be computed as: Now, on one hand, one could consider a hypergraph signal obtained from a set of signals X as in Definition 8, then from 18, the frontal slices of the shifted signal in node v 7 are given by which represents a linear combination of the signal tubes at neighboring nodes. Note that there is a sequential shift along the third dimension for each frontal slice of the shifted signal. This sequential shift along the third dimension is beneficial to signals that are naturally correlated along the third dimension, as in the case of hyperspectral information. Then, the shifting across the third dimension can exploit the correlation among different channels and benefit applications such as compression and denoising. On the other hand, one could consider a hypergraph signal from one dimensional signal x as in Definition 9, then from 18, the frontal slices of the shifted signal in node v 7 are where x i ∈ x. Note that the first frontal slice (y s ) correspond to the shifted signal proposed in the CP-based HGSP Algorithm 1: t-eigendecomposition for a 3-order tensor [42]. approach [9] but divided by two. Different from [9], the frontal slices 2 through N + 1 consider the cross-product of the signals at nodes connected to node 7 and the rest of the nodes in the graph.

IV. HYPERGRAPH SPECTRUM SPACE
Central to HGSP is the Fourier transform. In our formulation, the generalization comes directly from the tensor eigendecomposition (t-eigendecomposition). Given the symmetric shifting operator F s ∈ R N ×N ×N s , the t-eigendecomposition (Algorithm 1) is determined by where V = [ V 1 , . . . , V N ] ∈ R N ×N ×N s is an orthogonal tensor, Λ = diag(λ 1 , . . . , λ N ) ∈ R N ×N ×N s is an f-diagonal tensor whose frontal slices are diagonal matrices, and (λ j ∈ R 1×1×N s , V j ∈ R N 1 ×1×N 3 ) corresponds to the j-th eigen-pair of eigen-tuple and its corresponding eigen-matrix. The hypergraph Fourier transform of a hypergraph signal X s ∈ R N ×1×N s (t-HGFT) is then X F s = V T * X s , with the inverse hypergraph Fourier transform (t-IHGFT) given by X s = V * X F s . Since the tensor V forms an orthonormal set with cardinality N , i.e., V T * V = V * V T = I NNN s , perfect Fourier representation and recovery of the signal is achieved. See Fig. 3(b).

A. Hypergraph Frequency
Given the above definition of the hypergraph Fourier transform, the next step is to determine a notion of frequency, which is key for filtering. In GSP, the notion of frequency is provided by the eigenvalues and eigenvectors of the representing matrix F ordered based on the total variation (TV) [44]. The total variation and hence the frequency describes the rate of signal changes among neighbors. For t-HGSP, the total variation of the hypergraph Fourier basis can also be used to define a notion of frequency. To this end, we first define the total variation on hypergraphs around the concept of the hypergraph shifting and use this definition to determine an ordering of the eigen-tuples of the adjacency tensor, A s , that provides a notion of frequency.
Definition 10 (Total Variation on a Hypergraph): Let F s = A s . The total variation of a hypergraph signal X s on H is defined as the similarity measure between the hypergraph signal X s and its shifted version Y s computed as (k) . The normalization ensures numerical stability when filtering as it prevents excessive scaling of the shifted signal [44]. The norm · is computed according to the t-product Definition 4.
Since the TV describes how oscillatory the signal is on the hypergraph, it leads to the concept of low and high frequencies on the hypergraph Fourier spectrum, V, as shown next.
Theorem 1: Let F s = A s and λ i and λ j be two eigen-tuples of F s with corresponding eigen matrices V i and V j . Given that F s is symmetric, it has real eigen-tuples. If the eigen-tuples are ordered as λ i > λ j , then the total variations of their eigenvectors Proof: See Appendix A. Theorem 1 provides an ordering of frequency according to the eigen-tuples λ j 's. As a result, the hypergraph Fourier basis of the adjacency tensor A s which are the eigen-pairs (λ j , V j ) are ordered from lowest to highest frequency if λ 1 ≥ λ 2 ≥ · · · ≥ λ j ≥ · · · ≥ λ N . Recall that tubal scalars are ordered according to Definition 3.
In parallel to GSP theory, we can also define the Laplacianbased total variation and provide an ordering for the hypergraph Fourier basis of the Laplacian tensor L s . Hence, we define the Laplacian-based total variation on a hypergraph as Then, the total variation of a Fourier basis vector V j is In contrast to the Fourier basis obtained from the adjacency tensor, the eigenvectors of the Laplacian shifting operator are ordered from lowest frequency to highest as λ 1 ≤ λ 2 ≤ · · · ≤ λ j ≤ · · · ≤ λ N . In simple graphs, the eigenspace of the graph Laplacian L = VΛV T provides a similar notion of frequency as in classical Fourier analysis. The graph Laplacian eigenvectors associated with smaller eigenvalues vary slowly across the graph and the eigenvectors associated with larger eigenvalues oscillate more rapidly and are more likely to have dissimilar values on vertices connected by an edge [5]. For the hypergraph Fourier transform defined here, in Fig. 6, we visually compare the variation among three of the column vectors inV (1) which is the first frontal slice ofV. Note that there is higher variation among neighboring nodes as the value of the eigenvalue increases. Given the notion of frequency, we then introduce the definition the bandlimited signals on a hypergraph as follows.

Definition 11 (Band-limited Signal):
A hypergraph signal X s ∈ R N ×1×N s is a K-band-limited if its t-HGFT, X F s , has tubal scalars (x F s ) l,1 = 0 for all l > K where K ∈ 1, 2, . . . , N. The smallest K is called the bandwidth of X s .
Note that hypergraph band-limited signals are not required to be low-pass or smooth, in which the eigen-tuples are sorted from low to high frequency as explained before. Instead, the band-limited definition is equivalent to restricting the number of nonzero tubal scalars in the hypergraph Fourier spectrum X F s [29] with no specific ordering on the eigen-tuples.

V. DISCUSSIONS AND INTERPRETATIONS
In this section, we discuss the benefits of using t-product as the core operation for the proposed HGSP framework. We also establish the connection and advantages over classical GSP frameworks, focusing on the aspects that help to interpret the notion of hypergraph frequency.

A. Interpretation of Hypergraph Spectrum Space
To analyze the hypergraph spectrum space we consider not only the connection of HGSP with DSP and GSP but also the connection of the t-product with the Fourier transform.

1) Connection of the T-Product With the Fourier Transform:
To better understand the hypergraph spectrum of a signal, we consider the connection of the t-product with the discrete Fourier transform (DFT) [39]. As observed in Algorithm 1, the t-eigendecomposition is efficiently computed by traditional matrix eigendecomposition of the frontal slices of the shifting operator F s in the Fourier domainF s . Note that, as depicted in Fig. 7, a better interpretation of the hypergraph spectrum Fig. 7. The hypergraph spectrum space computed from the Laplacian tensor L s of the hypergraph in Fig. 5(top-left) is analyzed here. Better interpretability is achieved by considering the connection of the t-product with the DFT [39], which also provides a notion of frequency along the third dimension. Then, each frontal slice of the t-eigendecomposition in the Fourier domain depicted in this figure is efficiently computed by traditional matrix-decompositionL . Then, for each frontal slice, the notion of lower to higher frequency is provided by the increasing ordering of the eigenvalues as explained before in Section IV-A, and a notion of lower to higher frequency along the third dimension is given by the DFT of the tubal scalars.
can be obtained by analyzing the t-eigendecomposition in the Fourier domain. Then, for each frontal slice, the notion of frequency is provided by ordering the eigenvalues as explained in Section IV-A. Additionally, a notion of frequency along the third dimension is provided by the DFT of the tubal scalars.
Note that one of the important properties of the DFT is that when the signal is real and even, which is the case of the tubal scalars in the symmetrized tensors, its Fourier spectrum is also real and even. This can also be observed in Fig. 7 where the DC component is represented by k = 1, and the rest of frequencies, from lowest to highest, are given by the pair of frontal slices k = p + 1 ⇔ k = 2(N + 1) − p for p = 1, . . . , N, with k = N + 1 ⇔ k = N + 2 being the highest frequency. Here, we see that by the symmetrization of the shifting operator, we guarantee that its Fourier domain representation is real and even which also leads to an eigendecomposition with real eigenvalues, which is characteristic of undirected graphs and hence a desired property for undirected hypergraphs.

2) Connection With DSP and GSP:
To obtained an intuitive interpretation on the hypergraph frequency, we consider its relationship with DSP and GSP. In DSP, the discrete Fourier transform (DFT) of a signal s = {s n : [6]. The discrete frequencies ω k = 2πk/N, k = 0, 1, . . . N − 1 are related to the degree of variation of the spectral components. To observe this, we consider the definition of the total variation in DSP: where s [n−1] N is the one-time shifted signal and this shifting operation can be represented by multiplying by the circulant matrix: induces an ordering from lowest to highest frequencies on the eigenvalues as λ 0 , λ 1 , λ N −1 , λ 2 , λ N −2 , · · · , with the lowest frequency corresponding to λ 0 and the highest frequency corresponding to λ N/2 for even N or λ (N ±1)/2 for odd N , which is the conventional frequency ordering in DSP [45]. Hence, higher frequency components change faster over time, which implies larger total variation. Note that, the shift matrix A c can be seen as the adjacency matrix of a cycle line graph [6]. Thus, we now consider the total variation and frequency in GSP, replacing A c by the normalized shifting operator of any graph as where F norm = 1/|λ max |F. Then, λ m is larger than graph frequency λ l if the TV(v m ) > TV(v l ) since larger variation indicates faster change is the signal values between neighboring nodes on the graph, which represents higher frequencies. For undirected graphs, the eigenvalues are real numbers, and the frequency decreases with the increase of the increasing the eigenvalues as demonstrated for the general case of hypergraphs in Section IV-A. If the graph is undirected, the eigenvalues are complex and the frequency changes following a similar pattern as explained before for DSP. Now, we consider the generalization to HGSP by replacing the shifting operator by an M -th order tensor F s . Then, the generalized TV on hypergraphs in (21) measures the variation on the eigen-tensors components as As explained before, given that we consider undirected hypergraphs and use symmetrized shifting operators, the eigen-pairs are real in both the spatial domain λ j , V j and the Fourier domainλ j ,V j . Then, as a generalization of GSP and DSP, λ m is larger than hypergraph frequency λ l if the TV( V m ) > TV( V l ) since larger variation indicates faster changes in the signal values between neighboring nodes on the hypergraph, which represents higher frequencies. To illustrate this, consider the hypergraph example in Fig. 8 where, clearly, there is higher variation between the eigen-tensor V j and its shifted version F norm s * V j as the eigenvalues increase. Note that, we analyze this in the Fourier domain given that the ordering of the eigen-tuples (Definition 3) is better analyzed by considering the connection of the t-product with the Fourier transform.

B. Connections to Other Frameworks
1) Graph Signal Processing: HGSP aims at being a more general framework than traditional GSP, enabling the processing of signals in high-dimensional graphs that capture polyadic relationships among the nodes. Thus, the proposed HGSP framework is a generalization of traditional GSP which is a special case of HGSP with M = 2. Consider this particular case (M = 2) in the core elements of the proposed framework: r Algebraic descriptor: GSP considers the case in which each edge connects exactly two nodes while HGSP considers hyperedges that can connect more than two nodes. A graph is a uniform hypergraph with maximum edge cardinality M = m.c.e(H) = 2 represented by an 2nd-order adjacency tensor A ∈ R N ×N with entries a p 1 ,p 2 = a p 2 ,p 1 = 1 if the node v p 1 is connected to node v p 2 and zero otherwise. Note that the symmetrization operation defined in Section III is applied for M > 2 which excludes the case of simple graphs given that they are already symmetric according to Definition 2. Thus, this 2nd-order adjacency tensor A ∈ R N ×N is the binary adjacency traditionally used at the core of GSP. r Spectral Properties: the t-eigendecomposition [42], being a generalization of traditional matrix algebra, reduces to traditional matrix eigendecomposition when M = 2. Hence, the t-HGFT is the same as the GFT when M = 2. HGSP being a generalization of GSP not only holds the same benefits as the traditional GSP framework but also provides additional advantages. First, the adjacency tensor can, without ambiguities, encode high-order relationships that are widespread in many real-word applications [1], [13], [14], [15]. Second, the high-order shifting operation can either jointly process multiple signals (Definition 8) or consider cross-node relationships (Definition 9), modeling the join effects of nodes within a hyperedge.
2) Hypergraph Signal Processing: Recently, the CP tensor decomposition was used to define a hypergraph signal processing framework [9]. This section compares the core elements of the CP-based HGSP framework [9] with the ones in our framework and states the benefits of using t-product tensor decompositions for HGSP. As shown in Fig 9(a), the HGSP framework in [9] uses the orthogonal-CP decomposition of the shifting operator F computed as F represents either the Adjacency tensor A or the Laplacian tensor L which are super-symmetric tensors. As depicted in Fig 9(b), [9] defines the spectrum of the signal x ∈ R N on the hypergraph asx = Vx where V = [f 1 , . . . , f N ] T . Then for the hypergraph signal defined as the hypergraph Fourier spacex can be computed as M − 1 times the Hadamard product ofx, that isx =x · · · x T . The drawback of this approach is that tensors in general do not have an exact CP symmetric orthogonal decomposition [37] as shown in the example in Fig. 3(a), and finding the rank R of a tensor is NP-hard. Furthermore, since usually the CP-rank R < N, the authors in [9] add additional vectors f r by using zero coefficients λ r for r = R + 1, . . . , N. These N − R additional vectors f r should meet the conditions of: (1) orthogonality f r ⊥ f i for r = i, (2) normalization with |f r | = 1, and (3) f T r x = 0, ∀r > R. However, f r vectors satisfying these conditions can only be added when the signal is band-limited with bandwidth W ≤ R. Thus, as illustrated in Fig 9(c) and in the example in Fig. 3(a), since x lives in R N , the hypergraph Fourier transform defined in [9] maps the signal to a subspace of R N when R < N. Hence, perfect recovery on the inverse HGFT of x ∈ R N is not guaranteed unless the signal is band-limited with bandwidth W ≤ R. Note that one might consider conditions different from [9] in order to obtain a loss-free HGFT when using the CP symmetric orthogonal decomposition. On the one hand, one could set the rank to N , which would likely lead to higher errors on the tensor decomposition since the rank is not optimal, but a loss-free HGFT would be obtained. On the other hand, one might ignore the third condition on the additional f r vectors, optimize the rank to find the best tensor approximation, and then add N − R additional vectors f r associated to zero λ r coefficients which are not necessarily orthogonal to the signal. For the latter, the notion of frequency would be lost on these additional bases, given that they are all assigned to zero eigenvalues.
In the framework introduced here, the t-eigendecomposition always provides an exact tensor factorization of the shifting operator. Particularly, as shown in Fig. 9(d), when this tensor is symmetrized F s , the t-eigendecomposition is composed by an orthogonal tensor of eigen-matrices V and an f-diagonal tensor of eigen-tuples Λ. This representation allows for a loss-free generalization of the basic elements in a HGSP framework: the hypergraph Fourier transform (t-HGFT), in Fig. 9(e), and the inverse hypergraph Fourier transform (t-iHGFT), in Fig. 9(f). Moreover, the HGFT and iHGFT, defined here, are more robust than those obtained in the CP-based approach, given that any t-product tensor decomposition algorithm could be used and lead to reproducible results. In contrast, the performance of the CP-based approach relies on selecting a suitable tensor decomposition algorithm. Then, the advantages of adopting t-eigendecompositions to define an HGSP framework with respect to previous approaches are compelling -(1) t-product decompositions preserve the intrinsic structures of tensors; (2) the high-dimensional nature of signals is preserved; (3) the orthogonal eigenbasis spans the whole signal space which allows for loss-free Fourier transforms [42] and the symmetric properties of the tensors together with the Fourier domain computation of t-products allow for computationally efficient calculations.
In the following sections, essential tools in DSP are formulated within the context of the proposed t-HGSP framework and the aforementioned advantages are demonstrated through simulations on different types of data.

VI. HYPERGRAPH FILTERS
Filtering is an essential tool used in signal processing applications as measurement noise is always present [46]. Filtering in standard DSP amplifies or attenuates the contributions of each of signal's frequency components asŷ(f ) =q(f )x(f ) whereq(f ) is the transfer function of the filter. Multiplication in the Fourier domain corresponds to convolution in the time domain, i.e, y(t) = (q x)(t) [5]. Notably, these fundamental concepts have a natural parallel formulation in the proposed hypergraph Fourier transform -frequency filtering can be directly generalized as are, respectively, the tubal scalars at frequency λ l of the output signal, Y s , the input signal, X s , and the filter response, Q s , in the frequency domain. When taking the inverse Fourier transform of Y s , the tubal scalars of Y s are given by which can be written in tensor-tensor product notation as In GSP, graph spectral filtering is used to implement wellknown filtering techniques designed via optimization formulations such as Gaussian smoothing, bilateral filtering, total variation filtering, anisotropic diffusion, and nonlocal means filtering [32]. Many of these filters, in fact, arise as solutions of regularized ill-posed inverse problems such as denoising, inpainting, and super resolution. Hypergraph filters can also be designed via optimization approaches [9]. As an example, let Y s be the second moments of the observed noisy signal. To enforce a priori information that the second moments of the clean signal X s are smooth with respect to the underlying hypergraph, we include a regularization term as where d = γe 1 , and D[d] is a diagonal tensor with d repeated down the diagonal and Φ( X s , F) is a function that measures the smoothness of the signal X s in the hypergraph. Well-known optimization problems with types of regularizations such as the Tikhonov regularization have been previously formulated and solved for third order tensors in the literature [42]. A numerical experiment demonstrating this approach is included in Section III of the Supplemental Material.

VII. SAMPLING ON HYPERGRAPHS
Sampling is one of the fundamental concepts in DSP. Standard sampling theory relies on concepts of frequency domain analysis, shift-invariant signals, and band-limitedness. A common problem in networks, usually described by graphs, is to determine which nodes play the most important role. Graph signal sampling then becomes an essential tool for a GSP framework. Traditional sampling theory and spectral graph theory have been combined to generalize Nyquist sampling principles on graphs [29], [47], [48], [49], [50], [51], [52], [53]. Suppose that we want to sample Q tubal scalars from the hypergraph signal X s ∈ R N ×1×N s , then the sampled signal X s (S) ∈ R Q×1×N s with Q < N and S = {s 1 , s 2 , . . . , s Q } ⊂ V (H) denotes the sequence of sampled nodes. For the sampling of these Q tubal scalars, the sampling operator is defined as a 3 rd -order tensor Ψ ∈ R Q×N ×N s whose first frontal slice is given by and the other frontal slices are all zeros. Now, to recover X s ∈ R N ×1×N s , lets define an interpolation operator, Φ ∈ R N ×Q×N s , such that where X * s ∈ R N ×1×N s approximates X s either exactly or approximately [44]. For perfect recovery Φ * Ψ must be equal to the identity tensor which is, in general, not possible since the multi-rank (Definition 7) of Φ * Ψ, ρ(Φ * Ψ) ∈ R 1×1×N s , has frontal entries ρ(Φ * Ψ) (k) ≤ Q < N. However, perfect recovery can be achieved for band-limited signals considering the following property.
Lemma 1: If X s is a K-band-limited signal then where In this lemma, it is implied that K frequency components carry all the information of the signal; hence the GSP theory on the perfect recovery of band-limited signals [29] can be generalized to hypergraphs by the following theorem.
Theorem 2: Let V [K] be the tensor formed by K lateral slices of V, i.e. the orthonormal set { V 1 , V 2 , . . . , V K }, and let the sampling operator Ψ ∈ R Q×N ×N s satisfy for i = 1, 2, . . . , 2N + 1. Since ρ(Ψ * V [K] ) (k) ≤ min(Q, K), the sample size Q has to be greater than or equal than the bandwidth K. The interpolation operator Φ = V [K] * U ∈ R N ×Q×N s with U * Ψ * V [K] = I K×K×N s , achieves perfect recovery, i.e. X s = Φ * Ψ * X s for any K-bandlimited signal.
Proof: See Appendix C. As in GSP, from Theorem 2, we observe that an arbitrary sampling operator may not lead to perfect recovery even for band-limited hypergraphs signals [29]. To have perfect recovery, the sampling operator Ψ must be a qualified sampling operator by satisfying the multirank condition (36). The optimization of the sampling pattern such that Ψ is qualified sampling operator and robust to noise is part of our future work. It is important to note that, for hypergraphs, a sampling operator Ψ that satisfies the multirank condition (36) is not limited to the sampling of tubal scalars as in 32, but instead a qualified sampling operator Ψ could be designed such that the sampled signal X s (S) is formed by the linear combination of the signals in X s . This type of sampling where only a linear combination of the original signal can be obtained is present in many applications such as compressive spectral imaging (CSI) systems which aim at capturing large volumes of spatio-spectral information of a scene of interest from a set of noisy undersampled observations [54], [55].

A. Sparse Hypergraph Signal Representation
Sparse signal representations assume that most of the information can be captured by just a few basis functions in an orthonormal basis set [56]. Moments and cumulants are commonly used in the statistical analysis of random variables [57]. However, storing empirical moment tensors can be prohibitively expensive. Here, we demonstrate that the proposed hypergraph Fourier Transform methodology can be successfully applied to obtain a sparse representation of a hypergraph signal carrying empirical moments (Definition 9). To this end, we consider the "Cameraman" image of size 128 × 128 which is first tiled into non-overlapping blocks of size 16 × 16. As in [9], for each block, a hypergraph with a 3rd-order adjacency tensor, A, is obtained by using the image adaptive neighborhood hypergraph (IANH) model [41]. The hypergraph signal, X s , was built from the one dimensional grayscale signal according to Definition 9, capturing its 2-D properties. The sparse representation of the hypergraph signal is then obtained by keeping only the largest C unique frequency coefficients in the hypergraph spectrum domain, X F s . We compare our results with the recently proposed CP-based HGSP framework [9]. Given that finding the rank R is an NP-hard problem and the signal is not necessarily bandlimited, R = N was used in this experiment for the symmetric orthogonal CP approximation [35]. Fig. 10(a) illustrates the significant performance gains attained by the proposed t-HGSP framework over competing state-of-the-art methods.
As in [9], we also test the efficiency of t-HGSP sparse signal representation over the MovieLens data sets [58] where each movie data point has rating scores and tags from viewers. For this numerical experiment, the scores of movies are considered as signals and the graph model is constructed based on the tag relationships, i.e., two or more movies are connected in a hyperedge if they have similar tags [9]. For convenience, we used a subset of 400 movies and we set M = 3. With the graph and hypergraph models, the sparse hypergraph signal representation is obtained using the same method explained above. Results depicted in Fig. 10(b) demonstrate that the proposed t-HGSP framework outperforms competing state-of-the-art methods.
Additionally, we consider a hyper-spectral point cloud sensed by NASA's Goddard's LiDAR, Hyperspectral and Thermal (G-LiHT) airborne imaging system [26], [59]. A hypergraph with a 3rd-order adjacency tensor A was then built based on the similarities of the elements in the point cloud (N = 882). The hypergraph signal X s was built from the set of 1-D signals in the hyperspectral point cloud according to Definition 8. The sparsified signal was obtained as before. To compare our results to the CP hypergraph-based approach [9], [35] and GSP [45], a sparse signal representation for each 1-D signal is obtained independently. The simple graph [45] is also built based on similarities. Fig. 10(c) illustrates large performance gains attained by the t-HGSP framework over state-of-the-art methods which demonstrates the benefit of jointly processing multiple signals.

B. Spectral Clustering
Clustering is key in many applications such as social network analysis, computer vision, and computational biology [60]. In [K] =V [K] (:, :, 1) ∈ R N ×K be the first frontal slice of the tensorV [K] that contains the K eigen-matrices in the Fourier domain. 4 : For i = 1, . . . , N, let y [K] . 5: Cluster the N nodes with features vectors (y i ) i=1,...,N in R K with the k-means algorithm into K clusters {C 1 , . . . , C K }. 6: Use the clustering result to partition the hypergraph into S 1 , . . . , S K such as S i = {v j |y j ∈ C i }, i.e. the j-th node belongs to the i-th partition if y j ∈ C i . recent years, spectral clustering has aroused more and more attention due to its good clustering performance and solid theoretical foundation [61]. However, standard spectral clustering methods consider only pairwise connections between nodes and there are many applications in which interactions involving more than two nodes occur. Hence, it is fundamental to extend spectral clustering to hypergraphs. The definition of a suitable spectral space is key when developing a hypergraph spectral clustering. Li et al. [62] and Ahn et al. [60] proposed a hypergraph spectral clustering approach based on the eigenspace of what they called the hypergraph processed similarity matrix. Zhang et al. [9] introduced a hypergraph spectral clustering approach based on the hypergraph Fourier space given by the symmetric orthogonal CP decomposition. In contrast, we generalize graph spectral clustering techniques to our t-HGSP framework by means of the hypergraph Fourier space of the Laplacian tensor (Algorithm 2). As in [9], we used the zoo dataset [63] to create an M = 3 hypergraph and compare our results with those from the hypergraph similarity method (HSC) in [60] and CP-based HGSP [9]. The average Silhouette [64] of nodes was used to asses clustering quality among the different methods. This metric is particularly useful if the ground truth labels are not known. The silhouette value for each point is a measure of how similar that point is to points in its own cluster, when compared to points in other clusters. The silhouette value ranges from -1 to 1. A high average silhouette indicates a good clustering. Given that the ground truth labels are known, we also used the Normalized Mutual Information (NMI), weighted F1 score, and Jaccard index to compare the clustering performance [65]. In addition to the already mentioned advantages of our framework, in Figs. 10(d) and 12(a), it is shown that the proposed hypergraph Fourier space of the Laplacian tensor is not only suitable for hypergraph spectral clustering but also outperforms the prior art. Additionally, to demonstrate the performance of the proposed framework on higher-order (M > 3) hypergraphs, we consider different subsets of two real-world dataset, the coauthorship networks, cora and DBLP, as depicted in Fig. 11 . In this hypergraph, the nodes represent documents and all documents co-authored Nodes represent documents and all documents co-authored by an author are in one hyperedge [66]. The labels classify nodes (papers) into categories such as "algorithms" and "data mining". Hyperedges' cardinality is color-coded: (red) |e| = 4, (green) |e| = 3, and (blue) |e| = 2. by an author are in one hyperedge [66]. The labels classify nodes (papers) into categories such as "algorithms", "data mining", and so on. For the higher-order orthogonal symmetric CP decomposition, we used the algorithm proposed in [35]. As observed in Fig. 12(b)-(d), the proposed approach outperforms HSC [60] and CP HGSP [9] when high-order hypergrahs are consider. More details about the spectral clustering algorithm, metrics, and other results can be found in the supplemental material.

C. Denoising
Hypergraph filters can be used to denoise a hypergraph signal. Let z = x + n ∈ R N be a one dimensional noisy signal with noise n and Z s ∈ R N ×1×N s its corresponding hypergraph signal. Assuming that the original signal is smooth on the graph, a simple approach to denoising is by spectrum shrinkage where filtering reduces to hard-thresholding of the spectral coefficients [67] such that  where λ c is a cut-off frequency. Then, the denoised signal is X s = Q s * Z s . Since the frequency response of a clean signal is expected to be dominated by low-frequency components, the high-frequency components are most likely to be generated by noise. This noise is thus attenuated by hard-thresholding the transform coefficients [67]. To test, the performance of our approach, we conducted a denoising experiment on the "Cameraman" image. The experimental setup is similar to the one in the compression experiment, except for noise addition and the use of overlapping windows. The same hard-thresholding filter design was used for GSP and CP-based HGSP denoising. We contaminated the signal with Gaussian noise. The cut-off frequency was optimized for the best performance in all the cases. The proposed hypergraph framework outperforms CPbased HGSP and GSP denoising as shown in Table I. Even though almost no improvement over GSP is observed on the performance metrics, as shown in Fig. 13, the details of the image are better recovered by the proposed approach. A bigger section of the image is included in Fig. 14 of the Supplemental material.

IX. CONCLUSION
A novel tensor-based t-HGSP framework using t-product decompositions is introduced. The proposed framework not only generalizes traditional GSP to high-dimensional hypergraphs but more importantly is loss-free, preserves the dimensionality of the signals, and is less computationally complex when compared to the state-of-the-art. The core elements of the new t-HGSP framework have been laid down, including the concept of a hypergraph signal, hypergraph shifting, frequency analysis, and band-limited signals. Additionally, fundamental tools such as sampling and filtering were introduced and demonstrated experimentally. Notably, the proposed framework outperforms state-of-the-art by considering high-order properties of one-dimensional signals or jointly processing a set of signals, as demonstrated in the denoising and sparse signal representation experiments. Additionally, natural hypergraphs such as co-authorship networks can be exploited for applications such as classification and clustering without being a map to a low-dimensional graph. However, the key to the success of hypergraph-based methods is having a meaningful hypergraph that may only be readily available for some applications. Hence, our current work is focused on learning optimal hypergraph topologies from data that could further boost the performance of the proposed t-HGSP tools.
Given the proposed t-HGSP framework, many opportunities emerge to continue exploring HGSP and its applications. For instance, further interpretations on the t-HGSP spectrum could be drawn by considering its connection with manifolds and the connection of GSP with other areas such as information theory. Moreover, other tensor decompositions such as HOSVD [68], tensor-train [69], and O-SVD [70] could be potentially used to define other HGSP frameworks with different properties. Our future work is focused on optimizing sampling patterns for different sampling operator configurations, hypergraph neural networks, scalability, and generalizing hypergraph's shifting operators to tensors that are not symmetric under the t-product algebra.

APPENDIX A PROOF OF THEOREM 1
Let F s = A s and λ i and λ j be two eigen-tuples of F s with corresponding eigen matrices V i and V j . Given that F s is symmetric, it has real eigen-tuples. If the eigen-tuples are ordered as λ i > λ j , then the total variations of their eigenvectors satisfy TV H ( V i ) < TV H ( V j ).

APPENDIX B PROOF OF LEMMA 1
If X s is a K-band-limited signal then where V [K] = [ V 1 , . . . , V K ] ∈ R N ×K×2N +1 and ( X F s ) [K] ∈ R K×1×2N +1 correspond to the K nonzero elements of the hypergraph signal in the frequency domain X F s .
Proof: Since X s is band-limited, then (x F s ) l,1 = V T l * X s = 0 when l > K. Then, we have that APPENDIX C PROOF OF THEOREM 2 Let V [K] be the tensor formed by K lateral slices of V, i.e. the orthonormal set { V 1 , V 2 , . . . , V K }, and let the sampling operator Ψ ∈ R Q×N ×N s satisfy for i = 1, 2, . . . , 2N + 1. Since ρ(Ψ * V [K] ) (k) ≤ min(Q, K), the sample size Q has to be greater than or equal than the bandwidth K. The interpolation operator Φ = V [K] * U ∈ R N ×Q×N s with U * Ψ * V [K] = I K×K×N s , achieves perfect recovery, i.e. X s = Φ * Ψ * X s for any K-bandlimited signal. Proof: To proof the theorem, we show that P = Φ * Ψ is a projector operator and that the range of Φ is the space of K-bandilimited signals as in [44]. Since ρ(Ψ * V [K] ) (k) = K for i = 1, 2, . . . , 2N + 1 and ρ(U * Ψ * V [K] ) (k) = K for i = 1, 2, . . . , 2N + 1, then ρ(U) (k) = K for i = 1, 2, . . . , 2N + 1. Thus, the interpolation operator Φ = V [K] * U ∈ R N ×M ×N s spans the set of K-bandlimited signals. To prove that P is a projector, we show that P 2 = P * P = P as P * P = Φ * Ψ * Φ * Ψ, where (a) follows from U * Ψ * V [K] = I K×K×N s . Notice that when M = K, U is the inverse of Ψ * V [K] ; and when M > K, U is the pseudo-inverse of Ψ * V [K] as in the traditional matrixbased sampling theory.