AU2020103329A4 - A based on recursive least squares online distributed multitask graph filter algorithm. - Google Patents

Abstract

: An online distributed multitask graph filter algorithm based on recursive least squares (RLS) is used to design node-variant graph filter and promote the cooperation between filter coefficients. By introducing auxiliary variables, the algorithm is optimized by alternating direction method of multipliers (ADMM) method. The algorithm mainly includes the following steps: generating the filtered graph signal of a single node; generation of global objective function; the solution of distributed online ADMM. The invention is applicable to the case where the filter coefficients are different but the similarity is high in the graph filter with variable nodes, which can promote the cooperation of filter coefficients, and use RLS algorithm to improve the speed of filter estimation. The invention can be applied to the scene where the observation data exhibit non-Euclidean structures, such as sensor, transportation, communication, social or biological networks. In the simulation part of this patent, we test both single task and multi task situations, and introduce a practical example, i.e. the reconstruction of temperature in the United States. The simulation results show that the algorithm has a certain competitiveness compared with the existing distributed graph signal methods in the estimation speed and accuracy.

Description

Editorial Note 2020103329 There is 8 pages of Description only.

1. Background and Purpose In many current interesting applications, such as sensor, transportation, communication, social or biological networks, observation data often exhibit non-Euclidean structures, which makes the traditional signal processing technology unable to analyze them. As a consequence, in recent years, people are committed to developing new graphic signal analysis methods, thus forming the research field of graph signal processing (GSP). GSP extends the classical discrete-time signal processing tools to the signals in the discrete domain, where the vertices are connected by graphs. A characteristic of graph signal is that the signal domain is not a metric space. For example, in biological networks, the vertices may be proteins, genes, enzymes, etc., and the existence of edges means that these molecules undergo chemical reactions. In order to cope with graph signals, GSP relies on two ingredients: the graph shift operator (GSO) on one hand, which accounts for the topology of the graph, and the graph Fourier transform (GFT) on the other hand, which allows to represent graph signals in the graph frequency domain. Built upon the definition of the GFT, graph filters play a central role in processing graph signal spectra by selectively amplifying or attenuating frequency components. Various architectures of graph filters have been proposed, including finite impulse response (FIR) and infinite impulse response (IIR) filters. From a perspective of scalability, and considering energy constraints and band-limited communications that may be encountered in large networks of distributed nodes such as sensor networks, significant efforts have been made recently to derive distributed graph filters. These filtering procedures allow each node to exchange only local information with its neighboring nodes. Much of the GSP research has focused on static graph signals, that is, signals that need not evolve with time. However, a wide spectrum of network-structured problems requires adaptation to time-varying dynamics. Sensor networks, social networks, vehicular networks, communication networks, and power grids are some typical examples. Prior to the more recent GSP research, many earlier works on distributed networks have addressed problems dealing with this challenge by developing processing strategies that are well-suited to data streaming into graphs. Most of the existing distributed graph signal processing methods are based on LMS algorithm, which has the advantage of simple calculation and stable numerical value. However, the convergence speed of graph signal processing based on LMS algorithm may be slow due to the non-energy preserving property of shift matrix. In addition, for the node-variant Graph Filter, the filter coefficients are similar when the graph signal is smooth. We can consider multitasking in graph signal processing to promote the cooperation among filter coefficients of nodes. The purpose of this patent is to solve the problem of slow convergence rate caused by traditional methods in the design of distributed graph filter and how to cooperate to promote estimation when the filter coefficients are similar. A more flexible graph filtering model is considered, where each node in the graph seeks to estimate a local node-varying graph filter. This allows us to exploit more degrees of freedom in the filter coefficients to better model graph signals. In this patent, we are interested in online learning of linear graph models for representing streaming graph signals in a distributed manner. We focus on RLS algorithm, since RLS strategy can accelerate the convergence speed of graph filter design, and solve the disadvantage of slow convergence caused by LMS strategy. In this work, we are interested in online distributed learning of linear graph models without assumption of band-limited processes. We use graph filter models in the time-vertex domain where there is no need to decompose the graph shift operator. The formulated optimization problem relies on minimizing a global cost consisting of the aggregate sum of individual costs at each vertex. Our present invention proposes a based on RLS online distributed multitask graph filter algorithm. In addition, we also consider the application of multitasking to the graph filter design, which can promote the cooperation of filter coefficients when the filter coefficients are different but similar. Simulation results show the efficiency of the proposed algorithms and validate the theoretical models.

2. Graph Filter 2.1 Graph Signal Processing Tools

Let us consider a graph g= (N,, ),which is connected, undirected and weighted.

={1,2,..., N} and E denotes the set of nodes and the set of edges, respectively. If node n

and i are connected, then (n,i)e .The (n,i)-th entry wn of adjacency matrix

W={1wjN NxN represents the weight of the relationship between the vertices n and i.

The collection of neighbor nodes of any particular node n isrepresentedas AN,={i:(n,i)eJ.

Likewise, matrix S denotes the graph shift operator, whose entry s,1 i is non-zero only if

n = iorie n,. Generally, the shift matrix is chosen as Laplace matrix or adjacency matrix. A

signal over a graph G is defined as x=[x,...,xN E RNwhose n-th element x,

represents the signal sample of node n .

2.2 Graph Filter Step 1: Generationfilteredgraph signal. Classical graph filters are called (node-invariant) GFs, which is defined as: M-1 HA h°S'" (1) m=O

where h° ={h"I} represents the vector of coefficients, which is the same for all nodes, and S is the shift matrix. It's obvious that GFs are polynomials of the graph-shift operator. In this patent, we focus on the node-variant Graph Filter, which is defined as:

H m diag(h('))Sm (2)

where h('") representing the vector of coefficients, which is different for each node. Filtered graph signal y(i) is generated from the input graph signal x(i)asfollows:

y(i)=Jdiag(h('))S'"x(i-m)+v(i), i >M-1 (3) m=O

where v(i) represents the i.i.d. zero-mean noise at time i. Noting that here we consider a more general model that embeds the time dimension, that is, the m-hop spatial shift S'" is now performed in the m time slots. Step 2: Preservingshift signal matrix. If we keep the following matrix composed of the shift signals:

X(i) A [x(i -1), Sx(i - 2),..., S M -'x(i - M)] (4)

Only one shift at time i is needed to generate the filtered graph signal y(i). Let P(i)

denote the Nx M matrix given by: P(i) A [x(i), Sx(i -1),..., S M -1 x(i - M +1)] (5) then the model (3) can be replaced by: y(i)= P(i)diag(h ")+ v(i) (6) Step 3: Generating thefiltered graph signal of a single node. Model (6) can be split into the form of each sample y, (i) at node n :

yn(i)= p,(i)h"+ vn(i) (7)

where h" denotes the filter coefficient vector at node n collected into h(m),i.e., [h,], =[h'")]n,

and pnT(i) isthe n-throwof P(i) givenby:

p,(i) col{[x(i)],[Sx(i -1)],...,[S M 'X(-M +1'An (8)

3. Generation of Global Objective Function Step 1: multitask and RLS strategies in the graphfilter We propose an RLS based estimator to track the unknown filter coefficient vector while enforcing similarity between neighbors' weight vectors. Let us consider the data model (6), and without losing generality, suppose {y(i),x(i),v(i)} are zero-mean jointly wide-sense stationary random processes. The estimator at time t is the optimal solution of the following optimization problem: N T N Minimizeh,,hN ZZ '(yn(t)-p n(t)ho°) 2 +p ||h-h| (9) n=1 t=1 n=1 iE

where h,,0<A<1,§>0 represents the true filter coefficients, the forgetting factor of RLS

algorithm, and the regularization coefficient of similarity between nodes and neighbors, respectively. Step 2: usingADMMmethod to solve global costfunction. In order to solve problem (7) in the form ofADMM, we can transform it into the form of ADMM by introducing auxiliary variables and some linear constraints. Therefore, problem (7) can be rewritten as follows:

Minimizeh,. hN 1 '.y t ...p t n=1 t=1 n=1 iE~l

s.t. hn= x,,n =1,., N, (10) h, =vn =1,., N, i=1,.., I

where x, and v.,,Vn 1,..., N,i E An are auxiliary variables.

3.The Distributed Online ADMM Step 1: form the augmented Lagrange. we introduce an extended ADMM algorithm called distributed online ADMM algorithm. In order to solve problem (10) in the form of ADMM, by introducing Lagrange multiplier 7, and u , as well as a positive constant p ,we can form augmented Lagrange of (7) as follows:

ZP (th,, , ., , yp.,)n=1,.,Nji=1_,.,|N |

T 2 II ,(y(t) p (t)T h,) + 1 -V12 (1) n=1 t=1 n=1 ie A.

2 +Zy (h h X))+ u((hi --x| N -v.1| n=1 n=1 iE , 2n=1 n is

. Then, the update steps of ADMM to solve problem (17) are as follows:

hkI = arg min2p(h,xk" 3V Yk " ) (12)

{xAVk+}= arg min2, (hk*kIxvk,,uk) (13)

7n =7 + p (h* - xk (14)

k+iI k ki- _ Vk- p =p,+ pI h -V (h , (15)

Step 2: calculate matrices and vectors about RLS. Since we use RLS algorithm in the cost

function, each node first calculates the correlation matrix Rn(T) andvector ,(T).

T t R,(T) =n-'p(t)p Q) (16) t=1

r(T)= AT,(tP.(t) (17)

Step 3: update hn(T). Combining (11) with (12), we obtain that the update of h in the

ADMM as follows: N T N N h* -ag m InI A-(y (t) p (t)Th") 2 ±+Z'T(h"-x)+± uh n=1 t=1 n=1 n=1 ic IV"

+ h -xk + P1h, - vi 1 2n=1 2n=1 i. ,,

which can be split into a single node form

T 2 h, = arg minZh -tn n h ) Iy(h -x) +)T h, t=1 (19)

Here, we use the Online Algorithm with Varying T, i.e., only one ADMM update iteration is performed in each slot. The classical ADMM algorithm needs to be iterated step by step to get the optimal solution. Therefore, we should run ADMM iteration several times for each time T, which could lead to a huge amount of computation. This processing is similar to some existing adaptive algorithms, in which gradient descent is performed only once per slot. Specifically, we replace k with T-1 in equation (17), and get updates that are suitable for varying time T:

1 h,(T)=I (2r(T) -,(T -1) nRj )±p ~lv±II)1 (20) - pi,,(T -1)+px,(T -1)+p , v, (T -1)) ieo0

Step 4: broadcasts h.(T) . This step is a broadcast step, that is, after calculating filter

coefficients of each node h, (T), the node transmits h, (T) to its neighbor nodes.

Step 5: update auxiliary variable. Similar to updating h, (T) , combining (9) with (11), we

obtain that the update of in the ADMM as follows:

{x±,v~}argmin,~ x,~-v 2 N n=1 En 2 =1 n=1 iE N I~Ihg

+ 22 - h n =1 2n=1 LAN

which can be split into a single node form:

{x*{k1 = arg min 1x 1 -7 Tx U(

2 p k I PAT'I(22) +L1h - xn + : 1hi - n 2 ~ 2 jK ,2 then use the online method to update, i.e., we replace k with T-1 in equation (22), and get updates that are suitable for varying time T:

1 x,(T)= 2, v (T -1)+ y,(T -1)+ph(T) (23)

1 v .(T) = (2px (T -1)+ p, (T -1)+ phk*' (T -1)) (24) 2,8 + p

Step 6: broadcasts v,,(T). Similar to step 4, the node transmits v ,(T) to its neighbor

nodes. Step 7: update Lagrange multiplier The update of 7, and u, are also decomposed across nodes:

7k* = h* k+ X (25) k+1 k i+ + p (h? (26)

+ pn~= pn -vn

and use the online method to update:

7n,(T) = yn (T - 1) + p ( hn (T) - xn(T)) (27)

p.,i(T)=,,i(T-1)+ p( hi (T)- v.,(T)) (28)

4. NUMERICAL RESULTS In this section, we present some numerical results to evaluate the performance of the proposed RLS based multi task graph filter design algorithm, i.e. algorithm 1. Firstly, we simulate the single task scenario and multi task scenario respectively when the input data are i.i.d., and then consider a practical application scenario to illustrate the effectiveness of the algorithm. Numerical results verify the correctness of the theoretical analysis, which shows that the algorithm has faster convergence speed and better convergence performance. 4.1 single task case

We consider the single task scenario, that is, the graph filter coefficients h' of all nodes are

the same .The processes (x(i),v(i)} are assumed to be zero-mean Gaussian with:1)

=E[x T (i)x(i)]=diag{o ,}; 2) R, = E[v (i)v(i)]= diag{o,}n ; and 3) x(i) and

2 2 v(i) are independent of each other. The variances o, and oa, are generated from the uniform

distributions U2(1,1.5) and U(0.1,0.15), respectively. The graph filter order was set to L=3

and the ideal coefficients h' are generated from the uniform distributions U(0,1).Thisdata

model is applied to a sensor network with N= 60 nodes. The shift matrix S is set as the normalized

W adjacency matrix, i.e., S = . At this point, all the eigenvalues of S are smaller than

1. Therefore, the energy of the shifted signal Smx decreases as m increases. In the simulation, we compare the diffusion LMS strategy, the preconditioned LMS (PLMS) strategy and our proposed strategy. Simulated results were averaged over 500 Monte-Carlo runs.

For the diffusion LMS algorithm, we set the step size parameter as p = 0.05 . For the

preconditioned LMS algorithm, we set step size parameter as p =0.05. And for our proposed

strategy, we set parameters as A= 0.98,$= 0.9,p = 0.1. The network MSD performance of the three algorithms is reported in Fig.1. It can be observed from the Fig.1 that the convergence speed of our algorithm is faster than the other two algorithms, and the performance is better.

4.2 multitask case We consider the multitasking scenario when the input data is i.i.d., that is, the graph filter coefficients of nodes are different but similar. This data model is also applied to a sensor network, both with N = 60 nodes. The shift matrix was the normalized adjacency matrix, that is,

S= W . Graph signal x(i) and noise v(i) settings are the same as in the single task 1.1A. (W)

scenario. The graph filter order was set to L=3 , and we set parameters as

0.98,p =0.9,p =1.3. The ideal coefficients h' are generated as follows:

cos0 +1 h*=h"+r sinO6.-0.5 (29)

,sin,, +0.7)

Ok=2rc(k-1)/N+rc/8 (30)

We compare DLMS algorithm, non-cooperative LMS algorithm and our proposed algorithm. For DLMS algorithm and non-cooperative LMS algorithm, we set the step size parameter as

p =0.08. For our proposed algorithm, we set parameters as A = 0.98,p = 0.9,p = 0.1. The

network MSD performance of the two algorithms is reported in Fig.2. It can be observed from the Fig.2 that the proposed strategy still has fast convergence speed and good performance in the multitask scenario.

4.3 reconstruction on U.S. Temperature Dataset In this example, we considered A data set, in which the hourly temperature measurements of N = 109 stations in the United States in 2010 with T = 8759 hours were collected. An undirected graph illustrated in Fig.3, is constructed by using a 7 nearest neighbors' method, which relies on geographical distances. The graph signal at each vertex n represents the temperature observed at the n-th station.

The data is divided into two parts, the first Ttain= 6570 hours is the training set, and the last

hours is the test set. The purpose of this experiment is to learn a graph filter that minimizes the reconstruction error on the training set, that is, TamN M 2

min Yk(i)-1h k[S'x(i- m+1)], (31) i=1 k=1 m=1

where y(i) is the real ground temperature at time i,and x(i) represents the observed value after

sampling,i.e., x(i)=diag(1s)y(i). 1s is the set indicator vector, whose n-th entryisequal

to one, if node n is sampled, otherwise it is zero. S represents a fixed sample set and is represented by the red nodes in Fig.3. The graph filter order was set to L = 4. We run different algorithms on the training set to learn the graph filter coefficients, such as the multitask diffusion LMS, the multitask diffusion PLMS and our proposed algorithm. In the last 120 hours of the test set, we provide the real and reconstructed temperatures at the unobserved nodes (the black circle in Fig.3). The experimental results are shown in Fig.4. It can be seen from the Fig.4 that our algorithm has good reconstruction performance.

Editorial Note 2020103329 There is 3 pages of Claims only.

Claims (4)

1. The claims defining the invention are as follows: A based on recursive least squares online distributed multitask graph filter algorithm.

   (1) Graph Filter Design Process Step 1: Generationfilteredgraph signal. We focus on the node-variant Graph Filter, which is defined as: M-1 H=Zdiag(h('))S' (1) M=O

   where h('") and S representing the vector of coefficients and shift matrix, respectively. Filtered graph signal y(i) is generated from the input graph signal x(i) as follows:

   y(i)=Jdiag(h'))S'"x(i-m)+v(i), i >M-1 (2) m=O where v(i) represents the i.i.d. zero-mean noise at time i. Step 2: Preserving shift signal matrix. We keep the following matrix composed of the shift signals: P(i) A [x(i), Sx(i -1),..., S M -1 x(i - M +1)] (3) then the model (2) can be replaced by: y(i)= P(i)diag(h(m))+ v(i) (4) Step 3: Generating thefiltered graph signal of a single node. Model (4) can be split into the form of each sample y,(i) at node n : yn(i)= p,'(i)h + v(i) (5)

   where h, denotes the filter coefficient vector at node n collected into h('"),i.e., [h"],m =[h')],

   and p7T(i) isthe n-throwof P(i) givenby:

   M -x(i-M+1)],} (6) pn(i) col{[x(i)],[Sx(i -1)],...,[S

   (2) Generation of Global Objective Function Step 1: multitask and RLS strategies in the graphfilter We propose an RLS based estimator to track the unknown filter coefficient vector while enforcing similarity between neighbors' weight vectors. The estimator at time t is the optimal solution of the following optimization problem: N T N

   Minimizeh,,...,hN I A -p n=1 t=1 n=1 ieA,

   where h,,0<A<1,,8>0 represents the true filter coefficients, the forgetting factor of RLS

   algorithm, and the regularization coefficient of similarity between nodes and neighbors, respectively. Step 2: usingADMMmethod to solve global costfunction. In order to solve problem (7) in the form ofADMM, we can transform it into the form of ADMM by introducing auxiliary variables and some linear constraints. Therefore, problem (7) can be rewritten as follows:

   N T N 2 Minimizeh1 ,...,hN I-'(y"(t) - p"(t) T h,)2 + I x -v,|12 n=1 t=1 n=1 ieA/N

   S. t. hk=x,,n==1,...,9 N, (8) h, = v, ,n = 1,..., N, i= 1,... 91,,

   where x, and v,,Vn=1,...,N,ic Efn are auxiliary variables.

   (3) The Distributed Online ADMM Step 1:form the augmented Lagrange. In order to solve the cost function (8), by introducing

   Lagrange multiplier 7, and u,, , as well as a positive constant p ,we can form augmented

   Lagrange of (7) as follows:

   ,({ht xk , .,, y,,p, }n=1,...,Nji=1_,,Nj| y-, 2+px x-v (y)(t)p(t)Th,) 2 (9) n=1 t=1 n=1 icAf, N N ~N 2 N

   +Zr (h -xn)±+ n1 n=1 iEIV, uZn, -vn,)+ n=1 ||h - x 2 =1 iENi h - | Step 2: calculate matrices and vectors about RLS. Since we use RLS algorithm in the cost

   function, each node first calculates the correlation matrix R,(T) and vector ri (T).

   T R(T) = A 2p (t)p,(t) (10)

   r (T)= AT',(tP.(t)(11)

   Step 3: update h,(T). We use the solution of classical ADMM method, that is, fix other

   variables and update the hn(T) variable. hn(T) is the filter coefficients of each node at time T,

   which is calculated as follows:

   1 h((T)±= pl1 2R,(T) + p (1 +|Q2 )(2r(T)- 7, (T -1) 1,|(1

   - pi,, (T -1)+px(T -1)+p 0 v (T -1))

   Step 4: broadcasts h,(T). This step is a broadcast step, that is, after calculating filter

   coefficients of each node h(T), the node transmits h,(T) to its neighbor nodes.

   Step 5: update auxiliary variable. Similar to updating hn(T), we fix other variables and

   update auxiliary variables x,(T) and v,,(T) respectively.

   x, (T) = - 1 2$81 v `(T - 1) + y/,(T - 1) +ph, (T) (13) 2p|MIVn 1 |+ P " n

   vn,(T) = 2 (2p8x (T -1)+,u,,(T -1)+pk+l(T- )+ (14) )6 +p

   Step 6: broadcasts v.,(T). Similar to step 4, the node transmits v (T) to its neighbor

   nodes. Step 7: update Lagrange multiplier The update of 7I and u are also decomposed across nodes:

   7,(T)= yn(T -1)+ p( hn (T))- x(T)) (15)

   p.,i (T)=pn (T -1) + p(h, (T) - vn (T)) (16)

   Fig1. Network MSD performance for single task case

   Fig.
2. 2 Network MSD performance for multitask case

   Fig.
3. 3 Graph topology for the U.S. temperatures dataset.

   Fig.
4. 4 True temperatures and reconstructed ones at an unobserved node.