{'fraction_non_alphanumeric': 0.030244201491774803, 'fraction_numerical': 0.005977316848881169, 'mean_word_length': 4.281975175391257, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Graph reconstruction from the observation of diffused signals

Signal processing on graphs has received a lot of attention in the recent years. A lot of techniques have arised, inspired by classical signal processing ones, to allow studying signals on any kind of graph. A common aspect of these technique is that they require a graph correctly modeling the studied support to explain the signals that are observed on it. However, in many cases, such a graph is unavailable or has no real physical existence. An example of this latter case is a set of sensors randomly thrown in a field which obviously observe related information. To study such signals, there is no intuitive choice for a support graph. In this document, we address the problem of inferring a graph structure from the observation of signals, under the assumption that they were issued of the diffusion of initially i.i.d. signals. To validate our approach, we design an experimental protocol, in which we diffuse signals on a known graph. Then, we forget the graph, and show that we are able to retrieve it very precisely from the only knowledge of the diffused signals.


I. INTRODUCTION
Signal processing on graphs has recently gained a lot of interest, and recent years have seen the appearance of numerous new techniques based on its principles. Although it was orignially seen as an extension of classical signal processing to more complex domains -thus generalizing the Fourier approach -, applications now emerge in various fields such as wavelets on graphs, that find application in fields such as brain imagery [12]; or filtering, applied to big data analysis [16]. Other examples are the recent work of Dong et. al [7] that have studied the political orientations of the Swiss cantons from voting data; or the application of Abry et. al [1] to photographic prints classification.
One common aspect of all the existing applications is that they rely on the Graph Fourier Transform (GFT) of the graph, which is a cornerstone of signal processing on graphs. This particular transform is performed using the eigenvectors of a particular matrix that is dependent on the graph used to represent the support of the signals. However, such a matrix is not always available. In some even more complex -and frequent -cases, the graph has no real physical existence. Example of this are when sensors are thrown in a field, or are used to measure the brain activity (EEG). There is an obvious link between the signal components observed in the various signals, but there is no intuitive notion of graph.
To circumvent this problem, methods have recently been developed to allow one to recover a graph structure from the observation of signals. The idea behind this is to provide a structure that would naturally represent the relationships between the signal components. Even in the case when the graph does not exist, such techniques should be able to provide a graph that captures these relationships, and would provide a support for the application of graph signal processing tools.
One possible approach was proposed in the Ph.D. dissertation of Costard [5]. Her work is based on previous work from Dempster [6], that estimates the covariance matrix of signals with the constraint that some values of its inverse -called the precision matrix -are null. This latter work has then allowed Wermuth [19] to show that the non-null entries of this precision matrix correspond to edges of a graph representing the studied signals. Based on this result and on the additional works in [15], [2] on model selection, Costard [5] reviews existing methods for recovering a sparse precision matrix. Additionally, she introduces a new method called ABIGlasso, based on Bayesian inference techniques initialized using the Graphical lasso [9]. Related practical algorithms for estimating a sparse inverse covariance matrix were developed since then (eg. [10]).
Another existing approach to retrieve a graph from the observation of signals is the one introduced in the recent work of Dong et. al [7]. The notion of smoothness of a signal on a graph was originally defined in [17] as a quantity representing the total variation of a signal with respect to adjacent nodes. Based on the assumption that adjacent nodes in a graph should share similar signal values -i.e signals should be smooth -Dong et. al propose a graph reconstruction technique that enforces this smoothness constraint.
To the best of our knowledge, these two approaches are the main efforts that have been made to answer the problem of graph reconstruction from the observation of signals. However, for the approach of Costard, the author shows that her work is only applicable to a small number of variables [5]. Concerning the approach of Dong et. al, they use the assumption that signals should be smooth on the recovered graph. In this paper, we propose a new method to recover a graph from signals, with different assumptions on the observed signals; that is to say they are issued from the diffusion of initially i.i.d. signals in the graph. Our method consists of serching a solution in a system of linear constraints, which is of polynomial complexity [3] and can be done efficiently via the interior point method [13].
In this article, we first introduce in Section II some definitions and properties that are needed for a complete understanding of our method. Then, in Section III, we introduce the setup that enables us to validate our approach, presented in Section IV. Finally, experimental evaluation is provided in Section V, and we discuss our method and future work in Section VI.

II. DEFINITIONS AND PROPERTIES
In the context of signal processing on graphs, a graph is a mathematical structure used to represent the support of the observed signals. Particular instances of graphs can be used to represent the classical domains. For example, a chain graph can model the support of time-defined signals, and a grid can model the supporting pixels of an image.
With this definition, a signal on a graph can be seen as a value attached to every node in the graph: A convenient way to represent a graph is through its adjacency matrix W, defined as follows: Definition 3 (Adjacency matrix): In this article, we consider graphs with the following properties: From the definition of the W, we can define the diagonal degrees matrix D of the graph, representing the degree of every node of V: Definition 4 (Degrees matrix): Using the definitions of W and D, we can finally define the normalized Laplacian Ł of the graph. This is an operator analogous to the classical Laplacian operator [17], that models the derivative of one diffusion in the graph: Definition 5 (Normalized Laplacian): is called the diffusion matrix. By construction, and with our assumptions on G, Ł is a real symmetric matrix, and can thus be diagonalized into a basis of eigenvalues Λ Ł = diag(λ 1 . . . λ N ) and associated eigenvectors X Ł = {χ 1 ; . . . ; χ N }. From [4], we know that the eigenvalues λ 1 ≤ λ 2 ≤ . . . ≤ λ N of Ł have some interesting properties. One of them is that they are all in [0; 2], with λ N = 2 for and only for bipartite graphs (i.e graphs for which there exist two sets V 1 and V 2 such that . As a consequence, the eigenvalues of T are in [-1; 1].

III. PROBLEM SETUP
In this paper, we propose a method to recover a graph G from the observation of signals diffused on it, with the assumption that they were initially i.i.d. Hereafter, we denote as X the N × M matrix of i.i.d. signals, and Y = T K X the matrix of diffused signals, for a given value of K. Thus, our problem becomes as follows: From the only observation of Y, and under the assumption that Y is issued from the diffusion of initially i.i.d. signals, is it possible to recover T, and thus G?
More specifically, to validate our algorithm in Section IV, we create an experimental protocol consisting of five steps: 1) We randomly generate a graph G, and compute the associated diffusion matrix T. In real applications, this graph is unknown, and is what we want to learn using our algorithm.
We use the algorithm presented in Section IV to recover an estimate T for T from only knowledge of Y.

5) We recover G from T.
To study the reconstruction error, we use the root mean square error (RMSE) between the real and recovered adjacency matrices: Definition 6 (Root mean square error): In this work, we only consider graphs such that the eigenvalues of T have pairwise distinct absolute values. Note that this restricts the class of considered graphs to nonbipartite graphs that have one eigenvalue equal to 1 and one equal to −1.

IV. GRAPH RECONSTRUCTION PROCESS
In this section we detail our algorithm to recover an estimate G for a graph G, from only knowledge of Y, as described in Section III.

A. Covariance of the diffused signals
Let Σ Y be the covariance matrix of Y. With the assumption that the signals Y are issued from the diffusion of initially i.i.d. signals, we can make the following development: with the last simplification coming from -with the i.i.d. assumption -E[XX ⊤ ] = I N and from the observation that T is symmetric. If we denote by Σ Y the empirical covariance of Y, it is well known that Σ Y converges to Σ Y as M tends to infinity. A consequence of this property is that we can obtain the eigenvectors of T from the diagonalization of Σ Y -since the eigenvectors of all powers of a matrix are identicalprovided that M is big enough.
However, this is not the case for the eigenvalues of T.
Since Σ Y is equal to T 2K , we have the property that the eigenvalues Λ T 2K of T 2K are all positive. Therefore, we can only obtain the absolute values of the eigenvalues of T -hereafter denoted Λ | T| -through the 2K-root of T 2K .
In the remainder of this paper, we will assume that Σ Y = Σ Y . The impact of the noise when retrieving the covariance matrix will be briefly analyzed in Section V, and will be the main part of our future work.

B. Recuperation of the signs of Λ T
At this point, we have a way to recover the eigenvectors X T of T and the absolute values Λ | T| of its eigenvalues, from only knowledge of the signals Y, provided that M is big enough. Now, we want to be able to recover a vector s of signs (−1 or +1), such that Λ T = diag(s)Λ | T| are eigenvalues such that T = X T Λ T X T ⊤ obeys some chosen properties: • Values on the main diagonal of T are all equal to 0. • Off-diagonal values are either positive or null. • T is symmetric. These three properties are by construction observed for any diffusion matrix in the class of graphs that we consider. Therefore, we will define a set of constraints that enforce these properties.
The problem in retrieving s directly is that the values we want s to take are non-continuous, making the problem too difficult to solve. As a consequence, we will not put any constraint on the domain of s, and will consider that it takes its values in R N .
Based on the observation that the trace of T should be equal to 0, and using the equality between the trace of a matrix and the sum of its eigenvalues, we create two sets of constraints: 1) N equality constraints : These constraints enforce the main diagonal of T to be null, when considering s as its eigenvalues. 2) N (N −1) 2 inequality constraints: with triu(N ) being a function returning the indices of the triangular upper part of a N × N matrix. These constraints enforce the off-diagonal elements to be either positive or null. Moreover, considering only the triangular upper part enforces the symmetry of the recovered T. Since a null vector of dimension N is a trivial solution to the constraints, we use the observation that the eigenvalue associated to the constant-sign eigenvector should be equal to 1 [17]. Therefore, we add a final constraint to prevent this case: with O N−1 being a null vector of size N − 1.
By constraining the cells of T through (6) and (7) -using its eigendecomposition T = X T diag(s)X T ⊤ -we make s converge toward a vector of values that should tend to the same sign as the expected eigenvalues. Finally, we only keep the signs of the eigenvalues, and inject them in the previously obtained Λ | T| : with ⊙ being the elementwise product of vectors, and sign(s) being a function returning a vector in which the i-th component is equal to −1 when s i < 0, and +1 otherwise. We can note that, even though the problem is linear, we have a quadratic number of constraints in 2). We will show in Section V that we can limit ourselves to only a fraction of these constraints, thus reducing the complexity while still obtaining the same results. A solution to this set of linear constraints can be efficiently obtained using the interior point method [13].
Remark 1: The restriction of our method to graphs having identical absolute eigenvalues, as described in Section III, can be understood here. As a matter of fact, if T has eigenvalues of opposite signs, then when retrieving Σ Y ∼ T 2K , these signs have disappeared. A consequence of this observation is that the diagonalization of Σ Y will not necessarily give the same eigenvectors as those of T, since two distinct eigenspaces were merged into one.

C. Summary of the algorithm
In this section, we have proposed an algorithm to retrieve T from the observation of diffused signals Y, under the assumption that they were initially i.i.d. Our method works in five steps, summarized as follows: 1) We compute Σ Y , the empirical covariance matrix of

5) Since
T is a simple normalization by the degrees of W, then it is null where and only where W is null. As a consequence, we fix a small ε -to circumvent machine precision limitations -and retrieve W from T as follows: ∀i, j ∈ {1; . . . ; N } : W i;j = 0 if T i;j < ε 1 otherwise .
(10) As we will show in Section V, this algorithm works well for the class of considered graphs. However, we make in this first work some assumptions on the number of observed signals, and on the distinguishability of the eigenvalues of the covariance matrix. Also, the value of K is currently fixed in our work, and is common for every signal. We think that we can get rid of this additional parameter, and will study this hypothesis in future work.

V. EXPERIMENTAL EVALUATION
In this section, we first study the reconstruction performance -and the associated execution time -of our method for various graph orders, using the settings introduced in Section III. Then, we briefly study the impact of the value of M -and thus of the noise on the covariance matrixon the successfulness of our algorithm.
To perform our tests, we make some experimentation choices, such as the use of an Erdős-Rényi model [8], with edge probability P = 0.3, for G; and the choice of a normal distribution N (0, 1) for X. Similar results are observable for different graph definitions -such as random geometric graphs -and choices of distribution.
All tests in this document are performed using the CVX package [11] with solver SDPT3 [18] for MATLAB [14], with default parameters.

A. Reconstruction performance of the method
We have seen in Section IV-B that the number of inequality constraints describing the matrix to reconstruct was quadratic in the dimension of the signals. For this reason, we limit ourselves a randomly selected subset of them. By doing so, we slightly increase the number of possible solutions.
However, it appears that this class of solution vectors has still a very low cardinality -in fact it is most of the time unique -and that it is sufficient to guide our solver by constraining only a small number number of off-diagonal cells to be positive. Figure 1 shows the mean RMSE between the initial W and the reconstructed estimate W, for various ratios of kept constraints, i.e for various values of α in α N (N −1) 2 .
Here, we use an oracle to find the best ε possible in (10), although the gap with the real values of T is generally high and it is quite easy to distinguish noise from data. We make this choice because we want to study the reconstruction performance of our method. In real applications, ε should impose a degree of sparsity on the reconstructed matrix.
We can see from this table that our algorithm works quite well, and succeeds in reconstructing the original graph in 100% cases in our tests when we keep a number of inequality constraints over 1 N . We can also see that even with a smaller number of constraints, the method is still robust, but some errors start to appear, mainly for smaller graphs. However, the first column shows that keeping some of these inequality constraints are essential to guide the solver toward the correct solution.
Because removing some constraints does not overly increase the error, we study the impact it has on the execution time of our algorithm. Figure 2 depicts the time required by our algorithm to reconstruct the graph, as a function of the graph order, for the same ratios of kept constraints. As expected, removing some inequality constraints in (7) considerably reduces the execution time of our algorithm. As a matter of fact, we can remove a nonlinear number of constraints while keeping the same performance, thus making the number of inequality constraints linear. A consequence of this is that our method is now based on finding a solution to a linear number of linear constraints. Therefore, it can be applied to graph reconstruction with high dimension signals.

B. Impact of the noise on the reconstruction
At this point, we have made one strong assumption that the number M of observed signals is big enough to consider Σ Y = Σ Y . Although it is true in the limit case of an infinite number of signals, in real applications there is some noise in the empirical covariance matrix. Figure 3 shows the impact of the value of M on the RMSE, for randomly generated graphs of 15 nodes.
In addition to these curves, it is interesting to remark that the RMSE when performing our tests is generally either 0 or very close, or around 0.5. Although the mean error decreases as the M increases, the previous observation shows that noise RMSE as a function of the number of signals, for graphs of 15 nodes, when using Σ Y instead of Σ Y , with a ratio of inequality constaints α = 1 N . Mean for 1000 randomly generated graphs.
when retrieving Σ Y is the main problem of our method, and will therefore be the main focus of our future work. One promising direction comes from the observation that most of the error is concentrated in the eigenvectors of Σ Y and not in its eigenvalues. Therefore, we might be able to use the fact that the first eigenvector of T should be of a certain form -from the observation that the components of the first eigenvector of Ł are values that are function of the graph degrees -to realign the eigenvectors of Σ Y properly.

VI. CONCLUSIONS
In this paper, we have introduced a scalable method to retrieve a graph from the observation of signals diffused on it. We have experimentally shown that we could reduce the number of inequality constraints in our problem, thus making our method as complex as finding a solution to a linear number of linear constraints.
The main limitation of our work is currently the impact of the noise on the reconstruction of the covariance matrix of the observed signals. For this reason, this will be the core aspect of our future studies. Additionally, we want to get rid of the parameter K, describing the number of diffusions of the observed signals, in our algorithm. This would enable us to consider signals observed at different points in time in an unified framework.