Bayesian Topology Learning and noise removal from network data

Learning the topology of a graph from available data is of great interest in many emerging applications. Some examples are social networks, internet of things networks (intelligent IoT and industrial IoT), biological connection networks, sensor networks and traffic network patterns. In this paper, a graph topology inference approach is proposed to learn the underlying graph structure from a given set of noisy multi-variate observations, which are modeled as graph signals generated from a Gaussian Markov Random Field (GMRF) process. A factor analysis model is applied to represent the graph signals in a latent space where the basis is related to the underlying graph structure. An optimal graph filter is also developed to recover the graph signals from noisy observations. In the final step, an optimization problem is proposed to learn the underlying graph topology from the recovered signals. Moreover, a fast algorithm employing the proximal point method has been proposed to solve the problem efficiently. Experimental results employing both synthetic and real data show the effectiveness of the proposed method in recovering the signals and inferring the underlying graph.


Introduction
Many applications including internet of things (IoT) networks, sensor networks, social networks, gene networks, customer consumption data in utility companies, financial data, regional temperature data, and brain computer interface measurements result in rapidly growing data volumes. IoT networks are a rapidly emerging technology where smart sensors are connected to easily enable data collection anywhere and anytime [1]. Figure 1 illustrates how an IoT network may involve collecting large volume of multi-variate time series data which are connected in some abstract senses.
Storing and analyzing these huge data volumes is challenging due to the need for high computation and power resources. The emerging field of graph signal processing (GSP) [2,3] simplifies the analysis of large data volumes by the use of graph theory, where each graph node represents a component of the system. Several applications can be analyzed by graphs, which can capture the underlying topology among different entities of the network. For example, a network of thermal sensors which measures temperatures in different regions, can be modeled by a graph.
In some cases, the graph topology is not known a priori. Thus, one of the desired goals of the GSP framework is estimating the underlying topology given a set of measured data. Some researches concentrated on directed topology estimation [4][5][6][7][8][9][10][11][12][13][14], where specific process models have been used which are suitable for limited applications. However, 1 3 the Gaussian model is a ubiquitous signal model which can be utilized for a wide range of applications. By using data sets generated by Gaussian distributions, the solution of a topology inference problem usually leads to an undirected graph structure. The "covariance selection" proposed by Dempster [15], was one of the pioneering works to capture the connectivity underling the Gaussian measurements. Later on, Banerjee et. al. [16] proposed an ℓ 1 regularized optimization problem to find a sparse precision matrix which is the inverse of the covariance matrix. The precision matrix carries information about the partial correlations between the random variables. Friedman and his colleagues have proposed a fast algorithm, called graphical Lasso, to implement inverse covariance estimation [17]. There are other works investigating the inverse covariance matrix estimation with different approaches and different rates of convergence [18,19]. However, none of these approaches can leverage the signal model to infer the underlying graph topology. In other words, they focus on pairwise relationships of entities to find their structure instead of applying the signal model information.
Besides, these approaches are very sensitive to noise since they are based purely on the measurement values.
A state-of-the-art machine learning approach for estimating a graph topology from measurements generated by a Gaussian Markov Random Field processes has been proposed in [20]. The main limitations of this algorithm are twofold. First, no low cost implementation is proposed that scales efficiently with the number of nodes. Second, the effect of noisy measurements on the performance evaluation has not been considered. Kalofolias has proposed an efficient and scalable algorithm to learn a graph from noise-free observations where it is assumed that the node degrees are positive (they can not be zero) [21]. A generalization of the methods in [20] and [21] for different Laplacian and structural constraints has been proposed in [22]. Using graph signal frequency domain analysis, the underlying network topology is inferred from spectral templates in [23]. The authors assumed that there exists a diffusion process in the graph shift operator (GSO 1 ) that can generate the given observations and the GSO is jointly diagonalizable with the observations covariance matrix. This method was first applied to capture the underlying topology when the graph signals are generated from a GSO filter output fed by a white signal input. Then, this method was extended to the colored input, generating nonstationary diffusion processes 2 . Some other research works with similar ideas are presented for different stationary and non-stationary processes in [24][25][26][27][28]. Dictionary learning [29,30] and transform learning [31] have also been used for inferring the graph topology. In [29][30][31], a specific relation between the Laplacian matrix and the dictionary atoms has been sought and hence these algorithms are applicable when we have some knowledge about signal representation in Fig. 1 Illustration of learning a network graph structure from IoT data 1 The GSO is a matrix which captures the graph's local topology and the graph Fourier transform is defined using its eigenvectors. 2 A process is wide sense stationary if its first moment is constant over the graph vertex set and its covariance matrix is invariant with respect to the localization operator [24].
The inexpensive sensors found in a wide variety of applications distort the observed data samples. The resulting measurement errors in the data analysis are conventionally modeled as noise [42] and it is generally desirable to remove this noise in many applications [43,44], which generally produced a more accurate learned topology. In a recent research [45], Liu and her colleagues proposed a noisy data cleansing algorithm to remove unwanted distortion from Industrial IoT (IIoT) sensor data in manufacturing. They showed the efficiency of their proposed method on the multi-variate timeseries data collected from a real IIoT based four stage compressor manufacturing plant. In [46], the effect of distortion on the trading data has been discussed, where the stock market's impact on this data is characterized by a graph [47]. There are also other studies in the GSP framework considering the noise effect on the observations, e.g. [48], but they assumed that they have knowledge about the underlying graph and did not tackle graph topology inference and noise removal at the same time. While it is a naive strategy to find the graph topology first and then remove the noise from the graph signal, this approach has a low performance due to the accumulation of errors in topology learning and signal estimation.
To solve this type of problem in a topology learning framework, an optimization problem is usually formulated with at least two objective terms. The first term is always a data fidelity term controlling the energy of the signal recovery error. The second term takes care of the graph signal smoothness over the learned topology. A smooth signal has smooth variations on the underlying topology, or in other words, connected nodes in the estimated structure have similar signal values. To overcome over-fitting, the data fidelity term is regularized by a smoothing term and an exhaustive search is applied to find the regularization parameter. In many applications, the relationship between the noise variance and the regularization parameters may be held directly, e.g. image restoration [49]. Thus, estimating the noise variance not only provides a cleaner version of the signal for further analysis, but also helps to overcome the overfitting problem. All things considered, the estimation of noise variance is the main step to design the filter for smoothing the signal and denoising the given measurements.
It is usually assumed that the measurement noise has a zero mean Gaussaian distribution model which is independent of the desired signal. Thus, finding the noise variance is of interest to design a filter for signal denoising. The noise variance estimation is not only useful from a signal recovery perspective, it is also important from a graph learning point of view. Since finding the underlying structure from a noise-free data set is more accurate than from noisy graph signals, signal denoising can contribute to estimating a more accurate graph topology. In this respect, Chepuri et al. [50] proposed an approach to find the connectivity graph and remove noise from measurements. They proposed a non-convex, cardinality constrained, boolean optimization problem along with a convex relaxation. They have assumed that the number of edges is known in advance and the proposed method scales with the desired number of edges. Another approach has been proposed in [51], which estimates the topology and removes the noise from the measured signals, simultaneously. They adopted a factor analysis model for the multi-variate signals and imposed a Gaussian prior on the latent variables that control the graph signals. By applying a Bayesian framework, the maximum a posteriori (MAP) probability estimate of the latent variable is investigated. Considering signal smoothness over the underlying graph, this procedure leads to an optimization problem over the graph topology and graph signals, simultaneously. A limitation of this approach is that an exhaustive search is applied to find the regularization parameter or noise variance which is not very suitable, especially for the situation in which the measurement noise varies over time.
Our contribution Given a set of noisy multi-variate signal measurements, we propose a new algorithm to perform the underlying topology learning and graph signal noise removal. This graph topology characterizes the affinity relationship among the multi-variate signals. We have proposed a minimum mean square error estimation (MMSE) approach for signal recovery and topology learning and unlike [51], we estimate the regularization parameters analytically rather than by using a grid search. We show that the regularization parameters are linked to the noise variance which is usually unknown, and hence, finding the noise variance helps adjusting the parameters precisely. Besides, it is the main role in designing a filter to denoise the graph measurements. In other words, to compare the current paper with other similar ones, especially [51], two specific points must be considered; (1) understanding the relation between amount of noise and the regularization parameters in the denoising procedure, and (2) a rigorous strategy to estimate the noise variance instead of assuming that it is known in advance. The MMSE procedure is utilized to estimate the optimal Laplacian matrix, whose eigenvalue matrix is the precision matrix of the Gaussian Markov Random Field process. This work is an extension of our previous paper [52], in which we studied the problem of joint graph signal recovery and topology learning using an off-the-shelf optimization toolbox to implement the algorithm. However, in this version, we propose a fast algorithm and prove its convergence analytically. Moreover, we provide simulation results on real world data sets in this paper. To compare the results with the state of the art methods, we applied some measures which evaluate the performances from two distinct perspectives which are signal recovery error and graph topology estimation accuracy. Therefore, we can compare our proposed method against other methods that employ signal denoising strategies, like the one proposed in [22]. A simple graph filter is applied based on the estimated noise variance which can be used in many multi-variate measurements applications, including industrial internet of things IIoT applications. In the preprocessing stage in some IIoT applications, it is necessary to denoise the measurements, remove outlier data, detect anomalies, and estimate the missing values of some data. By the proposed method, the knowledge of the underlying data structure can be applied to implement all of these tasks. As an application of our idea in the area of IIoT networks, we apply our proposed method to estimated the missing values of sensor readings for a power consumption dataset.
The rest of the paper is organized as follows; in Sect. 2, some preliminaries about graph theory and signal processing on graphs are presented. Section 3 formulates the graph topology learning problem via a Bayesian framework and proposes a general algorithm for implementation, called Bayesian Topology Learning (BTL). In Sect. 4, a method is proposed to implement BTL efficiently via solving a proximal point algorithm, called BTL-PPA. The proof of convergence is also presented in Sect. 5. Section 6 presents experimental results obtained from the proposed method for temperature sensor networks and IIoT applications by using both real and synthetically simulated data. Finally, the Sect. 7 concludes the paper. Throughout the paper, lowercase regular letters, lowercase bold letters, and uppercase bold letters denote scalars, vectors, and matrices, respectively. The rest of the notation is presented in Table 1.

Background on graph signal processing
Let G = (V, E, ) to be a graph with the vertices v i ∈ V , the edge set E , and the undirected weight matrix . Each edge (v i , v j ), 1 ≤ i, j ≤ N has a weight of w ij in the corresponding entry of ∈ ℝ N×N + . Thus, w ij = 0 implies no connection and w ij = w ji > 0 quantifies the strength of connection between the node v i and the node v j . For simplicity of presentation, we assume that the diagonal entries of the weight matrix are zero, i.e. there is no self loop. The adjacency matrix  stores only the existence of the edges, regardless of their weights. In other words, if there is a connection or edge between the vertices v i and v j , the adjacency matirx entity a ij = 1 and zero otherwise. The degree matrix is defined as The degree d i is the sum of the edge weights connecting node i to its neighbors. Then, = diag( ⋅ N ) , where N is the all ones N × 1 vector. The combinatorial and normalized graph Laplacian matrices are also defined as = − and norm = N − − 1 2 − 1 2 , respectively, where N is the identity matrix of size N × N . Since the Laplacian matrix is real and symmetric, its eigendecomposition is written as follows where and ∈ ℝ N×N are eigenvalue and eigenvector matrices, respectively, and (⋅) T denotes the matrix transpose operator. For the normalized Graph Laplacian, all eigenvalues are between 0 and 2. Since the Laplacian matrix can uniquely characterize the graph structure, the graph topology learning problem can be formulated as a problem of graph Laplacian matrix learning. The k'th graph signal is represented as below An important concept in the GSP framework is signal smoothness with respect to the intrinsic structure of the graph. The local smoothness around vertex i at time k is given as [2] where N i is the neighborhood of node i. Then, to quantify the global smoothness (GS), we have [2] Figure 2 compares the smoothness of a signal living in different topologies.

Bayesian Topology Learning
Bayes' rule relates the probability of an event to the prior knowledge that we have around it and updates the belief. Here, in the same way, it is assumed that we have some knowledge about the latent variable which factorizes the graph signals and also relates to the underlying graph structure. By Bayesian inference, we investigate the posterior probability of this latent variable and use its probability distribution function to find the graph Laplacian matrix and denoise the observed graph signals. The signal is smooth with respect to the G 1 , less smooth with respect to G 2 and likely to be non-smooth with respect to G 3 . The edge weights are all ones. These figures are generated by GSPBOX [53] Assume that we are given a data matrix of size N × K , containing noisy graph signal measurements in its columns. The rows of the matrix correspond to the vertices of the underlying graph and each column/measurement is as follows Let us adopt a multi-variate Gaussian distribution for the measurement noise [k] , given by the following probability density function (pdf ) where ∈ ℝ N is all zero vector. 3 To find the underlying graph topology, first, we use the factor analysis model to leverage a representation matrix that can be linked to the graph Laplacian/topology directly. In this model, each graph signal is represented by , ∀k , where the unobserved latent variable [k] ∈ ℝ N controls each graph signal via the eigenvector matrix [54] and [k] ∈ ℝ N is the mean of graph signal [k] . For simplicity and without loss of generality, hereafter, it is assumed that ∀k, [k] = . As discussed in [51], the motivation of using the factor analysis model is that each graph signal can contribute to the underlying graph structure estimation. Similar to [51], it is assumed that  (12) is simplified as follows where † = † T . The posterior density of is obtained by applying the Bayes rule as follows and the Minimum Mean Square Error (MMSE) estimator of the latent variable is defined as below [55] where the covariance matrix is If increases, the distributions in (11) and (12) have larger uncertainties and thus the estimation of (16) is less accurate and has a larger MSE. Due to the linear relationship between and , i.e. = , the graph signals can be estimated by ̂ = , which can be simplified as follows where and will be estimated later. Equivalently, (18) can be represented in a matrix form as ̂ = ( + ) −1 which is similar to the one presented in [51], where ̂ = ( + ) −1 for a parameter . In [51], is the regularization parameter corresponding to the smoothness of graph signals over the topology. In other words, the graph signals are smoothed by the graph filter ( + ) −1 and hence = can be called the filter coefficient or the regularization parameter or noise variance, interchangeably. If approaches zero, ̂ and are going to be identical and for larger , the effect of graph Laplacian matrix on filtering the signal is larger. In other words, we have and thus for the true denoised version of , i.e. ̂ = , we have = ̂ . Given a fixed error , if is large, then ̂ is sufficiently smooth hence ̂ will likely be small. This is the case where we make a big effort to denoise the observation . If is small, the opposite will be true. So this shows that given a noisy observation (hence the amount of noise added to the clean signal is "fixed"), different estimation approaches (or equivalently adjustment in [51]) will lead to different denoising effect. However, has not been investigated analytically in [51] and it was estimated by a grid search.
To make the importance of noise variance estimation more clear, we generate a graph with N = 50 nodes with K = 100 graph signals, i.e. , via the procedure explained in the simulation Sect. 6.1. Then, the graph signals are contaminated by Gaussian noises with difference variances to provide . Given and the known Laplacian matrix , we try to denoise the graph signals via (18) by two different filters, i.e. employing the correct variance and an incorrect estimated variance which is half of the true value. This procedure is repeated for 10 different variances and the Fig. 3 The effect of denoising the graph signals when using a correct noise variance estimator normalized mean square error of the true signals and the estimated ones are computed. The results in Fig. 3 show how we can improve our measurements by the filter in (18), especially when we have a better estimation of the noise variance. Since the Laplacian matrix estimation is performed in the next step based on these measurements, using denoised graph signals certainly helps in a better topology estimation.
To estimate the parameters and † , an expectation maximization procedure is used which proceeds by optimizing the following log-likelihood function where (⋅) denotes the expectation and ̂ † , ̂ , and ̂ are the estimators of † , , and . Hereafter, for brevity of notation, denotes all the parameters, i.e. ∶= , , † = ( , ) . By using (10)- (12) and (17) and some manipulations, (20) is rewritten as follows (see Appendix 1 for more details) where = 1 K̂ ̂ T is the empirical covariance matrix of the graph signals and the pseudo-determinant | ⋅ | is used due to the singularity of and 0 . The pseudo-determinant of a matrix is the product of its non-zero eigenvalues. To find the optimal parameters, Q must be optimized with respect to each argument iteratively up to a convergence. By taking the derivative with respect to , we have where the optimal is found by a numerical method, e.g. Newton-Raphson algorithm. Then, Q is maximized with respect to the Laplacian matrix as follows where the first three constraints ensure a valid Laplacian matrix and the last one avoids trivial solution for c > 0 . The second term of the objective function promotes the signal smoothness over the estimated topology [as explained in (4)]. The optimization problem in (23) is not easy to implement due to the first and third terms. To solve the issue with the pseudo determinant term, we propose the following minimization problem where det (⋅) stands for the determinant and = 1 N T . The equivalence of logdet ( + ) and log| | has been justified in [22]. The matrix inverse term which needs high computational power can also be replaced by a less computationally expensive term. Theorem 1 of [56] proposed upper and lower bounds for the trace of symmetric positive definite matrices inverse as and thus for a fixed , minimizing ‖ ‖ 2 F results in the minimization of Tr ( + ) −1 . To sum up, we propose to solve the following minimization problem for the Laplacian matrix estimation argmin − logdet ( + ) + Tr( ) + Tr ( + ) −1 where ‖ ‖ 2 F can also be considered as the control term for the distribution of off-diagonal elements, i.e. the edge weights of the estimated graph. Since (27) is a convex optimization problem, any off the shelf convex solver may be used, e.g. YALMIP [57]. However, in the next section, a proximal point algorithm is proposed to solve (27) efficiently. To conclude this section, the three steps of the proposed method are presented in Algorithm 1.

The efficient implementation
In this section, we discuss milestones to implement the proposed method with a fast and efficient algorithm. The proposed algorithm to solve (27) follows these steps: first by Proposition 2, the last three constraints in (27) are rewritten in the form of inner product of two matrices, helping to form a simpler Lagrangian function. Then, the Lagrangian is iteratively maximized with respect to the dual variable and then minimized with respect to the primal variable, i.e. the graph Laplacian matrix. To do the minimization step, we apply a proximal point algorithm. The maximization can be done by the Newton or quasi-Newton methods where the derivative and Hessian with respect to the dual variable are obtained by Lemma 2 and 1. Theorem 1 finds the stopping criterion guaranteeing the convergence of iteration between the maximization and minimization problem.

Proposition 1 Since is a symmetric matrix, the summation over all entries of the i'th row (or i'th column) is rewritten as where i is an N × N matrix in which the i'th column is the all one vector and the rest is zero.
Proposition 2 Having, ⋅ = , the constraint ij ≤ 0 for i ≠ j in the problem (27) can be replaced by Proof The weight matrix is a matrix with zero diagonal entries as long as we have assumed that there is no self loop. If all edges are positive (and so ij ≤ 0 for i ≠ j ), we have ‖ ‖ 1 = ‖ ‖ 1 and If there is only one edge with negative weight (or correspondingly ∃i ≠ j , ij > 0 ), ‖ ‖ 1 ≠ ‖ ‖ 1 . Note that if all edges are negative and thus ij ≥ 0 and ii ≤ 0 , this property is proved in the same way. However, this is not the case here as long as it is assumed that ∈ S N + , forcing the diagonal entries of to be non-negative. ◻ By applying the Propositions 1 and 2, the minimization problem of (27) is rewritten in the following form (27) argmin − logdet ( + ) + Tr( ) Hereafter, we interchangeably use the inner product of two matrices instead of the trace operator, i.e. Tr( 1 ⋅ 2 ) = ⟨ 1 , 2 ⟩ . It is assumed that the feasible set F = { ∈ S N + , B( ) = } is not empty and then due to the convexity of (30), any local optimal solution is also global optimal. The Lagrangian function over the primal variable and the dual variable = [ 1 , … , N , N+1 , N+2 ] is given as below Since the objective function of (30) is convex and there exists a strictly feasible point (in F ), Slater's condition holds and there is a strong duality or saddle point property [58]. Strong duality means that it is possible to find the optimal of the dual problem instead of the primal one. Thus, can be estimated via the following optimization problem where However, it is difficult to implement (33) directly. Therefore, we propose to use the proximal algorithm, which is an standard tool for solving constrained and nonsmooth minimization problems [59]. The proximal operator of a closed proper convex function g ∶ R n → R ∪ {+∞} is represented by g ∶ R n → R n and its scaled version with parameter is given as below One of the tools for solving the general optimization problem ̂ = argmin g( ) is the proximal point algorithm, also called proximal iteration, given by where t is the iteration index. In this algorithm, k and g k converge to the set of minimizers of g and its optimal value, respectively [60]. The proximal operator provides a smoothed version of the function by adding the regularized term. The proximal operator can also be computed by [59] where ∇ is the gradient operator and M g ( ) is the Moreau-Yosida regularization. The Moreau-Yosida regularization of g( ) in (33) with parameter > 0 is defined as [59,61,62] where (a) = follows from Von Neumann-Fan minimax theorem [63,64]. Like the proximal operator, the Moreau-Yosida regularization provides a smooth version of the function. Moreover, M g ( ) and g( ) have the same set of minimizers (30) argmin and thus we solve (38) instead of (34), equivalently. In particular, we propose to find the minimum of M g ( ) by applying the proximal point algorithm in (36).

Finding ;
The adjoint operator of B , called hereafter B a , is obtained as follows where denotes the partial derivative with respect to and (a) follows the inner product property. By using (31) and some simple manipulations, we have To simplify (38), we introduce the change of variable as and then ( ; ) can be rewritten as follows where To find the minimizer of J ′ , ; and simplify (42), the following Lemma from [65] will be used. 1. = 1 − 2 and 1 2 = 2. + is continuously differentiable and its derivative at for every ∈ S N is given as where • denotes the Hadamard product and ∈ S N is defined as follows

3.
Proof See [65]. ◻ To find the minimizer of J ′ , ; , the derivative with respect to ′ is set to zero as follows and by some simple manipulations, the solution is � = + � ( ; ) where � = 1+2 . Then (42) is simplified as follows (see Appendix 2)

Find M g ( ) via (38)
Lemma 2 The derivative of ( ; ) with respect to is By taking the derivative of ∇ ( ; ) with respect to and applying (45), the following corollary follows.

Corollary 1 The Hessian of ( ; ) with respect to is
Using the first and second order derivatives of ( ; ) with respect to , the unconstrained maximization problem of (38) can be solved by Newton or quasi-Newton methods, like L-BFGS. Let opt denote argmax ( ; ) and by using (38),

Lemma 3 The derivative of M g ( ) with respect to the graph Laplacian matrix is
Proof The result follows by taking the derivative of (49) with respect to and applying Lemma 1, part (3), similar to the proof of Lemma 2. ◻ Considering (36), (37), and (52), the graph Laplacian matrix is estimated via the following iteration rule where (t) stands for the iteration index. All steps to update the Laplacian matrix are summarized in Algorithm 2. It is also possible to use t rather than to show the possibility of adjustment in each iteration for a faster convergence of the objective variable . To find t in each iteration, the exact or backtracking line search can be applied [58,59].

Convergence analysis
Assume that the eigendecomposition of ( ; ) is represented as = T . Hence, we have where the scalar function + (⋅) is a convex function and also bounded with bounded eigenvalues. From theorems B and C in [66], it follows that + ( ) is Lipschitz continuous. Thus, where l c is the Lipschitz constant and (t) = (t) ; (t) opt . 3 2 .

Lemma 4 The scalar function + (x) is Lipschitz continuous with the Lipschitz constant
Proof See Appendix 4. ◻

Lemma 5 If is a real symmetric matrix with the eigenvalue matrix
Proof See Appendix 5. Proof The proof follows readily by applying the eigenvalue decomposition (44) and using Lemma 4 and 5. ◻

Theorem 1 The proposed BTL-PPA converges when the following stopping criterion is set for dual variables update
where the optimum primal and dual variables are (t * ) and (t * ) , respectively.

Numerical results
The proposed algorithm is tested for simulated and real data for different scenarios. For topology learning performance, the results of our algorithm are compared to those of three existing algorithms: GL-SigRep [51], CGL [22], and the learning sparse graph algorithm in [20], called LSG here. In synthetic data simulation part, although the performance of the LSG algorithm is similar to the BTL-PPA, it has a higher computational complexity compared to the other three algorithms.
The results for signal recovery performance are only compared to those of GL-SigRep, because the other two methods have no policy for signal representation and recovery. To implement GL-SigRep, an optimal selection of the regularization parameters via exhaustive search is applied and thus the appropriate ratio is found in order to maximize the performance for the algorithm in [51].

Synthetic data
The synthetic data is drawn from a Gaussian Markov Random Field process. • Normalized Mean Squared Deviation of graph topology estimation: , where ̂ denotes the estimated Laplacian matrix, • Normalized Mean Squared Error of signal reconstruction: , where TP, FP and FN are the numbers of true positives, false positives, and false negatives, respectively. Also, precision is the number of correctly recovered edges to the number of reconstructed edges in the estimated graph and recall is the number of correctly recovered edges to the number of edges in the ground-truth graph. This performance measure solely takes into account the support of the recovered graph while ignoring the weights, • Normalized Mutual Information: NMI( ,̂ ) = 2MI( ,̂ )

H( )+H(̂ )
, where MI is the mutual information between the set of edges in the estimated graph and the true graph and H( ) is the entropy based on the probability distribution of edges, i.e. the probability of zeros and ones in the matrix [70].  Figure 4 shows that BTL-PPA outperforms GL-SigRep for different noise variances and its NMSE are lower than those of GL-SigRep. Considering the NMSD comparison, the BTL-PPA topology learning algorithm works very well compared to other algorithms and its topology learning capability is robust in different noise powers. The experiments are run for difference values of noise variances, for a fixed signal power. Hence, the horizental axes can also be linked to different signal to noise (SNR) ratio. As shown in Fig. 4, when the noise increases, the NMSD for CGL algorithm decreases which seems to be counter-intuitive at first glance. However, it is necessary to note that for the current chosen threshold to remove weak edges, we get more edges leading to a lower NMSD (due to less variance in edge weights) and at the same time, as shown in Fig. 5, F-measure is decreasing due to an increasingly lower precision. Figure 5 also corroborates the strong performance of the proposed algorithm from the perspective that how much percentage of edges are learned correctly.

Signal inpainting in IoT applications
In the IoT framework, a challenging area of research is preprocessing of data after sensing and before providing it to further processing (for more detail, please see [1] Figure 3). One of these preprossesing steps is to estimate the missing values, i.e. signal inpainting. These missing values may be due to faulty recording, signal transmission to the central unit or storage. Assume the sensor readings are recorded in the vector = M ; U , where M is the portion of the signal that is correctly recorded and known while U includes missing values. Following the procedure explained in [71] but with our designed filter in (18), the signal is estimated as follows where M is the cardinality of M .
To apply the proposed methods on a relevant IIoT data, the "Household Power Consumption" data set has been downloaded from [72]. This data set is an individual household electric consumption data set collected via submeters in 5 years (2006-2010), while we only use a subset of it here. First K = 100, 000 time samples have been saved from the first block of data and then it is downsampled to 1000 samples. The data set contains 7 measurements, including the active power, reactive power, voltage, intensity and three energy submeter recordings from different places in a house. The intensity sensor readings are chosen to be tested here and then all multipliers of 50 in its time instances were deleted to model missing values, i.e. y k fork ∈ {50, 150, … , 950, 1000} are the missing values. Now, this dataset is given to the Algorithm 2 to find the underlying graph Laplacian matrix along with the estimated noise variance. The missing values are estimated by applying (57). Figure 6 shows the result of signal inpainting via our proposed filter. For a better illustration, only the missing values and the estimated ones have been shown. The main reason for the good performance is that the designed filter considers each dimension of the multi-variate signal as an entity of the entire network which is related to others based on a structure and hence it tries to use this knowledge in the filtering process. In other words, based on the measured signals, the algorithm learns the structure of the sensor readings first and then tries to use this data structure to interpolate the missing values.

Temperature data
In this experiment, the daily temperatures of N = 48 states of the USA mainland are stored for the years 2011 to 2014, i.e. K = 1461 [73]. Here, the graph signals are average daily temperatures and the underlying topology describes the temperature relation between states. We do not have access to the ground-truth topology, but a geographical based graph may be proposed to compare the results. A graph is considered where the nodes are the states and an edge weight between two states is computed by the Gaussian RBF of the physical distance.
First, the temperature data of 2011 is used to learn the underlying topology. Figure 7 compares different graph learning algorithms with respect to their capability to capture the underlying connections. It is shown that the proposed BTL-PPA can detect most of the edges. Figure 8 shows the topology learned by our proposed method over the real map of the United States. It can be inferred that a state weather influences that of the neighboring states, which is also corroborated by physics of temperature propagation. There are more edges in the far right side of the map due to the higher density of nodes, representing east coast regions. Then, from right to left of the map, there is a bit discontinuity among edges which can be representing the blockage of weather propagation by Appalachian mountains. The smaller number of connections in the middle of the figure can also be due to the Rocky mountain chain, also mentioned in [6]. At the same way, a few edges in the right side of the California corroborates the effect of Sierra Nevada mountain range. The numerical results to compare different algorithms' performances are given in Table 2.
In the second experiment with real temperature data, the data set is divided into two groups of data, the training set and the testing set. The first half of the data, i.e. the averaged daily temperature of 2011 and 2012 k = 1, … , 731, 4 is utilized as the training data to learn the underlying topology. Then we consider this topology as the ground-truth graph and the remaining data, i.e. the data from the year 2013 and 2014, is used to estimate the topology and compare to the ground-truth to check the consistency of the learning algorithms. In other words, a cross validation procedure is done to verify how the learned topology from given training data is different from the one learned by the test data. Table 3 shows the results for the cross validation scenario. The proposed BTL-PPA algorithm has the best consistency results among all algorithms with respect to all performance measure.

Conclusion
Many information networks involve multiple interacting entities for which finding the topology connecting these entities connections is somehow important for real-world applications. In the graph topology inference framework, we tried to estimate the structure of the underlying data from multi-variate measurements. In other words, we investigated an algorithm to explore the link between the signal model and the graph topology. In this paper, a factor analysis model was used for signal representation in the graph domain and a Bayesian inference method was applied to learn the Laplacian matrix (which can Fig. 8 The learned graph topology via BTL-PPA algorithm from real temperature data in 2011: the background map is the USA mainland  uniquely represent the graph topology) and to estimate the noise variance at the same time. To formulate the problem, we used a Bayesian framework and proposed a minimum mean square estimation approach to denoise the measurements. Finally, a convex optimization problem over the graph Laplacian matrix was proposed and solved via a proximal point method to estimate the topology from a denoised version of graph signals. The experimental results corroborate the performance of the proposed algorithm for a wide range of sensor networks and IoT applications.