Publication in the conference proceedings of EUSIPCO, Lisbon, Portugal, 201

Djuric, Petar M.

Wang, Yunlong

English

ZENODO

ON SEQUENTIAL ESTIMATION OF LINEAR MODELS FROM DATA WITH CORRELATEDNOISEYunlong Wang and Petar M. Djuric´Department of Electrical and Computer EngineeringStony Brook University, Stony Brook, NY 11794Email: {yunlong.wang, petar.djuric}@stonybrook.eduABSTRACTIn this paper, we consider the problem of Bayesian sequentialestimation on a set of time invariant parameters. At everytime instant, a new observation through a linear model isobtained where the observations are distorted by spatiallycorrelated noise with unknown covariance, whereas in time,the noise samples are independent and identically distributed.We derive the joint posterior of the parameters of interestand the covariance, and we propose several approximationsto make the Bayesian estimation tractable. Then we proposea method for forming a pseudo posterior, which is suitablefor settings where estimation over networks is applied. Bycomputer simulations, we demonstrate that the Kullback–Leibler divergence between the pseudo posterior and aposterior obtained from a known covariance decreases as theacquisition of new observations continues. We also providecomputer simulations that compare the proposed method withthe least squares method.Index Terms— Bayesian inference, distributed estimation,unknown covariance, pseudo posterior1. INTRODUCTIONEstimation over cooperative networks has been widelystudied in the literature (e.g., [1, 2]), where the agentsestimate the state of nature in a distributed manner byexchanging information with neighbors. In [3, 4, 5],we address consensus-based distributed estimation of linearmodels within the Bayesian framework.In addressing consensus-based estimation methods overnetworks, it is important to reformulate the formation ofthe optimal posterior to be a function of summation ofcertain statistics. In [4] and [6], it is shown that by usingthis strategy, the original problem can be converted to aproblem of distributed summation, which is suitable for theconsensus method [7]. In [8], the authors consider a problemThis work was supported by NSF under Awards CCF-1320626 andECCS-1346854.where the observation noises are spatially correlated withknown covariance matrix. In this paper, we are layingthe grounds for distributed estimation to scenarios wherethe observation noise is spatially correlated with unknowncovariance, whereas in time, the noise is independent andidentically distributed. We focus on the required processingby an agent at every time step and show why it cannot obtainthe exact posteriors of the unknowns. Then we resort toa suboptimal approach that can be used to approximate theindividual agent’s belief to a suitable form. This belief canthen readily be used in settings of distributed estimation.More speciﬁcally, we study sequential estimation where atevery time instant an agent gets observations from severalsensors about a vector of time invariant parameters througha linear model. In order to get the posterior of the parameterof interests, the agent needs to marginalize out the unknowncovariance. The marginalization, however, is computationallyintractable. Therefore, we deﬁne a pseudo posterior as anapproximation of the posterior and propose a method for theagent to approximately reach the optimal Bayesian result.This approach lends itself readily to processing over networksof agents.The paper is organized as follows. In the next section westate the problem. In Section 3, we present the mathematicaldevelopment of the proposed method. With the result inSection 3, we propose the method for estimation in Section 4.Section 5 provides simulation results, and Section 6 containsconclusions.2. PROBLEM FORMULATIONWe address sequential Bayesian estimation of a vector oflinear parameters θ ∈ RK×1 by a single agent. Thesolution of this problem is a key to resolving the problemof distributed Bayesian estimation, where each agent receivesmeasurements from its sensors, with the measurements beingdistorted by time independent but spatially correlated noise.In Fig. 1, we display a general scenario of distributedestimation, where each agent has M sensors providing theagent with spatially correlated observations. The agentsexchange information for cooperative estimation of theunknowns of interest.Fig. 1. A schematic diagram of a distributed system, whereAn denotes the nth agent, S(m)n represents the mth sensor ofagent An, and y(m)n,t denotes the measurements provided fromS(m)n to An at time instant t.Consider that at each time instant t ∈ N+, an agent receivesdata yt generated by a linear model of the formyt = Htθ +wt, (1)where Ht ∈ RM×K is a matrix known by the agent, and wtdenotes the observation noise modeled as a Gaussian randomvector with zero-mean and covariance Σ ∈ RM×M , which istime invariant but unknown. It is also assumed that M ≥ Kand the observation noise wt is white in time.We use the model in (1) and deﬁne the agent’s prior on Σ asan inverse Wishart distribution, W−1(Λ0, ν), given byp(Σ) ∝ |Σ|− ν+M+12 exp(−12tr(Λ0Σ−1)), (2)where Λ0 is a scale matrix, ν represents degrees offreedom, and tr(·) denotes the trace of the matrix inside theparentheses. Let the prior of θ be Gaussian denoted byp(θ|Σ) ∝ |C−10 |12 exp(−12(θ − θ0)C−10 (θ − θ0)),(3)where  denotes matrix transpose. We set C−10 =H0 Σ−1H0.From the Bayes’ rule, the posteriors of θ and Σ can beobtained byp (θ,Σ|It) ∝ p (It|θ,Σ) p (θ|Σ) p (Σ) , (4)where It refers to all the information up to time t, whichincludes yτ , Hτ , for all τ ∈ {1, 2, · · · , t}. We propose thatthe agent formulates its belief about θ as the marginalizedposterior of θ given by,βt = p(θ|It)=∫Σ0p(θ,Σ|It)dΣ, (5)where Σ  0 denotes that Σ is positive deﬁnite.We also deﬁne a benchmark for this problem. We assumethat there exists a genie agent which knows both It and Σ.Its belief is given by,β(opt)t = p(θ|It,Σ). (6)In this work, our aim is to derive expressions for sequentialupdate of the posterior βt and the estimate of θ.3. MATHEMATICAL DEVELOPMENT OF THEPROPOSED METHODIn this section, we ﬁrst derive the posterior held by the agentand then we propose the approximations to make the problemtractable for distributed estimation of θ.3.1. Posterior of individual agentUsing the model in Section 2, the expression (1) suggests thatthe likelihood of the data at time t can be written byp(It|θ,Σ) ∝ |Σ|− t2 exp(−12t∑τ=1wτ Σ−1wτ)(7)with wτ = yτ −Hτθ denoting the observation noise.Using the Bayes’ rule, one can show that the posterior of θ isgiven byp(θ,Σ|It) ∝ p(It|θ,Σ)p(θ|Σ)p(Σ)= |Σ|− ν+t+M+12 |C0|− 12 exp(−12Lt), (8)where Lt is deﬁned byLt =t∑τ=1wτ Σ−1wτ + tr(Λ0Σ−1)+(θ − θ0)C−10 (θ − θ0). (9)We show in the appendix that Lt can be reformulated into thefollowing quadratic form:Lt = (θ − θ˜t)C−1t (θ − θ˜t) + St (10)with Ct and θ˜t given byCt =(t∑τ=0Hτ Σ−1Hτ)−1, (11)θ˜t = C−1t(t∑τ=1Hτ Σ−1yτ +C−10 θ0), (12)where St is deﬁned bySt = (θ0 − μt) (C0 +Mt)−1 (θ0 − μt) + tr(Λ0Σ−1)+t∑τ=1(yτ −Hτμt)Σ−1(yτ −Hτμt). (13)In (13), μt represents the maximum likelihood estimate of θin the form ofμt =(t∑τ=1Hτ Σ−1Hτ)−1( t∑τ=1Hτ Σ−1yτ), (14)and Mt denotes the covariance matrix of μt deﬁned by,Mt =(t∑τ=1Hτ Σ−1Hτ)−1. (15)We write p(θ,Σ|It) = p(θ|Σ, It)p(Σ|It), and from theexpression of Lt in (10), we can get the optimal belief held bythe genie agent β(opt)t = p(θ|Σ, It). The belief is expressedby a multivariate Gaussian distribution given byβ(opt)t = N(θ˜t,Ct). (16)By using (8) and (10), we can show that after integrating outθ, the marginalized posterior of Σ has the following form:p(Σ|It) ∝ |C0C−1t |−12 |Σ|− ν+t+M+12 exp(−12St), (17)where St is deﬁned in (13).3.2. Approximation of the posteriorThe expression for μt in (14) suggests that μt is a functionof Σ, and thus, the posterior of Σ in (17) is not an inverseWishart distribution. This means that the integrationp(θ|It) =∫Σ0p(θ|Σ, It)p(Σ|It)dΣ, (18)is not tractable.However, from (17) and (13), Σ has an inverse Wishartdistribution if the following three approximations are valid.First,(C0 +Mt)−1 ≈ C−10 , (19)which is valid for large t. The reason is that when t is large,the elements ofMt become very small in comparison to thoseof C0. Second,|C0C−1t | ≈ rt, ∀t ∈ N, (20)where rt is approximately a constant. In the special case whentheHts are identical or proportional to each other for differentts, rt is truly a constant.The third approximation is that when t is large, the estimateof Σ held by an agent becomes close to the true value of Σ.Then, in (14), we substitute Σ with Σ̂t−1, and we have theapproximation μt ≈ μ̂t whereμ̂t =(t∑τ=1Hτ Σ̂−1t−1Hτ)−1( t∑τ=1Hτ Σ̂−1t−1yτ), (21)which means the agent can be viewed as if it knows the truevalue of Σ in calculating μt.Even with these approximations, the joint posterior of θ andΣ is not a Normal Inverse Wishart distribution. Then wepropose that the agent uses an additional approximation forgenerating its pseudo-posterior, i.e.,p(θ|It) ≈ β̂t, whereβ̂t = p(θ|It,Σ = Σ̂t), (22)and where Σ̂t denotes the latest estimate of Σ held by theagent at time instant t, which can be, e.g., be the maximumlikelihood estimate. However, in this paper, we propose thatthe agent uses the MMSE estimate. This estimate is the meanof the marginalized posterior of Σ, i.e.,Σ̂MMSE,t =∫Σ0Σ p(Σ|It) dΣ. (23)With the above three assumptions, we can analytically solvethe integration in (23) and obtain Σ̂t byΣ̂t =1νt −M − 1( t∑τ=1(yτ −Hτ μ̂t)(yτ −Hτ μ̂t)+Λ0 +H0(θ0 − μ̂t)(θ0 − μ̂t)H0), (24)where νt = ν+ t and μ̂t is deﬁned in (21). Furthermore, thisestimate approaches Σ̂MMSE,t with time.The expressions (16) and (22) show that the pseudo posteriorof agent An of θ has the following form:β̂t = N(θ̂t, Ĉt), (25)whereĈt =(t∑τ=0Hτ Σ̂−1t Hτ)−1, (26)θ̂t = Ĉ−1t(t∑τ=1Hτ Σ̂−1t yτ + Ĉ−10 θ0), (27)with Σ̂t being deﬁned in (24).4. THE PROPOSED METHODIn this section, we summarize the method for updating thepseudo posterior of θ employed by the agent. At each timeinstant t, it implements the following steps:Initialization: At t = 0, the agent forms its prior by (2) and(3), and initializes Σ̂0 by Σ̂0 = Λ0ν−M−1 . The steps belowdescribe the tth recursion.Step 1 The agent receives the data yt and the regressors Htfrom the sensors and calculates μ̂t by (21).Step 2 With μ̂t and (24), the agent updates its estimate Σ̂t.Step 3 The agent forms its pseudo posterior β̂t as aGaussian distribution with a mean θ̂t and covariance Σ̂tdeﬁned in (26) and (27).We point out that since an agent approximates its belief by aGaussian distribution, the distributed estimation problem canbe converted to a distributed summation problem [6].5. SIMULATIONIn this section, we provide computer simulations that showthe performance of our method in terms of convergence andnumerical comparisons with the least squares (LS) method.We performed two experiments. In the ﬁrst experiment, weimplemented the proposed and the LS methods with identicalrandom data in 1000 realizations. In each of the trials, θ =[3, 3, 2, 2], M = 10, K = 4, t ∈ {1, 2, · · · , 200}. Weset the elements of Ht ∈ RM×K to be independent randomvariables uniformly distributed on [3, 5]. Also, in every trial,we drew Σ from its prior, an inverse Wishart distributionwith Λ0 = 10IM (with IM ∈ RM×M denoting the identitymatrix) and ν = 12. We also set θ0 = [0, 0, 0, 0] andgenerated H0 in the same way as we did Ht, t > 0.As a performance metric for the different methods, we usedthe mean square deviation at time t, MSD(t), deﬁned asthe average value of ‖θ̂t − θ‖2 over 1000 implementations.To demonstrate the advantage of the proposed method, wecompared the MSD of the proposed method with that of theLS method and of the genie agent. The LS estimate of θ attime instant t is given byθ̂(LS)t =(t∑τ=1Hτ Hτ)−1( t∑τ=1Hτ yτ). (28)The MSD of the genie agent served as a benchmark. We notethat the estimate of θ of the genie agent is θ˜t given by (12).In the second experiment, we repeated everything except thatwe generated Σ with another scale matrix Q2, where Q ∈RM×M was with elements that were independent randomvariables uniformly distributed on [0, 5]. The intention ofthis experiment was to show the performance of the proposedmethod with highly correlated data.0 100 20010−210−1100101102tExperiment 20 100 20010−1100101tMSDtExperiment 1LS methodProposed methodFictitious genie agentFig. 2. Asymptotical performance of the proposed method.The results of the two experiments are shown in Fig. 2 (onthe left of experiment 1, and on the right of experiment 2),where we plotted the MSDs of the different methods. It canbe seen that the proposed method shows a faster convergencethan the LS method in terms of MSD. This difference is evenmore obvious when processing data with a higher correlation(as in experiment two).0 50 100 150 200100102104106108tDKLExperiment 1Experiment 2Fig. 3. The evolution of KL divergence.In Fig. 3, we plotted the Kullback–Leibler (KL) divergence,DKL(β(opt)t ||β̂t)), between the belief of the genie agent andthe belief of the agent that employs the proposed method.From the ﬁgure, it can be seen that the KL divergence keepsdecreasing as time evolves, but it does it rather slowly.6. CONCLUSIONIn this paper, we considered Bayesian estimation in thepresence of Gaussian noise with unknown covariance butthat is independent in time. We ﬁrst derived the posteriorof the agent and then based on that we proposed threeapproximations that allow for a closed form solution for thebelief update. We also presented an approach for recursiveestimation of the parameter vector of interest. This resultcan be used in parameter estimation over networks. Bycomputer simulations, we showed that the proposed methodoutperformed the least squares method in terms of meansquare deviation.7. APPENDIXHere we show the derivation from (9) to (15). First byexpanding the square terms in (9), we haveLt =t∑τ=0θHτ Σ−1Hτθ+t∑τ=1yτ Σ−1yτ + θ0 C−10 θ0 + tr(Λ0Σ−1)−2θ(t∑τ=1Hτ Σ−1yτ +C−10 θ0). (29)From (12), we have C−1t θ̂t =∑tτ=1Hτ Σ−1yτ + C−10 θ0,which implies that if we add and subtract θ̂t C−1t θ̂t in theabove equation, it will becomeLt = (θ − θ̂t)C−1t (θ − θ̂t)− θ̂t C−1t θ̂t+t∑τ=1yτ Σ−1yτ + θ0 C−10 θ0 + tr(Λ0Σ−1). (30)With the above equation, we can writeSt =t∑τ=1yτ Σ−1yτ + θ0 C−10 θ0 + tr(Λ0Σ−1)−θ̂t C−1t θ̂t. (31)Next, let S(1)t =∑tτ=1 yτ Σ−1yτ − μt M−1t μt with Mtbeing deﬁned in (15). From (14) and (15), we haveμt M−1t μt = μtt∑τ=1Hτ Σ−1yτ , (32)which impliesS(1)t =t∑τ=1(yτ −Hτμτ )Σ−1(yτ −Hτμτ ). (33)Similarly, we can deﬁne S(2)t = θ0 C−10 θ0 + μt M−1t μt −θ̂t C−1t θ̂t, where the last term can be expanded asθ̂t C−1t θ̂t = (C−10 θ0 +M−1t μt)(C−10 +M−1t )−1×(C−10 θ0 +M−1t μt)= μt M−1t (C−10 +M−1t )−1M−1t μt+θ0 C−10 (C−10 +M−1t )−1C−10 θ0−2μt M−1t (C−10 +M−1t )−1C−10 θ0. (34)With the result in (34), one can show that S(2)t can bereformulated asS(2)t = (θ0 − μt)(C−10 (C−10 +M−1t )−1M−1t)× (θ0 − μt)= (θ0 − μt) (C0 +Mt)−1 (θ0 − μt) . (35)Noting that St = S(1)t +S(2)t + tr(Λ0Σ−1), by (33) and (35)we have shown that the equation (10) holds.REFERENCES[1] I. D. Schizas, A. Ribeiro, and G. B. Giannakis, “Consensusin ad hoc wsns with noisy links part I: Distributed estimationof deterministic signals,” IEEE Transactions on SignalProcessing, vol. 56, no. 1, pp. 350–364, 2008.[2] F. S. Cattivelli and A. H. Sayed, “Diffusion LMS strategiesfor distributed estimation,” IEEE Transactions on SignalProcessing, vol. 58, no. 3, pp. 1035–1048, 2010.[3] P. M. Djuric´ and Y. Wang, “Distributed Bayesian learningin multiagent systems: Improving our understanding ofits capabilities and limitations,” IEEE Signal ProcessingMagazine, vol. 29, no. 2, pp. 65–76, 2012.[4] Y. Wang and P. M. Djuric´, “Sequential estimation of linearmodels in distributed settings,” in Proceedings of the EuropeanSignal Processing Conference, 2013.[5] Y. Wang and P. M. Djuric´, “Reaching bayesian belief overnetworks in the presence of communication noise,” in GlobalConference on Signal and Information Processing (GlobalSIP),2013 IEEE. IEEE, 2013, pp. 591–594.[6] L. Xiao, S. Boyd, and S.J. Kim, “Distributed average consensuswith least-mean-square deviation,” Journal of Parallel andDistributed Computing, vol. 67, no. 1, pp. 33–46, 2007.[7] R. Olfati-Saber and R. M. Murray, “Consensus problems innetworks of agents with switching topology and time-delays,”IEEE Transactions on Automatic Control, vol. 49, no. 9, pp.1520–1533, 2004.[8] Z. Weng and P. M. Djuric´, “Efﬁcient estimation of linearparameters from correlated node measurements over networks,”Signal Processing Letters, IEEE (to appear), 2014.

On Sequential Estimation of Linear Models From  Data with Correlated Noise

https://zenodo.org/record/54521/files/1569926831.pdf

On Sequential Estimation of Linear Models From Data with Correlated Noise

Abstract

Similar works

Full text

Available Versions

ZENODO