{'fraction_non_alphanumeric': 0.04083230985564515, 'fraction_numerical': 0.016660436802477707, 'mean_word_length': 3.9442928804013024, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 44, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Topological Signal Processing over Simplicial Complexes

The goal of this paper is to establish the fundamental tools to analyze signals defined over a topological space, i.e. a set of points along with a set of neighborhood relations. This setup does not require the definition of a metric and then it is especially useful to deal with signals defined over non-metric spaces. We focus on signals defined over simplicial complexes. Graph Signal Processing (GSP) represents a very simple case of Topological Signal Processing (TSP), referring to the situation where the signals are associated only with the vertices of a graph. We are interested in the most general case, where the signals are associated with vertices, edges and higher order complexes. After reviewing the basic principles of algebraic topology, we show how to build unitary bases to represent signals defined over sets of increasing order, giving rise to a spectral simplicial complex theory. Then we derive a sampling theory for signals of any order and emphasize the interplay between signals of different order. After having established the analysis tools, we propose a method to infer the topology of a simplicial complex from data. We conclude with applications to real edge signals to illustrate the benefits of the proposed methodologies.


I. INTRODUCTION
Historically, signal processing has been developed for signals defined over a metric space, typically time or space. More recently, there has been a surge of interest to deal with signals that are not necessarily defined over a metric space. Examples of particular interest are biological networks, social networks, etc. The field of graph signal processing (GSP) has recently emerged as a framework to analyze signals defined over the vertices of a graph [2], [3]. A graph G(V, E) is a simple example of topological space, composed of a set of elements (vertices) V and a set of edges E representing pairwise relations. However, notwithstanding their enormous success, graph-based representations are not always able to capture all the information present in complex interaction networks, as suggested in [4]- [7]. The main reason is that many complex interactions, like those occurring for example in metabolic networks, cannot be reduced to pairwise interactions but require a full description of the rich multiway relations that take place [4]. It is then necessary to go beyond graphs to fully capture these more complex interaction mechanisms.
Hypergraphs offer a framework that helps to overcome such limitations. In very general terms, a topological space is composed of a set V of elements along with an ensemble of multiway relations, represented by a set S containing subsets of various order of the elements of V. The structure H(V, S) is known as a hypergraph. In particular, a class of hypergraphs that is particularly appealing for its rich algebraic structure is given by simplicial complexes, whose defining feature is the inclusion property stating that if a set A belongs to S, then all subsets of A also belong to S. Restricting the attention to simplicial complexes is not a strong limitation, as they are general enough to include most cases of interest, while at the same time exhibit a rich algebraic structure that can be properly exploited. Learning models based on simplicial complexes have been already proposed in brain network analysis [6], neuronal morphologies [8], co-authorship networks [9], collaboration networks [10], [11], tumor progression analysis [12]. Indeed, the use of algebraic topology tools for the extraction of information from data is not new: The framework known as Topological Data Analysis (TDA), see e.g. [13], has exactly this goal. Interesting applications of algebraic topology tools have been proposed to control systems [14], statistical ranking from incomplete data [15], [16], distributed coverage control of sensor networks [17]- [19], wheeze detection [20]. One of the fundamental tools of TDA is the analysis of persistent homologies extracted from data [21], [22]. Topological methods to analyze signals and images are also the subjects of the two books [23] and [24]. The goal of our paper is to establish a fundamental framework to analyze signals defined over a simplicial complex. Our approach is complementary to TDA: Rather than focusing, like TDA, on the properties of the simplicial complex extracted from data, we focus on the properties of signals defined over a simplicial complex. Our approach includes GSP as a particular case: While GSP focuses on the analysis of signals defined over the vertices of a graph, our topological signal processing (TSP) framework considers signals defined over simplices of various order, i.e. node signals, defined over points, edge signals, defined over pairs of points, signals defined over triplets of points, etc. Relevant examples of signals defined over pairs of points are flow signals, like blood flows between different areas of the brain [25], data traffic over communication links [26], regulatory signals in gene regulatory networks. It is known for example that the dysregulation of these regulatory signals is one of the causes of cancer [27]. Examples of signals defined over triplets are co-authorship networks, where a triangle indicates a triplet of authors that have at least a common publication [9] and the associated signal value is the number of such publications. There are previous works dealing with the analysis of edge signals, like [28]- [31]. More arXiv:1907.11577v1 [eess.SP] 26 Jul 2019 specifically, in [28] the authors introduced a class of filters to analyze edge signals based on the edge-Laplacian matrix [32]. A graph-based semi-supervised learning method for learning edge flows was then suggested in [29]. Other works analyzed edge signals using a line-graph transformation [30], [31]. However, all these approaches essentially use tools derived from a graph-based description of the domain. Conversely, we want to show that, while in the analysis of node signals there is no loss of information in using tools derived from the associated graph (i.e., the simplicial complex of order one), when we analyze edge (flow) signals, we need to consider the underlying simplicial complex of order two, which contains vertices, edges, and triangles. Furthermore, from the inference point of view, while when we observe a set of node signals, it makes perfect sense to infer the topology of the graph capturing the relations among those signals, see e.g. [33]- [35], when we observe flow signals, we need to infer the topology of a simplicial complex of order two, and so on. We published some preliminary results of our work in [1]. Here, we extend the work of [1], deriving a sampling theory for signals defined over complexes of various order and proposing new inference methods, more robust against noise. In summary, our contributions are listed below: 1) we build a spectral simplicial theory, as an extension of spectral graph theory, based on topology-aware definition of suitable dictionaries to represent signals defined over simplicial complexes of various order; 2) we derive a sampling theory defining the conditions for the recovery of high order signals from a subset of observations, highlighting the interplay between signals of different order; 3) we propose inference algorithms to extract the structure of the simplicial complex from high order signals.
The paper is organized as follows. Section II recalls the main algebraic principles that will be the basis for the derivation of the signal processing tools in the ensuing sections. In Section III, we will first recall the eigenvectors properties of higher-order Laplacian and the Hodge decomposition. Then, we provide a simplicial complex spectral theory for designing unitary bases to represent edge signals. In Section IV we illustrate some methods for recovering the edge signal components from noisy observations. Section V provides theoretical conditions to recover the entire edge signal from a subset of samples. Some methodologies to infer the simplicial complex structure from noisy observations are then proposed in Section VI. Section VII assesses the performance of the proposed methodologies over synthetic and real data. Finally, Section VIII draws some conclusions.

II. REVIEW OF ALGEBRAIC TOPOLOGY TOOLS
In this section we recall the basic principles of algebraic topology [36] and discrete calculus [37], as they will form the background required for deriving the basic signal processing tools to be used in later sections. We follow an algebraic approach that is accessible also to readers without a specific background on algebraic topology.

A. Discrete domains: Simplicial complexes
Every k-simplex has exactly k + 1 faces. An abstract simplicial complex X is a finite collection of simplices that is closed under inclusion of faces, i.e., if σ i ∈ X , then all faces of σ i also belong to X . The dimension of a simplex is one less than its cardinality. Then, a vertex is a 0-dimensional simplex, an edge has dimension 1, and so on. The dimension of a simplicial complex is the largest dimension of any of its simplices. A graph is a particular case of an abstract simplicial complex of order 1, containing only simplices of order 0 (vertices) and 1 (edges).
If the set of points is embedded in a real space R D of dimension D, we can associate a geometric simplicial complex with the abstract complex. A set of points in a real space R D is affinely independent if it is not contained in a hyperplane; an affinely independent set in R D contains at most D + 1 points. A geometric k-simplex is the convex hull of a set of k + 1 affinely independent points, called its vertices. Hence, a point is a 0-simplex, a line segment is a 1-simplex, a triangle is a 2-simplex, a tetrahedron is a 3-simplex, and so on. A geometric simplicial complex shares the fundamental property of an abstract simplicial complex: It is a collection of simplices that is closed under inclusion and with the further property that the intersection of any two simplices in X is also a simplex in X , assuming that the empty set is an element of every simplicial complex. Although geometric simplicial complexes are easier to visualize and, for this reason, we will often use geometric terms like edges, triangles, and so on, as synonyms of pairs, triplets, we do not require the simplicial complex to be embedded in any real space, so as to leave the treatment as general as possible.
The structure of a simplicial complex is captured by the neighborhood relations of its subsets. As with graphs, it is useful to introduce first the orientation of the simplices. Every simplex, of any order, can have only two orientations, corresponding to the permutations of its elements. Two orientations are equivalent if each of them can be recovered from the other by an even number of permutations. A k- . . , v i k } of order k, together with an orientation is an oriented k-simplex and is denoted by Two simplices of order k, σ k i , σ k j ∈ X , are upper adjacent in X , if both are faces of a simplex of order k + 1. Two simplices of order k, σ k i , σ k j ∈ X , are lower adjacent in X , if both have a common face of order k − 1 in X . A (k − 1)-face σ k−1 j of a k-simplex σ k i is called a boundary element of σ k i . We use the notation σ k−1 j ⊂ σ k i to indicate that σ k−1 j is a boundary element of σ k i . Given a simplex σ k−1 j ⊂ σ k i , we use the notation σ k−1 j ∼ σ k i to indicate that the orientation of σ k−1 j is coherent with that of σ k i , whereas we write σ k−1 j σ k i to indicate that the two orientations are opposite.
For each k, C k (X , R) denotes the vector space obtained by the linear combination, using real coefficients, of the set of oriented k-simplices of X . In algebraic topology, the elements of C k (X , R) are called k-chains. If {σ k 1 , . . . , σ k n k } is the set of k-simplices in X , a k-chain τ k can be written as τ k = n k i=1 α i σ k i . Then, given the basis {σ k 1 , . . . , σ k n k }, a chain τ k can be represented by the vector of its expansion coefficients (α 1 , . . . , α n k ). An important operator acting on ordered chains is the boundary operator. The boundary of the ordered k- So, for example, given an oriented triangle i.e., a suitable linear combination of its edges. It is straightforward to verify, by simple substitution, that the boundary of a boundary is zero, i.e., ∂ k ∂ k+1 = 0.
It is important to remark that an oriented simplex is different from a directed one. As with graph, an oriented edge establishes which direction of the flow is considered positive or negative, whereas a directed edge only permits flow in one direction [37]. In this work we will considered oriented, undirected simplices.

B. Algebraic representation
The structure of a K-complex is fully described by the set of its incidence matrices B k , k = 1, . . . , K. Given an orientation of the simplicial complex X , the entries of the incidence matrix B k establish which k-simplices are incident to which (k − 1)-simplices. Formally speaking, its entries are defined as follows: If we consider, for simplicity, a simplicial complex of order two, composed of a set V of vertices, a set E of edges, and a set T of triangles, having cardinalities V = |V|, E = |E|, and T = |T |, respectively, we need to build two incidence matrices B 1 ∈ R V ×E and B 2 ∈ R E×T . From (1), the property that the boundary of a boundary is zero maps into the following matrix form The structure of a K-simplicial complex is fully described by its high order combinatorial Laplacian matrices, of order k = 0, . . . , K, defined as ... ...
It is worth emphasizing that all Laplacian matrices of intermediate order, i.e. k = 1, . . . , K − 1, contain two terms: The first term, also known as lower Laplacian, expresses the lower adjacency of k-order simplices; the second terms, also known as upper Laplacian, expresses the upper adjacency of k-order simplices. So, for example, two edges are lower adjacent if they share a common vertex, whereas they are upper adjacent if they are faces of a common triangle. Note that in graphs, vertices can only be upper adjacent, if they are incident to the same edge. This is why the graph Laplacian L 0 contains only one term.

III. SPECTRAL SIMPLICIAL THEORY
In this paper, we are interested in analyzing signals defined over a simplicial complex. Given a set S, a signal is defined as a mapping of the form The cardinality of the elements of S induces the order of the signal. Even though our framework is general, in many cases we focus on simplices of order up to two. In that case, we consider a set of vertices V, a set of edges E and a set of triangles T , of dimension V , E, and T , respectively. The signals over each complex of order k, with k = 0, 1 and 2, are defined as the following maps: Spectral graph theory represents a solid framework to extract features of a graph looking at the eigenvectors of the combinatorial Laplacian L 0 of order 0. The eigenvectors associated with the smallest eigenvalues of L 0 are very useful, for example, to identify clusters [38]. Furthermore, in GSP it is well known that a suitable basis to represent signals defined over the vertices of a graph, i.e. signals of order 0, is given by the eigenvectors of L 0 . In particular, given the eigendecomposition of L 0 : the Graph Fourier Transform (GFT) of a signal s 0 over an undirected graph has been defined as the projection of the signal onto the space spanned by the eigenvectors of L 0 , i.e. (see [3] and the references therein) Equivalently, a signal defined over the vertices of a graph can be represented as From graph spectral theory, it is well known that the eigenvectors associated with the smallest eigenvalues of L 0 encode information about the clusters of the graph. Hence, the representation given by (12) is particularly suitable for signals that are smooth within each cluster, whereas they can vary arbitrarily across different clusters. For such signals, in fact, the representation in (12) is sparse or approximately sparse. As a generalization of the above approach, we may represent signals of various order over bases built with the eigenvectors of the corresponding high order Laplacian matrices, given in (6)- (8). Hence, using the eigendecomposition we may define the GFT of order k as the projection of a korder signal onto the eigenvectors of L k , i.e.
so that a signal s k can be represented in terms of its GFT coefficients as Now we want to show under what conditions (15) is a meaningful representation of a k-order signal and what is the meaning of such a representation. More specifically, the goal of this section is threefold: i) we recall first the relations between eigenvectors of various order of L k ; ii) we recall Hodge decomposition, which is a basic theory showing that the eigenvectors of any order can be split into three different classes, each representing a specific behavior of the signal; iii) we provide a theory showing how to build a topology-aware unitary basis to represent signals of various order starting only from topological properties.

A. Relations between eigenvectors of different order
There are interesting relations between the eigenvectors of Laplacian matrices of different order [39], which is useful to recall as they play a key role in spectral analysis. The following properties hold true for the eigendecomposition of k-order Laplacian matrices, with k = 1, . . . , K − 1.
Proposition 1: Given the Laplacian matrices L k of any order k, with k = 1, . . . , K − 1, it holds: 1) the eigenvectors associated with the nonnull eigenvalues of B T k B k are orthogonal to the eigenvectors associated with the nonnull eigenvalues of B k+1 B T k+1 and viceversa; 2) if v is an eigenvector of B k B T k associated with the eigenvalue λ, then B T k v is an eigenvector of B T k B k , associated with the same eigenvalue; 3) the eigenvectors associated with the nonnull eigenvalues λ of L k are either the eigenvectors of B T k B k or those of B k+1 B T k+1 ; 4) the nonnull eigenvalues of L k are either the eigenvalues of B T k B k or those of B k+1 B T k+1 . Proof. All above properties are easy to prove. Property 1) is straightforward: because of (4). Similarly, for the converse. Property 2) is also straightforward: If v is an eigenvector of B k B T k associated with a nonvanishing eigenvalue λ, then Finally, properties 3) and 4) follow from the definition of k- and from property 1).
Remark: Recalling that the eigenvectors associated with the smallest nonzero eigenvalues of L 0 are smooth within each cluster, applying property 2) to the case k = 1, it turns out that the eigenvectors of L 1 associated with the smallest eigenvalues of B T 1 B 1 are approximately null over the links within each cluster, whereas they assume the largest (in modulus) values on the edges across clusters. These eigenvectors are then useful to highlight inter-cluster edges.

B. Hodge decomposition
Let us consider the eigendecomposition of the k-th order Laplacian The structure of L k , together with the property B k B k+1 = 0, induces an interesting decomposition of the space R D k of signals of order k of dimension D k . First of all, the property can be decomposed into two parts: one belonging to im(B k+1 ) and one orthogonal to it. Furthermore, recalling that the whole space R D k can always be decomposed where the vectors in ker(L k ) are also in ker(B k ) and ker(B T k+1 ). This implies that, given any signal s k of order k, there always exist three signals s k−1 , s k H , and s k+1 , of order k − 1, k, and k + 1, respectively, such that s k can always be expressed as the sum of three orthogonal components: This decomposition is known as Hodge decomposition [40] and it is the extension of the Hodge theory for differential forms on Riemannian manifold to simplicial complexes. The subspace ker(L k ) is called harmonic subspace since each s k ∈ ker(L k ) is a solution of the discrete Laplace equation When embedded in a real space, a fundamental property of geometric simplicial complexes of order k is that the dimensions of ker(L k ), for k = 0, . . . , K are topological invariants of the K-simplicial complex, i.e. topological features that are preserved under homeomorphic transformations of the space. The dimensions of ker(L k ) are also known as Betti numbers β k of order k: β 0 is the number of connected components of the graph, β 1 is the number of holes, β 2 is the number of cavities, and so on [40].
The decomposition in (20) shows an interesting interplay between signals of different order, which we will exploit in the ensuing sections. Before proceeding, it is useful to clarify the meaning of the terms appearing in (20). Let us consider, for simplicity, the analysis of flow signals, i.e. the case k = 1. To this end, it is useful to introduce the curl and divergence operators, in analogy with their continuous time counterpart operators applied to vector fields. More specifically, given an edge signal s 1 , the (discrete) curl operator is defined as This operator maps the edge signal s 1 onto a signal defined over the triangle sets, i.e. in R T , and it is straightforward to verify that the generic i-th entry of the curl(s 1 ) is a measure of the flow circulating along the edges of the i-th triangle.
Similarly, the (discrete) divergence operator maps the edge signal s 1 onto a signal defined over the vertex space, i.e. R V , and it is defined as div(s 1 ) = B 1 s 1 .
Again, by direct substitution, it turns out that the i-th entry of div(s 1 ) represents the net-flow passing through the i-th vertex, i.e. the difference between the inflow and outflow at node i. Thus a non-zero divergence reveals the presence of a source or sink node. If we consider equation (20) in the case k = 1, recalling that B 1 B 2 = 0, it is easy to check that the first component in (23) has zero curl, and then it may be called an irrotational component, whereas the third component has zero divergence, and then it may be called a solenoidal component, in analogy to the calculus terminology used for vector fields.
The harmonic component s 1 H is a flow vector that is both curl-free and divergence-free. Notice also that, in (23), B T 1 s 0 represents the (discrete) gradient of s 0 .

C. Topology-aware unitary basis
In this section, we propose a method for building a unitary basis to represent signals of any order, capturing the structure of the underlying topology.
The approach we follow here is a generalization of the approach we proposed in [41] for graphs. We start reviewing the approach of [41], for the sake of clarity, and then we generalize it. Let us start from the graph reported in Fig. 1(a). A function capturing the connectivity properties of a graph is the cut size. Consider first a partition of the vertex set V into a pair of sets The cut size is defined as is a set function and it is known for being a submodular function [42]. Hence, F 0 (A 0 , A 1 ) can be transposed onto a function defined over R V through the so called Lovász extension [42]. For undirected graphs, the Lovász extension of the cut size is [42]: with x 0 ∈ R V . A possible way to build a unitary basis for signals defined on the vertices of a graph is to solve the following optimization problem The objective function to be minimized is convex, but the above problem is non-convex because of the unitary constraint. To simplify the search of the basis, we can relax the objective function to become Substituting (27) in (26), we still have a non-convex problem. However, its solution is known to be given by the eigenvectors of L 0 . From this perspective, the GFT can be interpreted as the projection onto the basis found as a solution of the relaxed problem.
Let us consider now the partition of V in three sets (A 0 , A 1 , A 2 ), and the simplicial complex of order 2, sketched in Fig.1 (b). The cut size is now defined as where a ijk = 1 if (i, j, k) ∈ T and a ijk = 0 otherwise. Our main result, stated in the following theorem, is that the Lovász extension of the triangle-cut function gives rise to a measure of edge signal variation along triangles.
The function in (29) represents the sum of the absolute values of the curls over all the triangles. A suitable unitary basis for representing rotational edge signals can then be found as the solution of the following problem Similarly to what we did with graphs, we relax f 1 (x 1 ) and substitute it with Substituting this function back to (30), the solution is given by the eigenvectors of B 2 B T 2 . However, since the kernel of B 2 B T 2 may be large, and recalling Hodge decomposition, it is useful to take as basis the eigenevectors of L 1 , which contains the eigenvectors associated to the nonnull eigenvalues of B 2 B T 2 and of B T 1 B 1 , plus the eigenvectors in the kernel of L 1 . It is worth noticing, by property 2) of Prop. 1, that the eigenvectors associated to the nonnull eigenvalues of B T 1 B 1 are equal to B T 1 u 0 i , where u 0 i are the eigenvectors associated with the nonnull eigenvalues of L 0 .

IV. EDGE FLOWS ESTIMATION
Let us consider now the observation of a flow signal affected by additive noise. Our goal is to recover the solenoidal, irrotational and harmonic signal components from the noisy data. The observation vector is then where v 1 is noise. Let us suppose, for simplicity, that the noise vector is Gaussian, with zero-mean entries all having the same variance σ 2 n . The optimal estimator can then be formulated as the solution of the following problem Note that problem Q is convex. Then, there exists multipliers (note that Slater's constraint qualification is satisfied). The associated Lagrangian function is Exploiting the orthogonality property B 1 B 2 = 0, it is easy to get the following KKT conditions Multiplying both sides of condition (c) by B T 2 , and using the second equality in (d) and condition (b), we get This means that s 2 and λ 2 may differ only by an additive vector lying in the nullspace of Similarly, multiplying (c) by B 1 and using the first equality in (d) and condition (a), we obtain which implies that s 0 = λ 1 +c 1 , with c 1 such that B 1 B T 1 c 1 = 0. Thus, condition (c) reduces to which says, as expected, that we can derive the harmonic component by subtracting the estimated solenoidal and irrotational parts from the observed flow signal x 1 . To recover the irrotational flow s 1 irr from the 0-order signal s 0 we need to solve equation (a) in (35). Note that L 0 is not invertible. For connected graphs, it has rank V −1 and its kernel is the span of the vector 1 of all ones. However, since the vector b = B 1 x 1 is also orthogonal to 1, the normal equation L 0 s 0 = B 1 x 1 admits the nontrivial solution (at least for connected graphs): Similarly, the 2-order signal s 2 , solution of the second equation in (35), can be obtained as The irrotational, solenoidal and harmonic components can then be recovered as followŝ Note that the first two conditions in (35) imply that the variables s 0 ,s 2 in Q are indeed decoupled so that the optimal solutions coincide with those of the following problemŝ V. SAMPLING AND RECOVERING OF SIGNAL DEFINED OVER SIMPLICIAL COMPLEX Suppose now that we only observe a few samples of a k-order signal. The question we address here is to find the conditions to recover the whole signal s k from a subset of samples. To answer this question, we may use the theory developed in [43] for signals on graph, and later extended to hypergraphs in [44]. For simplicity, we focus on signals defined over a simplicial complex of order K = 2, i.e. on vertices, edges and triangles. Given a set of edges S ⊆ E we define an edge-limiting operator as a diagonal matrix D S of dimension equal to the number of edges, with a one in the positions where we measure the flow, and zero elsewhere, i.e.
where 1 S is the set indicator vector whose i-th entry is equal to one if the edge e i ∈ S and zero otherwise. We say that an edge signal s 1 is perfectly localized over the subset S ⊆ E (or S-edge-limited) if s 1 = D S s 1 . Similarly, given the matrix U 1 whose columns are the eigenvectors of L 1 , and a subset of indices F, we define the operator where Σ F = diag(1 F ). An edge signal s 1 is |F|-bandlimited over a frequency set F if F F s 1 = s 1 . The operators D S and F F are self-adjoint and idempotent and represents orthogonal projectors, respectively, on the sets S and F. If we look for edges signals which are perfectly localized in both the edge and frequency domains, some conditions for perfect localization have been derived in [43]. More specifically: a) s 1 is perfectly localized over both the edge set S and the frequency set F if and only if the operator F F D S F F has an eigenvalue equal to one, i.e.
where A 2 denotes the spectral norm of A; b) a sufficient condition for the existence of such signals is that |S| + |F| > E. Conversely, if |S|+|F| ≤ E, there could still exist perfectly localized vectors, when condition (44) holds. In the following we first consider sampling on a single layer by extracting samples from the edge signals and then, we consider multi-layer processing using samples of signals defined over different layers of the simplex.

A. Single-layer sampling
In the following theorem [43] we provide a necessary and sufficient condition to recover the edge signal s 1 from its samples s 1 S D S s 1 .
Theorem 2: Given the bandlimited edge signal s 1 = F F s 1 , it is possible to recover s 1 from a subset of samples collected over the subset S ⊆ E if and only if the following condition holds: Proof. The proof is a straightforward extension of Th. 4.2 in [43] to signals defined on the edges of the complex. In words, the above conditions mean that there can be no |F|-bandlimited signals that are perfectly localized on the complementary setS. Perfect recovery of the signal s 1 from s 1 S can be achieved as where Q S = (I −D S F F ) −1 . The existence of the above inverse is ensured by condition (45). In fact, the reconstruction error can be written as [43] where we exploited in the second equality the bandlimited condition s 1 = F F s 1 .
To make (45) holds true, we must guarantee that D S F F s 1 = 0 or, equivalently, that the matrix D S F F is full column rank, i.e. rank(D S F F ) = |F|. Then, a necessary condition to ensure this holds is |S| ≥ |F |. An alternative way to retrieve the overall signal s 1 from its samples can be obtained as follows. If (45) holds true, the entire signal s 1 can be recovered from s 1 S as follows where U F is the E × |F| matrix whose columns are the eigenvectors of L 1 associated with the signal bandwidth |F|. Remark: It is worth to notice that, because of the Hodge decomposition (19), an edge signal always contains three components that are typically band-limited, as they reside on a subspace of dimension smaller than E. This means that, if one knows a priori, that the edge signal contains only one component, e.g. solenoidal, irrotational, or harmonic, then it is possible to observe the edge signal and to recover the desired component over all the edges, under the conditions established by Theorem 2.

B. Multi-layer sampling
In this section, we consider the case where we take samples of signals of different order and we propose two alternative strategies to retrieve an edge signal s 1 from these samples. The first approach aims at recovering s 1 sol and s 1 irr by using both, the vertex signal samples s 0 A = D A s 0 , with A ⊆ V, and the edge samples s 1 S = D S s 1 . Hereinafter, we denote by F irr , F sol and F H the set of frequency indexes in F corresponding to the eigenvectors of L 1 belonging, respectively, to the irrotational, solenoidal and harmonic subspaces. Note that, if s 1 is |F|-bandlimited then also s 1 H , s 1 sol and s 1 irr are bandlimited with bandwidth, respectively, |F H |, |F sol | and |F irr |. Furthermore, given the matrix U 0 with columns the eigenvectors of L 0 , we define the operator F 0 F0 = U 0 Σ F0 U T 0 where |F 0 | denotes the bandwidth of s 0 . Then, we can state the following theorem.
|F sH | can be recovered by using N 0 ≥ |F 0 | samples from s 0 and N 1 ≥ |F sH | samples from s 1 . Proof. Please see the supporting material document.
Let us now consider the case where we want to recover s 1 from samples collected from signals of order 0, 1, and 2, i.e. s 0 , s 1 and s 2 . We denote by s 2 M = D M s 2 the vector of triangle signal samples with M ⊆ T . Furthermore, given the matrix U 2 with columns the eigenvectors of the secondorder Laplacian L 2 , we define the operator where |F 2 | denotes the bandwidth of s 2 . Then, assuming s 2 bandlimited, it holds s 2 = F 2 F2 s 2 . Theorem 4: Consider the second-order simplex X = {V, E, T }, the edge signal s 1 = B 2 s 2 + s 1 H + B T 1 s 0 and assume that: i) the vertex-signal s 0 , the edge signal s 1 and the triangle signal s 2 are bandlimited with bandwidth, respectively, |F 0 |, |F 1 | = |F 0 | + |F H | + |F 2 | − (c 1 + c 2 ) and |F 2 |, with c 1 ≥ 0 and c 2 ≥ 0 the number of eigenvectors in the bandwidth of s 0 and s 2 belonging, respectively, to ker(L 0 ) and ker(L 2 ); ii) all conditions D A F 0 b) any |F|-bandlimited edge signal with |F| = |F 0 |+|F H |+ |F 2 | − (c 1 + c 2 ) can be recovered by using N 0 ≥ |F 0 | samples from s 0 , N 1 ≥ |F H | samples from s 1 and N 2 ≥ |F 2 | samples from s 2 .
Proof. Please see the supporting material document.

VI. INFERENCE OF SIMPLICIAL COMPLEX TOPOLOGY FROM DATA
The inference of the graph topology from (node) signals is a problem that has received significant attention, as shown in the excellent recent tutorial papers [33], [34], [35] and in the references therein. In this section, we propose algorithms to infer the structure of a simplicial complex. Given the layer structure of a simplicial complex, we propose a hierarchical approach that infers the structure of one layer, assuming knowledge of the lower order layers. For simplicity, we focus on the inference of a complex of order 2 from the observation of a set of M edge (flow) signals X 1 := x 1 (1), . . . , x 1 (M ) , assuming that the topology of the underlying graph is given (or it has been estimated). So, we start from the knowledge of L 0 , which implies, after selection of an orientation, knowledge of B 1 . Since L 1 = B T 1 B 1 + B 2 B T 2 , then we need to estimate B 2 . Before doing that, we check, from the data, if the term B 2 B T 2 is really needed. Since, from (23), the only components that may depend on B 2 are the solenoidal and harmonic components, we first project the observed flow signal onto the space orthogonal to the space spanned by the irrotational component, by computing where U irr is the matrix whose columns are the eigenvectors associated with the nonnull eigenvalues of B T 1 B 1 . Then, denoting with X 1 sH = (x 1 sH (1), . . . , x 1 sH (M )) the signal matrix of size E × M , we measure the energy of X 1 sH by taking its norm X 1 sH F : If the norm is smaller than a threshold η, we stop, otherwise we proceed to estimate B 2 .
The first step in the estimation of B 2 starts from the detection of all cliques of three elements present in the graph. Their number is T = trace (L 0 − diag(L 0 1)) 3 /6. For each clique, we choose, arbitrarily, an orientation for the potential triangle filling it. The matrix B 2 can then be written as where b n is the vector of size E associated with the n-th clique, whose entries are all zero except the three entries associated with the three edges of the n-th clique. Those entries assume the value 1 or −1, depending on the orientation of the triangle associated with the n-th clique. The coefficients t n in (53) are equal to one, if there is a (filled) triangle on the n-th clique, or zero otherwise. The goal of our inference algorithm is then to decide, starting from the data, which entries of t := (t 1 , . . . , t T ) are equal to one or zero. Our strategy is to make the association that enforces a small total variation of the observed flow signal on the inferred complex, using (31) as a measure of total variation on flow signals. We propose two alternative algorithms: The first method infers the structure of B 2 by minimizing the total variation of the observed data; the second method performs first a Principal Component Analysis (PCA) and then looks for the matrix B 2 and the coefficients of the expansion over the principal components that minimize the total variation plus a penalty on the model fitting error.
Minimum Total Variation (MTV) Algorithm: The goal of this algorithm is to minimize the total variation over the observed data set, assuming knowledge of the number of triangles. The set of coefficients t is found as solution of where t * is the number of triangles that we aim to detect. In practice, this number is not known, so it has to be found through cross-validation. Even though problem P MTV is non-convex, it can be solved in closed form. Introducing the nonnegative coefficients , the solution can in fact be obtained by sorting the coefficients c n in increasing order and then selecting the triangles associated with the indices of the t * lowest coefficients c n . Note that the proposed strategy infers the presence of triangles along the cliques having the minimum curl along its edges. Hence, we expect better performance when the edge signal contains only the harmonic components, whose curls along the filled triangles is exactly null.

PCA-based Best Fitting with Minimum Total Variation (PCA-BFMTV):
To robustify the MTV algorithm in the case where the edge signal contains also a solenoidal component and is possibly corrupted by noise, we propose now the PCA-BFMTV algorithm that infers the structure of B 2 and the edge signal that best fits the observed data set X 1 , while at the same time exhibiting a small total variation over the inferred topology. The method starts performing a principal component analysis of the observed data by extracting the eigenvectors associated with the largest eigenvalues of the covariance matrix estimated from the observed data set. More specifically, the proposed strategy is composed of two steps: 1) estimate the covariance matrix C X from the edge signal data set X 1 sH and builds the matrix U sH whose columns are the eigenvectors associated with the largest eigenvalues of C X ; 2) model the observed data set as X 1 sH = U sH S 1 sH and searches for the coefficient matrix S 1 sH and the vector t that solve the following problem , and setting to 1 the entries of t[k] corresponding to the t * smallest coefficients, and 0 otherwise; is a non-negative coefficient controlling the trade-off between the data fitting error and the signal smoothness. Although problem P T S is non-convex, it can be solved using an iterative alternating optimization algorithm returning successive estimates of S 1 sH , having fixed t, and alternately t, given S 1 sH . Interestingly, each step of the alternating optimization problem admits a closed form solution. More specifically, at each iteration k, the coefficient matrix S 1 sH [k] can be found as Then, given S

VII. NUMERICAL RESULTS
In this section, we test the validity of our inference algorithms over both simulated and real data.
Performance on synthetic data: Some of the most critical parameters affecting the goodness of the proposed algorithms are the dimension of the subspaces associated with the solenoidal and harmonic components of the signal and the number of filled triangles in the complex. In fact, in both MTV and PCA-BFMTV a key aspect is the detection of triangles as the cliques where the associated curl is minimum. Hence, if the signal contains only the harmonic component and there is no noise, the triangles can be identified with no error, because the harmonic component is null over the filled triangles. However, when there is a solenoidal component or noise, there might be decision errors. To test the inference capabilities of the proposed methods, as a first example we consider the simplicial complex illustrated in Fig. 2(a), composed of N = 274 nodes and with a percentage of filled triangles equal to one third. The observed edge flow is the result of the projection in (52) and so it contains only the solenoidal and harmonic components.
In particular, the solenoidal flow was generated as a linear combination of the eigenvectors associated with the 9 smallest nonnull eigenvalues of the true B 2 B T 2 , whereas the harmonic part was a linear combination of the first 9 eigenvectors in the kernel of L 1 . The complex recovered using the MTV algorithm is reported in Fig. 2(b), assuming t * = 4471, equal to the correct number of triangles. It can be noticed that, although the percentage of filled triangles is low and then there is a non negligible solenoidal component, the recovered topology is quite similar to the true simplex.
To better investigate this aspect, in Fig. 3(a), we report the triangle error probability P e , defined as the percentage of incorrectly estimated triangles with respect to the number of cliques with three edges in the simplex, versus the signalto-noise ratio (SNR), when the observation contains only harmonic flows plus noise. We considered a simplex composed of N = 50 nodes and with a percentage of filled triangles with respect to the number of second order cliques in the graph equal to 50%. We also set M = 50, t * = 105 and averaged our numerical results over 10 3 zero-mean signal and noise random realizations. The harmonic signal bandwidth |F H | is chosen equal to 105, which is equal to the dimension of the kernel of L 1 . From Fig. 3(a), we can notice, as expected, that in the noiseless case the error probability is zero, since observing only harmonic flows enables perfect recovery of the matrix B 2 . In the presence of noise, the MTV algorithm suffers and in fact we observe a non negligible error probability at low SNR. However, applying the PCA-BFMTV algorithm enables a significant recovery of performance, as evidenced by the blue curve that is entirely superimposed to the red curve, at least for the SNR values shown in the figure. In this example, the covariance matrix was estimated over 10 5 independent observations of the edge signals. The optimal γ coefficient was chosen after a cross validation operation following a line search approach aimed to minimize the error probability. The improvement of the PCA-BFMTV method with respect to the MTV method is due to the denoising made possible by the projection of the observed signal onto the space spanned by the largest eigenvectors of the estimated covariance matrix.
To test the proposed methods in the case where the observed signal contains both the solenoidal and harmonic components, in Fig. 3(b) we report P e versus the SNR, for different values of the dimension of the subspace associated with the solenoidal part, indicated as |F sol |. From Fig. 3(b), we can observe that the performance of both algorithms MTV and PCA-BFMTV suffers when the bandwidth |F sol | of the solenoidal component is large, whereas the performance degradation becomes negligible when |F sol | is small. In all cases, PCA-BFMTV significantly outperforms the MTV algorithm, especially at low SNR values, because of its superior noise attenuation capabilities.
Performance on real data: The real data set we used to test our algorithms is the set of mobile phone calls collected in the city of Milan, Italy, by Telecom Italia, in the context of the Telecom Big Data Challenge [45]. The data are associated with a regular two-dimensional grid, composed of 100×100 points, superimposed to the city. Every point in the grid represents a square, of size 235 meters. In particular, the data set collects the number N ij of calls from node i to area j, as a function of time. There is an edge between nodes i and j only if there is a non null traffic between those points. The traffic has been aggregated in time, over time intervals of one hour. We define the flow signal over edge (i, j) as Φ ij = N ij − N ji . We map all the values of matrix Φ into a vector of flow signals x 1 . We observed the calls daily traffic during the month of December 2013. The data are aggregated for each day over an interval of one hour.
Our first objective is to show whether there is an advantage in associating to the observed data set X 1 a complex of order 2, i.e. a set of triangles, or it is sufficient to use a purely graph-based approach. In both cases, we rely on the same graph structure, whose B 1 comes from the data set, after an arbitrary choice of the edges' orientation. If we use a graphbased approach, we can build a basis of the observed flow signals using the eigenvectors of the so called edge Laplacian in [32], i.e. L low We call this basis U low 1 . As an alternative, our proposed approach is to build a basis using the eigenvectors of where B 2 is estimated from the data set X 1 using our MTV algorithm. We call this basis U 1 . To test the relative benefits of using U 1 as opposed to U low 1 , we run a basis pursuit algorithm with the goal of finding a good trade-off between the sparsity of the representation and the fitting error. More specifically, for any given observed vector x 1 (m), we look for the sparse vector s 1 as solution of the following basis pursuit problem [46]: where V = U 1 in our case, while V = U low 1 in the graphbased approach. As a numerical result, in Fig. 4 we report the sparsity of the recovered edge signals versus the mean estimation error x 1 − Vs 1 F considering as signal dictionary V the eigenvectors of either the first-order Laplacian or the lower Laplacian. We used the MTV algorithm to infer the upper Laplacian matrix by setting the number t * of triangles that we may detect equal to 800. As can be observed from Fig. 4, using the set of the eigenvectors of L 1 yields a much smaller MSE, for a given sparsity or, conversely, a much more sparse representation, for a given MSE. An intuitive reason why our method performs so much better than a purely graphbased approach is that the matrix L 1 has a much reduced kernel space with respect to L low 1 and the basis built on L 1 captures much better some inner structure present in the data by inferring the structure of the additional term B 2 from the data itself.
As a further test, we tested the two basis U 1 and U low 1 in terms of the capability to recover the entire flow signal from a subset of samples. To this end, we exploit the band-limited property enforced by the sparse representation, enabling the use of the theory developed in Section V.A. Starting with the representation of each input vector x 1 as x 1 = Vs 1 , with either V = U 1 in our case, or V = U low 1 in the graph-based approach, we used the Max-Det greedy sampling strategy in [43] to select the subset of edges where to sample the flow signal and then we used the recovery rule in (47) to retrieve the overall flow signal from the samples. The numerical results are reported in Fig. 5, representing the normalized recovering error of the edge signal versus the number N s of samples used to reconstruct the overall signal. We can notice how introducing the term B 2 B T 2 , we can achieve a much smaller error, for the same number of samples.

VIII. CONCLUSION
In this paper we have presented an algebraic framework to analyze signals residing over a simplicial complex, with special focus on flow (edge) signals. We proposed methods to analyze signals of various order and to infer the structure of the simplicial complex from the data. We proved that, in applications over real traffic data, the proposed approach can significantly outperform methods based only on graph representations. Further developments include both theoretical aspects, especially in statistical modeling over simplicial complexes, and the application to a vast number of cases, as for example to biological or co-authorship data-sets.

APPENDIX A PROOF OF THEOREM 1
We begin briefly reviewing the basic properties of submodular functions and their Lovász extension [42], [47]. We recall its definition hereafter.
Definition 1: Let G : 2 E → R be a set function with G(∅) = 0. Let x ∈ R E be ordered w.l.o.g. in increasing order such that x j > x i } for i > 0. Then, the Lovász extension f : R E → R of G, evaluated at x, is given by [42]: Note that f (x) is piecewise affine w.r.t. x, and G(S) = f (1 S ) for all S ⊆ E. An interesting class of set functions is given by the submodular set functions, whose definition is: Definition 2: A set function G : 2 E → R is submodular if and only if, ∀A, B ⊆ E, it satisfies the following inequality: A fundamental property is that a set function G is submodular iff its Lovász extension f (x) is a convex function [42, p.172]. This implies that when f (x) is minimized on [0, 1] E , we find out the set of minimizers of the set function G.
Let us now prove Theorem 1. Assume the faces of X to be positively oriented with the order of their vertices, i.e. for the k-simplex σ k i we have as the set of σ 2 ∈ T with exactly one vertex in each set A i , for i = 0, 1, 2. We assume that ∀σ 2 i ∈ τ (A 0 , A 1 , A 2 ), there exists a permutation π of its vertices σ 0 such that and σ 0 ij ∈ A j , for j = 0, 1, 2. This returns a linear ordering of the vertices of the partition such that for all i < j, if σ 0 k ∈ A i and σ 0 m ∈ A j , then σ 0 k precedes σ 0 m in the vertices ordering. Define now, for a given partition {A i } i=0,1,2 , the set function g A0,A1,A2 : E → {0, 1} as where σ 1 i− denotes a 1-order simplex with opposite orientation with respect to σ 1 i . Given the 2-simplex σ 2 i = [σ 0 i0 , σ 0 i1 , σ 0 i2 ], we denote by σ 2 i \σ 0 i k for k = 0, 1, 2, the 1-simplex obtained by removing the vertex σ 0 i k from σ 2 i . More specifically, assuming as example k = 0, the edge ] is the boundary of σ 2 i after excluding the two edges impinging on node i 0 . Then, we can associate with the edge set function in (59) the following set function defined on triangles: (60) Since we assumed that ∀σ 2 i ∈ τ (A 0 , A 1 , A 2 ), there exists a permutation π of its vertices such that i k ∈ A k , for k = 0, 1, 2, we need to verify that the set function in (60) is invariant under a vertex permutation π. This can be proved by showing that the following equality holds where sgn(π) = (−1) Np , with N p the number of permutations made on the vertices such that σ 0 i k ∈ A k , for k = 0, 1, 2. Expanding the left side of (61) and using the fact that g A0,A1,A2 (σ 1 i ) = −g A0,A1,A2 (σ 1 i− ), the equality in (61) can be easily proved. Let us now observe that the set function G(A 0 , A 1 , A 2 ) in (60) is defined only on the set τ (A 0 , A 1 , A 2 ). However, if we consider all triangles in T , thanks to the linear ordering of the vertices of the partition, it holds 2 k=0 (−1) k g A0,A1,A2 (σ 2 i \ σ 0 i k ) = 0 for all triangles having two or three edges in the same set A k . This implies that: so that (60) can be rewritten as (62) Then, it holds where |τ (A 0 , A 1 , A 2 )| denotes the dimension of the set τ . Let us now derive the Lovász extension of this cut function. Denoting with σ 2 njk = [σ 1 n , σ 1 j , σ 1 k ] the oriented 2-order simplex with edges of indices n, j, k ∈ {1, . . . , E}, we can make explicit the dependence of the function G on the edge indexes n, j, k, writing (62) as wherẽ G njk = g A0,A1,A2 (σ 1 n ) − g A0,A1,A2 (σ 1 j ) + g A0,A1,A2 (σ 1 k ).
Then,G njk is a set function defined only on the power set of {n, j, k}. We can now derive its Lovász extension g(x n , x j , x k ). First, assume x n ≤ x j ≤ x k , with σ 2 njk = [σ 1 n , σ 1 j , σ 1 k ]. Then, from Def. 1, we have C 0 = {n, j, k}, C 1 = {j, k} and C 2 = {k}. Therefore, from (65) and using (59), it holds: Then, from (58) we get: Let us now assume x j ≤ x n ≤ x k . Then, we have C 0 = {n, j, k}, C 1 = {n, k} and C 2 = {k}, so that it results: and By following similar derivations, it is not difficult to show that for x n ≤ x k ≤ x j , x j ≤ x k ≤ x n , x k ≤ x j ≤ x n and x k ≤ x n ≤ x j it holds g(x n , x j , x k ) = x n − x j + x k . Therefore, from (64) defining the edge signal x 1 ∈ R E , we can write the Lovász extension of G as or, equivalently, as where B 2 (i, j) are the edge-triplet incidence coefficients defined in (3). This completes the proof of Theorem 1.

APPENDIX B SUPPORTING MATERIAL
This document contains some supporting materials complementing the paper: "Topological Signal Processing over Simplicial Complexes" submitted to IEEE Transactions on Signal Processing. Section A contains Proposition 2, and its proof, that will be instrumental in proving Theorems 3 and 4, respectively, in Sections B and C.

A. Proposition 2
We need first to derive the relationship between the bandwidth of the edge signal s 1 irr = B T 1 s 0 and that of the vertex signal s 0 , as stated in the following proposition.
Proposition 2: Let s 1 irr = B T 1 s 0 be the irrotational part of s 1 with s 0 = F 0 F0 s 0 a |F 0 |-bandlimited vertex signal. Then, s 1 irr is a |F 0 | − c 1 -bandlimited edge signal with c 1 ≥ 0 the number of eigenvectors in the bandwidth of s 0 belonging to ker(L 0 ). Proof. Let us define U 0 F0 the N × |F 0 | matrix whose columns u 0 i , ∀i ∈ F 0 are the eigenvectors of the 0-order Laplacian where the last equality follows from the bandlimitedness of , with C 1 the set of indices of F 0 associated with the eigenvectors belonging to ker(L 0 ), we get by stacking in the columns of U F0−C1 the eigenvectors u i associated with non-zero eigenvalues. Therefore, equation (71) reduces to From the equality B T 1 B 1 u i = λ i u i , multiplying both sides by B 1 , we easily derive u 0 i = B 1 u i so that equation (73) is equivalent to Since s 1 irr = B T 1 s 0 , we can rewrite (74) as This last equality proves that s 1 irr is a |F 0 | − c 1 -bandlimited edge signal with c 1 = |C 1 | ≥ 0. where with F sH the set of frequency indexes in F corresponding to the eigenvectors of L 1 belonging to the solenoidal and harmonic subspaces. Assuming that both the conditions D A F 0 F0 2 < 1 and D S F FsH 2 < 1 hold true, the matrix G is invertible and from the inverse of partitioned matrices [48], we get Then, using the bandlimitedness property, so that s 0 = F 0 F0 s 0 , s 1 H = F F H s 1 H and s 2 = F 2 F2 s 2 , we get  Under the assumptions D A F 0 F0 2 < 1, D S F F H 2 < 1 and D M F 2 F2 2 < 1, the matrixḠ becomes invertible and from the inverse of partitioned matrices [48], we get This concludes the proof of point a). Let us prove next point b). From Proposition 2, s 1 irr is a |F 0 |−c 1 -bandlimited edge signal. To find the bandwidth of the solenoidal part s 1 sol = B 2 s 2 we can proceed in a similar way to the proof of Prop. 2. Using the bandlimitedness property of s 2 , so that s 2 = F 2 F2 s 2 , it results with U 2 F2 ∈ R T ×|F2| the matrix whose columns u 2 i , ∀i ∈ F 2 are the eigenvectors of the second-order Laplacian L 2 = B T 2 B 2 . From Proposition 1 at each eigenvector u 2 i with u 2 i / ∈ kern(B T 2 B 2 ) corresponds an eigenvector u i = B 2 u 2 i of B 2 B T 2 with the same eigenvalue. Let us write U 2 F2 as where the columns of U 2 C2 are the c 2 = |C 2 | ≥ 0 eigenvectors in the kernel of L 2 belonging to the bandwidth of s 2 . Therefore, it results where we used the equality B 2 U 2 F2−C2 = U F2−C2 . Then equation (88) reduces to Since it holds B T 2 u i = u 2 i , ∀u 2 i / ∈ ker(B 2 ), we easily get Then s 1 sol is a |F 2 | − c 2 -bandlimited edge signal, so that s 1 is |F|-bandlimited with |F| = |F 0 | + |F H | + |F 2 | − (c 1 + c 2 ). Let us now consider the system in (87). We get Hence, from (93), it follows that to perfectly recover the harmonic part of s 1 we need a number of edge samples N 1 at least equal to the signal bandwidth |F H |. Finally, from the first and last equation in (92), we conclude that to retrieve s 0 and the solenoidal signal s 1 sol we need, respectively, a number of node samples N 0 ≥ |F 0 | and a number of triangle signal samples N 2 ≥ |F sol |.