Distributed Spatial Filtering by a Two-Hop Consensus-Type Algorithm

In this letter, we discuss distributed spatial filtering (DSF) on networked systems to obtain signal values with a desired spatial frequency characteristic from those assigned to nodes by a distributed algorithm. We present a two-hop consensus-type algorithm for DSF based on an existing one-hop algorithm. We prove that the range of the filter characteristics the presented algorithm can achieve is broader than that for the existing algorithm by deriving a necessary and sufficient condition for achieving DSF. Simulation results show that our filtering algorithm and a new filter characteristic it provides are effective in distributed anomaly detection by sensor networks.

Distributed Spatial Filtering by a Two-Hop Consensus-Type Algorithm Shinsaku Izumi , Member, IEEE, Xin Xin , Senior Member, IEEE, and Taiga Yamasaki Abstract-In this letter, we discuss distributed spatial filtering (DSF) on networked systems to obtain signal values with a desired spatial frequency characteristic from those assigned to nodes by a distributed algorithm. We present a two-hop consensus-type algorithm for DSF based on an existing one-hop algorithm. We prove that the range of the filter characteristics the presented algorithm can achieve is broader than that for the existing algorithm by deriving a necessary and sufficient condition for achieving DSF. Simulation results show that our filtering algorithm and a new filter characteristic it provides are effective in distributed anomaly detection by sensor networks.
Index Terms-Control of networks, distributed control, filtering, sensor networks.

I. INTRODUCTION
I N the systems and control community, networked systems, constructed by connecting subsystems (or nodes), have been a major research topic. This is motivated by the fact that modern engineering applications, e.g., swarm robots and sensor networks, can be categorized into networked systems.
Herein, we focus on distributed spatial filtering (DSF) on networked systems. DSF is to transform signal values given for nodes into ones with a desired spatial frequency characteristic in a distributed way. Fig. 1 shows an example for high-pass filtering, where x i denotes the signal value given for node i and it is assumed that nodes are closer to each other as their indices are closer. The nodes transform the given signal values for reducing the amplitudes of the low-frequency components in a distributed way. The main advantage of DSF is that the spatial properties of signals to be processed can be used. For instance, the spatial frequency of temperature is generally low because it takes similar values at close locations. With this property, we can reduce the noise in temperature measurements from a sensor network through spatial low-pass filtering. Further, some anomalies (e.g., fires) increase or decrease temperature at specific locations and make the spatial frequency of the measurements high, and thus spatial high-pass filtering allows to detect the anomalies. Izumi et al. [1] proposed a DSF method by focusing on the relation between a consensus algorithm and signal processing on graphs [2]. Then, in [3], a more sophisticated DSF method was proposed and verified through experiments with a real sensor network. However, the DSF method in [3] has the drawback that the achievable filter characteristics are limited to ones described as polynomials with non-zero real roots. This limitation is severe because the degree of the describing polynomial corresponds to the degree of freedom in choosing the filter characteristic, while high-degree polynomials generally have complex roots.
We thus aim to develop a DSF method to overcome this limitation. Our contributions are as follows. First, we present an extended version of the DSF method in [3]. By focusing on the fact that the filtering algorithm in [3] is a one-hop consensus-type algorithm, we extend the existing algorithm to the two-hop version where each node uses information on its two-hop neighbors, i.e., the neighbors of its neighbors. We then prove that DSF is achieved using our algorithm even for the filter functions described as polynomials with complex roots. This overcomes the above limitation and makes the range of the achievable filter characteristics wider. Second, we demonstrate the effectiveness of our DSF method through its application to anomaly detection by sensor networks. We choose a filter function that cannot be handled by the existing algorithm [3]. We then show that our DSF method with this filter function achieves higher detection performance than that when using the existing algorithm.
Finally, we discuss related works. This letter is related to several topics: spatial frequencies for distributed systems [4], signal processing over graphs [5] and its applications [6], and distributed control with multi-hop communication [7], [8].
In [4], Bode integrals based on spatial frequencies were introduced for distributed systems. Mollaebrahim and Beferull-Lozano [5] investigated the linear transformations of signals over graphs, and Yi et al. [6] considered average consensus with the filtering of signals over graphs. In [7], [8], coverage and consensus with multi-hop communication were discussed, respectively. However, the purposes of these works are different from ours.
Notation: We describe the real number field and the set of positive real numbers using R and R + , respectively. Let I be the identity matrix. For the numbers x 1 , x 2 , . . . , x n ∈ R, we use diag(x 1 , x 2 , . . . , x n ) to represent the diagonal matrix with x 1 , x 2 , . . . , x n on its diagonal. We define [ The cardinality of the set S is denoted by |S|.
For an undirected graph, we represent its graph Laplacian by L; i.e., L is defined as the difference between the degree matrix and adjacency matrix of the graph.

II. EXISTING RESULTS ON DSF
We first review the existing results [3] on DSF.

A. Preliminary: Spatial Filtering of Graph Signals
Consider an undirected graph with n vertices, which is described by G = (V, E) for the vertex set V := {1, 2, . . . , n} and the edge set E. We suppose that a signal value is assigned to each vertex of G, and describe the signal values for n vertices using s ∈ R n . Then, (G, s) is referred to as a graph signal, where the edges of G specify the connections between the signal values.
For the graph signal (G, s), its Fourier transform (i.e., the graph Fourier transform) is given as follows. For the graph Laplacian L of G, we define its i-th smallest eigenvalue in the sense of the modulus by λ i (i ∈ V). Then, the graph Fourier transform f ∈ R n is written as for the orthogonal matrix V ∈ R n×n such that V LV = with := diag(λ 1 , λ 2 , . . . , λ n ). Because L is a symmetric matrix by its definition, we can always find a V that satisfies V LV = for a given L. The graph Fourier transform f provides an expression for (G, s) in spatial frequency domain. Specifically, f expresses the magnitude of the differences between signal values connected by the edges of G. In (1), λ i (i ∈ V) corresponds to the spatial frequency, and the component of λ i is denoted by the i-th element of f , where it is noteworthy that λ 1 , λ 2 , . . . , λ n are nonnegative real numbers because L is the graph Laplacian of the undirected graph G.
The graph Fourier transform (1) allows the spatial filtering of graph signals described bỹ wheres ∈ R n represents the signal values of a filtered graph signal and h : R + ∪ {0} → R is a function to characterize the filter. From (1) and the orthogonality of V, i.e., V = V −1 , (2) implies that the spatial filtering is to filter a signal transformed into one in the spatial frequency domain by (1) using h and to perform the inverse transform of (1). In this way, we can obtain ans such that the graph signal (G,s) has a spatial frequency characteristic specified by h.

B. Problem Formulation
Consider the networked system having n nodes. Node i (i ∈ V) is given as the discrete-time model where x i (t) ∈ R and [x j (t)] j∈N i ∈ R |N i | are the state and the input, respectively, and N i ⊆ V is the index set of nodes of which node i has access to the information. The function g : . .} → R specifies the behavior of the node. Using g rather than g i means that the behavior of all the nodes is determined by a common function, which allows to make scalable. We represent the network topology of the system using the undirected graph G defined above, where the vertex set V and the edge set E express the indices of n nodes and the connections between them, respectively. In this case, the set N i is given by can be regarded as one of graph signals described in Section II-A by considering G and x(t) as the graph and the signal values, respectively.
Based on this observation, our interest here is to construct a system that works as a desired spatial filter for graph signals. To this end, we regard the initial state x(0) and the final state x(∞) as the signal values of an original graph signal and those of the filtered one, respectively, and address the following problem.
Problem 1: For the networked system , suppose that a filter function h and an initial state x(0) are given. Find a function g (i.e., a distributed algorithm executed by n nodes) to produce an x(∞) such that the graph signal (G, x(∞)) has a spatial frequency characteristic specified by h.

C. Filtering Algorithm
The idea of [3] to solve Problem 1 is to design the function g so that the relation between x(0) and x(∞) is equivalent to the spatial filtering (2) of graph signals when regarding x(0) and x(∞) as s ands, respectively. Based on this, was provided, where 0 (t) ∈ R \ {0} and 1 (t) ∈ R are timevarying gains that satisfy 0 (t) = 1 and 1 (t) = 0 when t ≥ m for a positive integer m. The provided g and (3) yield a consensus-type distributed algorithm. It should be noted that this algorithm converges in m timesteps due to 0 (t) = 1 and of the graph signal (G, x(0)) for h. Then, the following result was obtained in [3].
Lemma 1: For the networked system , suppose that a filter function h is given and let g be given by (4). Then, achieves DSF of the graph signal (G, x(0)) for h if and only if h is a polynomial of the form with a frequency variable λ, where a 0 , a 1 , . . . , a m ∈ R and the roots of the polynomial must be non-zero real numbers. Lemma 1 indicates that DSF is achieved using (3) and (4) if and only if the filter function h is a real polynomial having non-zero real roots.

D. Problem to Be Considered
Although Lemma 1 gives a necessary and sufficient condition on h for achieving DSF, it is difficult to achieve some filter characteristics under this condition. An example is shown in Fig. 2. The thick blue line and the thin red line indicate h(λ) := 1/(1 + e −20(λ−0.2) ) and its polynomial approximation h 10 (λ) ≈ 0.000149λ 10 − 0.00453λ 9 + 0.0583λ 8 − 0.413λ 7 + 1.74λ 6 − 4.41λ 5 + 6.15λ 4 − 3.21λ 3 − 2.18λ 2 + 3.27λ + 0.0100 1 via a least square fitting, respectively (and the green circles are detailed in Section III-C). As demonstrated later, h 10 describes a filter characteristic that is useful in an application to sensor networks. However, we can numerically confirm that h 10 has complex roots and does not satisfy the condition in Lemma 1. Therefore, we seek the polynomial approximation with non-zero real roots by solving an optimization problem for fitting using the function "fminsearch" in MATLAB, and obtain h r 10 (λ) ≈ −(4.37 × 10 −13 )λ 10 + (2.60 × 10 −11 )λ 9 + (1.
However, the plot of h r 10 indicated by the dashed black line in Fig. 2 is quite different from that of the original h, compared with the case of h 10 .

III. DSF BY TWO-HOP CONSENSUS-TYPE ALGORITHM
In this section, we overcome the aforementioned difficulty by extending the results in [3].

A. Proposed Algorithm
Consider the networked system again. We suppose here that the model of node i (i ∈ V) is given by 1 The coefficients of the polynomials are rounded off for brevity.
where [y j (t)] j∈N i ∈ R |N i | is the input, y i (t) ∈ R is the output, and g 1 : Our approach to the modified problem is to extend the existing filtering algorithm given by (3) and (4) to the twohop version, i.e., a consensus-type algorithm that uses the information on two-hop neighbors on the graph G in addition to that on the (one-hop) neighbors. Based on this, we propose where 2 (t) ∈ R is a time-varying gain satisfying 2 (t) = 0 for t ≥ m. In the proposed algorithm given by (6)-(8), the output y i (t), transmitted to the neighbors of node i, contains the state x j (t) (j ∈ N i ), and thus the nodes can obtain the information on their two-hop neighbors.

B. Main Result
For the proposed algorithm given by (6)-(8), the following main result is obtained.
Theorem 1: For the networked system where the node model (3) is replaced by (6), suppose that a filter function h is given and let g 1 and g 2 be given by (7) and (8), respectively. Then, achieves DSF of the graph signal (G, x(0)) for h if and only if h is a real polynomial of the frequency variable λ with at most degree 2m and non-zero roots.
Proof: See the Appendix. Theorem 1 shows that DSF is achieved using the proposed algorithm given by (6)-(8) even if the filter function h is a polynomial with complex roots. This makes the range of the achievable filter characteristics wider. Here, the gains 0 (t), 1 (t), and 2 (t) are chosen as, for instance, (12)-(14) or (12)-(14) and (18) (see the Appendix). Further, Lemma 1 ensures that the existing one-hop algorithm given by (3) and (4) cannot handle h other than polynomials with non-zero real roots. In this sense, Theorem 1 implies that the extension to the two-hop algorithm is a key to handling the case of the complex roots.
From the above discussion, our method to solve Problem 1 with the design parameters g 1 and g 2 is as follows: (i) approximate the given filter function h by a polynomial with non-zero roots; (ii) give g 1 and g 2 using (7)  Final state x(5) given by the proposed algorithm and the collective signal values given by (2) for s := x(0). m := 5. Let us recall that h 10 is a polynomial with complex roots and thus cannot be handled by the existing algorithm given by (3) and (4). From h 10 , (6)- (8), and (12)-(14), the algorithm for the nodes is obtained, where we use the fact that the degree of h 10 is even. The resulting final state x (5), which is equal to x(∞) (see the Appendix), is depicted in Fig. 3, where the thick red lines represent the elements of x(5) and the blue circles and lines represent the vertices and edges of G, respectively. Moreover, the thin green lines indicate the collective signal values given by (2) for s := x(0), wheres i is the i-th element ofs. We see that x(5) ands are the same although the proposed algorithm is distributed but (2) is not. This demonstrates our results.
Remark 1: The degree of the polynomial approximation of the filter function h should be chosen according to the guideline given in [1,Remark 6]. Note here that we do not need to consider the constraint that the resulting polynomial has no complex roots, unlike the case of [1].

C. Application to Distributed Anomaly Detection by Sensor Networks
Next, we apply the proposed DSF method to distributed anomaly detection by sensor networks, i.e., finding sensor nodes such that the measurements are quite different from those from other nodes in a distributed way. As explained in Section I, an example of anomaly is a fire, which causes the differences between temperature measurements.
Our idea to achieve the anomaly detection is to remove only the component of the zero frequency at which the signal values for all the nodes are the same. By removing the zero frequency component, we can extract the differences of the signal values from a typical one in the group of the nodes, which allows the anomaly detection. In the following, we show that the proposed algorithm can handle even a specialized filter function to remove only the specific frequency component, and demonstrate that the anomaly detection can be achieved by the proposed algorithm and the filter function.
Consider the example in Section III-B again. We represent the eigenvalues λ 1 , λ 2 , . . . , λ 8 of the graph Laplacian L of the graph G depicted in Fig. 3 by the green circles in Fig. 2. These eigenvalues correspond to the spatial frequencies of which the associated graph signal has the components, and   Fig. 2 indicates that the filter function h 10 has the property of removing only the zero frequency component of graph signals with G. As seen above, the proposed algorithm given by (6) Fig. 4, where node 4 has to be detected. It turns out that the elements other than x 4 (5) are close to zero. Thus, by introducing a threshold, node 4 can be detected. The thin green lines in Fig. 4 represent the corresponding result when using the existing algorithm given by (3) and (4) for the filter function h r 10 in Section II-D. We see that the elements other than x 4 (5) are not sufficiently reduced compared with the result for h 10 . This implies that the anomaly detection with the existing algorithm for h r 10 is difficult because the choice of the threshold heavily depends on the final states of nodes other than those at which anomaly occurs.
We evaluate the detection performance of the proposed algorithm for h 10 in the following way. We suppose that anomaly occurs at a node j and j is randomly chosen from the set V with equal probability. The initial state x j (0) of node j is set as x j (0) := b+ +w, where b, , w ∈ R are a base value, measurement noise, and disturbances due to the anomaly, respectively. The noise is Gaussian noise with mean 0 and variance 0.5, and w takes ±10 randomly with equal probability. For the other nodes, let x i (0) := b + (i ∈ V \ {j}). Then, we say that the anomaly detection is successful if only |x j (m)| is larger than the threshold d ∈ R + . Under this setting, the success rates for 1000 trials are summarized in Table I, where by considering the differences in the algorithms and the filter functions, we choose the different thresholds d := 5 and d := 9 for the proposed algorithm and the existing algorithm, respectively. Table I indicates that the detection performance of the proposed algorithm is higher than that of the existing algorithm for all of b := 10, b := 20, and b := 30. The reason for this is explained as follows. As seen in Fig. 4, the existing algorithm for h r 10 cannot sufficiently reduce the states of the nodes other than that to be detected. As a result, the existing algorithm cannot handle the case where |x j (0)| is smaller than the other states. Therefore, the anomaly detection for w = −10 fails, demonstrated as the success rate 51.2 % for b := 10. In addition, the existing algorithm cannot handle the different values of b by the common threshold d because the remaining states depend on b. This leads to the differences in the success rates for b := 20 and b := 30.
Remark 2: From (6)-(8), (12)-(14), and (18), the proposed algorithm does not depend on the graph G. However, the resulting signal values depend on G because the eigenvalues of the graph Laplacian L, i.e., the distribution of the frequency components of the associated graph signal, change depending on G.

IV. CONCLUSION
In this letter, we improved an existing DSF method for networked systems to overcome the difficulty of the limited achievable filter characteristics. In particular, we extended the existing one-hop filtering algorithm to the two-hop version. We then derived a necessary and sufficient condition on the filter function for achieving DSF. The advantage of the proposed algorithm is that the filter functions described as polynomials with complex roots can be handled, unlike the existing algorithm.
Next, we give the proof of the "if" part. It is assumed that a real polynomial of λ with at most degree 2m and non-zero roots is given as the filter function h. Because the complex roots of real polynomials come in complex conjugate pairs [9], we denote the complex roots of h by α j ±β j i (j = 1, 2, . . . , m c ), where i := √ −1 and m c ∈ R + ∪{0} is the number of the pairs. Further, we let α j ∈ R \ {0} (j = m c + 1, m c + 2, . . . , μ) be the real roots of h, where μ := m c + m r for the number m r of the real roots. In the following, we assume that the degree of h is even. Then, m r is even from the above property of the complex roots of real polynomials.
With the above notations andμ := m c + m r /2, we set 0 (t) := a 0 if t = 0, 1 otherwise, where a 0 is the zero degree term of h, as shown in (5). From m ≥μ, these gains satisfy 0 (t) = 1, 1 (t) = 0, and 2 (t) = 0 for t ≥ m, where the inequality on m is derived from the facts that the degree of h is at most 2m and the number of the roots