<figure><img src="image_1.png" /></figure>

HAL
open science

# Interacting Kalman Filters for Linear Systems with Coupling Based on Empirical Covariances

Aneel Tanwani

► To cite this version:

Aneel Tanwani. Interacting Kalman Filters for Linear Systems with Coupling Based on Empirical Covariances. IEEE Conference on Decision and Control, Dec 2025, Rio de Janeiro, Brazil. hal-05394048

HAL Id: hal-05394048

https://laas.hal.science/hal-05394048v1

Submitted on 2 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Interacting Kalman Filters for Linear Systems with Coupling Based on Empirical Covariances

Aneel Tanwani

**Abstract**— The problem of designing estimators for stochastic linear systems with distributed observations is considered. Each observation process is associated to a node in an undirected graph, which is used to compute a local estimate at the node. The injection gain used in the filter at each node is obtained from the empirical covariance of all the estimates available at that node. This way the underlying idea comes from the theory of ensemble filtering and we analyze the evolution of the coupled covariances over the entire graph. After providing a detailed derivation of the evolution of covariance matrix, we observe that the drift term in the differential equation for the coupled covariances has some inherent stability structure in case of regular graphs, which leads to the fluctuations around the steady state. We provide an illustration of our algorithm on an academic example while comparing it with centralized and ensemble Kalman filters.

***Index Terms***—Distributed filtering; Ensemble filtering; Interacting filters.

## I. INTRODUCTION

Design problems over networks have important applications, and there is an interest to design algorithms for such problems which are computationally efficient and improve the performance of the overall system while taking the exchange of information over the network [1]. Such criteria are also important when studying the filtering problem with distributed sensors, where the objective is to integrate observations from different sources to enhance the performance by minimizing the mean square estimation error. This paper considers the problem of state estimation for stochastic systems in distributed setup, where our primary objective is to propose a computationally light filtering algorithm and analyze its performance compared to some known techniques.

Most of the literature on distributed filtering appears in the setting of discrete-time systems. Traditional centralized approaches [2], [3] are limited in scalability and robustness, leading researchers to explore decentralized and distributed frameworks [4]. Earlier results on distributed filtering date back to [5] which consider the problem of estimating a random variable and also provide conditions on the formation of sensors that allow asymptotic convergence of the estimates. Thereafter, different techniques have been adopted in the literature and in particular, there has been some focus on combining Kalman filters with appropriate

consensus algorithms. The papers [6], [7] provided some initial developments on consensus-based Kalman filters where sensor nodes exchange state estimates with neighbors to collectively refine their predictions. This idea was then generalized further in [8], where the authors take several parameters into consideration such as the consensus matrix, filtering gain, and the rate of communication for reaching consensus. The idea of using time-varying consensus weights in distributed filtering has been explored in [9], whereas the detectability and stability of error dynamics has been studied in [10]. Distributed filtering with more general probability distributions and application to Gaussian distributions with linear discrete-time systems is studied in [11]. An instance of this problem with randomly switching graphs (or consensus weights) was studied by the author in [12]. More recently, we have also seen the distributed optimization viewpoint [13] which reformulates the optimal estimation problem as a problem of distributed optimization [14].

Another significant research direction in filtering algorithms, which also involves combining different Kalman filters, is the class of ensemble filtering algorithms. These methods focus on designing suboptimal filters that require less computational effort compared to their optimal counterparts, and are partly inspired by Monte Carlo integration methods [15]. One of the earlier works introducing ensemble Kalman filters [16] motivated the use of this technique for large-scale applications in geophysical sciences. The method involves simulating the evolution of individual particles in an ensemble through differential equations, which are coupled using the empirical mean and empirical error covariance. Several reviews [17], [18] and books [19], [20] provide an overview of the progress in this field. However, most of these works provide limited details on the theoretical analysis of the proposed filtering techniques, and rigorous mathematical investigation of ensemble Kalman filters has only recently gained traction [21], [22], [23]. Ensemble filtering methods have been analyzed through the lens of mean-field models, which are described by stochastic differential equations. In this framework, the ensemble particles serve as approximations of these mean-field models. Recent review articles elaborating on this perspective include [23], [24]. The limiting behavior of these particles is governed by a McKean-Vlasov-type diffusion process, often referred to as the mean-field process. In the literature, different formulations of this mean-field process have been proposed – for instance, by introducing noise into the prediction and correction steps of the Kalman-Bucy process [22], or by

A. Tanwani is affiliated with CNRS - LAAS, University of Toulouse, CNRS, 31400 Toulouse, France (e-mail: aneel.tanwani@laas.fr).

This work is sponsored by the project CYPHAI, financed by ANR-JST CREST program with grant number ANR-20-JSTM-0001.
defining it as a non-diffusion equation that is optimal in the measure transportation sense [25].

Our focus in this article is on designing a computationally efficient filtering algorithm with observations distributed over a network. In doing so, we borrow ideas from the literature on ensemble filtering so that the problem of distributed filtering can be solved with lighter computations. We consider undirected graphs where the nodes have access to an independent observation process. To keep the filtering algorithm light, the estimate is computed without solving the Riccati differential equation at any of the nodes. Each node therefore computes a suboptimal estimate of the state using the empirical covariance obtained from the estimates of the neighboring nodes. The resulting algorithms do not have the symmetric structure that makes the analysis tractable in ensemble filtering. We provide some initial investigations into studying the evolution of empirical variances in the distributed setup. It is also shown that in case of regular graphs (where each node has the same number of neighbors), we see certain invariance structure in the evolution of covariances which is important for analyzing stability later on.

## II. FILTERING BASED ON INTERACTIONS : AN OVERVIEW

We consider dynamical systems modeled by linear stochastic differential equations of the form

$$dx_t = Ax_t dt + G d\omega_t \quad (1)$$

where $(x_t)_{t \ge 0}$ is an $\mathbb{R}^n$-valued diffusion process describing the state. Let $(\Omega, \mathcal{F}, P)$ denote the underlying probability space. It is assumed that, for each $t \ge 0$, $(\omega_t)_{t \ge 0}$ is a zero mean $\mathbb{R}^m$-valued standard Wiener process with the property that $E[d\omega_t d\omega_t^\top] = I_m dt$, for each $t \ge 0$. The matrices $A \in \mathbb{R}^{n \times n}$ and $G \in \mathbb{R}^{n \times m}$ are taken as constant, and the process $(\omega_t)_{t \ge 0}$ does not depend on the state. The solution of (1) is interpreted in the sense of [26, Theorem 5.2.1], and is a stochastic process $(x_t)_{t \ge 0}$ which is adapted to a certain filtration associated with the initial condition and the Wiener process $(\omega_t)_{t \ge 0}$.

The centralized output measurement associated to the process (1) is of the form

$$dy_t = Hx_t dt + d\nu_t \quad (2)$$

where $H \in \mathbb{R}^{p \times n}$ is a constant matrix, with $(A, H)$ being observable, and $(\nu_t)_{t \ge 0}$ is a zero mean $\mathbb{R}^p$-valued standard Wiener process. The conventional filtering problem, with initial state having Gaussian distribution, deals with constructing a mean-square estimate of the state $x_t$, denoted by $\hat{x}_t$ so that $E[|x_t - \hat{x}_t|^2 | (dy(s))_{s \le t}]$ is minimized. The optimal estimate which achieves this minimum value is $E[x_t | (dy(s))_{s \le t}]$ and is computed recursively using a Kalman-Bucy filter.

### A. Vanilla ensemble Kalman filter

The classical solution to the filtering problem for system class (1) with observation process (2) requires us to solve a

Riccati differential equation to compute the injection gains used in the differential equation for the estimate. For a state process evolving in $\mathbb{R}^n$, this optimal filter involves solving $(n^2 + n)/2$ differential equations for the covariance process $(P_t)_{t \ge 0}$ and $n$ differential equations for the first moment of the estimate $(\hat{X}_t)_{t \ge 0}$. This can be quite cumbersome for large values of $n$ especially when we take into consideration the additional operations involved in computing these processes. For this reason, there has been extensive research for other methods to implement optimal filter, or its approximation. One such technique is based on the use of the ensemble, or feedback particle, filters. The basic idea of the particle filters is to simulate a collection of particles through stochastic differential equations which are coupled to each other through joint statistics of the population.

For the evolution of particles in the ensemble, we first describe a diffusion process of McKean-Vlasov type differential equation.

$$dS_t := AS_t dt + G d\bar{\omega}_t + Q_t H^\top V^{-1} [dy_t - (HS_t dt + d\bar{\nu}_t)] \quad (3a)$$

$$Q_t := E[(S_t - \hat{S}_t)(S_t - \hat{S}_t)^\top | (dy(s))_{s \le t}] \quad (3b)$$

$$\hat{S}_t := E[S_t | (dy(s))_{s \le t}], \quad (3c)$$

where the processes $\bar{\omega}_t$ and $\bar{\nu}_t$ are independent copies of the processes $\omega_t$ and $\nu_t$ respectively, that were introduced in (1) and (2). The motivation for working with this particular diffusion process comes from the fact that the conditional mean $\hat{S}_t$ and the conditional covariance $Q_t$ satisfy the same differential equations as the ones obtained for their optimal counterparts, that is, it can be checked that

$$\begin{aligned} d\hat{S}_t &= A\hat{S}_t dt + Q_t H^\top V^{-1} (dy_t - H\hat{S}_t dt) \\ \dot{Q}_t &= AQ_t + Q_t A^\top + GG^\top - Q_t H^\top V^{-1} H Q_t. \end{aligned}$$

From the mean-field process (3), we describe an ensemble of $M$ particles. For $\ell = 1, \dots, M$, the particle $S_t^\ell$ is described by the following stochastic differential equation conditioned upon the observation process :

$$dS_t^\ell = AS_t^\ell dt + Gd\omega_t^\ell + K_t^M (dy_t - HS_t^\ell dt - d\nu_t^\ell) \quad (4)$$

where, for each $\ell$, $\omega_t^\ell$ and $\nu_t^\ell$ represent the independent copies of the processes $\omega_t$ and $\nu_t$ respectively. The term $K_t^M$ denotes the empirical gain computed as follows:

$$K_t^M = Q_t^M H^\top V^{-1} \quad (5)$$

in which $Q_t^M$ denotes the empirical covariance, which is computed using the empirical mean $\hat{S}_t^M$, that is,

$$Q_t^M = \frac{1}{M} \sum_{\ell=1}^{M} (S_t^\ell - \hat{S}_t^M)(S_t^\ell - \hat{S}_t^M)^\top \quad (6)$$

$$\hat{S}_t^M = \frac{1}{M} \sum_{\ell=1}^{M} S_t^\ell. \quad (7)$$

One clearly sees that the computational effort in implementing the ensemble filters (4) involves solving $Mn$ differential equations.
## B. Distributed setup

When looking at the filtering problem in distributed setup, the measurements associated with system (1) are obtained from a set of $N$ sensors which are distributed in their localization and denoted by $dy_i$. For the sake of simplicity, in this article, we take

$$dy_i = Hx dt + d\nu_i, \quad i = 1, \dots, N. \qquad (8)$$

That is, for each node, $(dy_i(t))_{t \ge 0}$ describes an $\mathbb{R}^p$-valued continuous-time observation process. In equation (8), $\nu_i$ is a standard Wiener process, taking values in $\mathbb{R}^p$, and $E[d\nu_i(t) d\nu_i(t)^\top] = V dt$, with $V \in \mathbb{R}^{p \times p}$ assumed to be positive definite. Note that, even though the matrix $H$ is the same for each observation process and $\nu_i$ have the same distribution for each $i = 1, \dots, N$, the processes $d\nu_i$ and hence, $dy_j$, for $i \neq j$, are independent from each other.

The sensor nodes are connected via a graph $G = (\mathcal{V}, \mathcal{E})$, where $\mathcal{V} = \{1, \dots, N\}$ is the set of graph nodes, and $\mathcal{E}$ contains all the edges defined by a subset of the pairs $(i, j)$, $i \neq j$, $i, j \in \mathcal{V}$. We assume that the graph is undirected and connected. The neighbors of a node $i \in \mathcal{V}$ are denoted by $N_i$ and we adopt the convention that $i \in N_i$. The cardinality of the set $N_i$ is called the degree of the node $i$, and is denoted $d_i$. If each node in the graph has the same degree, we call it a regular graph. If $j \in N_i$, it simply means that node $i$ and $j$ can exchange information with each about their own local estimates. With this information structure in place, we are interested in the following problem:

$$\min_{\hat{x}_i} E[|x(t) - \hat{x}_i(t)|^2 | (dy_i(s), x_j(s)) \text{ with } s \le t, j \in N_i].$$

## C. Discussion about possible solutions

With the observation processes given in (8), each sensor node can compute an optimal state estimate $\hat{x}_i$ which minimizes the mean square error $E[|x(t) - \hat{x}_i(t)|^2 | (dy_i(s))_{s \le t}]$ without any exchange of information from the other nodes in the graph. However, due to exchange of information from other nodes, we expect each node to improve their estimate. Let us consider the following example to see that the increase in number of observations of a hidden random variable decreases the variance of the estimation error.

**Example 1.** Let $x$ denote some unknown constant, and $y_i$ provides a noisy measurement of $x$, that is, $y_i = x + \beta_i$ where $\beta_i$ is normally distributed with mean zero and variance $\sigma_\beta$. If we use one single measurement, the error variance in the estimate is $\sigma_\beta$. With $N$ such measurements, the optimal mean square estimate would be to take the mean of all the measurements that is $\hat{x} = \frac{1}{N} \sum_{i=1}^N y_i$, and the error variance decreases by a factor $\frac{1}{N}$.

Going back to the distributed filtering problem setup, the above example suggests that if the graph connections allow, the optimal solution may be obtained by collecting all the information at individual nodes. This is possible for complete

graphs where each node is connected to every other node in the graph. The resulting solution in such a case resembles the centralized solution. However, in general, with distributed sensor nodes, the nodes can only compute the estimate by exchanging information with their neighboring nodes, and they do not necessarily have the complete knowledge of the observations available at all the other nodes in the graph.

**Comparisons with ensemble filter:** In ensemble filtering algorithm (4), each particle updates its estimate using the empirical covariance which is computed from the knowledge of all the other particles in the ensemble. However, there is only one observation process which is accessible to all the particles. One can look at this algorithm as a distributed filtering algorithm over complete graphs (with one slight subtlety that the same observation process is shared by all the particles in ensemble filtering, whereas equation (8) describes $N$ independent processes, one for each node). From that viewpoint, the ensemble filtering algorithm actually provides a coupling rule for updating the estimate at each node using empirical covariance of the available population, and this particular coupling rule has not appeared in the literature on distributed filtering with generic graph structure. In next section, we investigate the utility of this empirical coupling rule for general graphs.

## III. ANALYSIS WITH EMPIRICAL COUPLINGS

Motivated by the discussions in the previous section, our goal is to study the problem of distributed filtering for connected and undirected graphs, where the individual nodes implement an estimation algorithm using the empirical covariance of the population comprising only its neighbors. This way the computational burden on each node is significantly reduced as the agents no longer have to solve the Riccati differential equations. Compared to the setting of ensemble filter (4), the difference arises because each node has possibly different neighbors and therefore the empirical covariance computed at each node is different. In other words, the injection gains used by each node are different whereas they are identical in (4). Due to this break in symmetry, we can no longer work with the arguments used in the analysis of usual ensemble filters, and need to investigate the new setup more carefully.

### A. Filtering algorithm

The algorithm that we propose for distributed filtering is thus described as follows: At each node $i \in \mathcal{V}$, we simulate the process $\hat{x}_i$ described by the differential equation

$$d\hat{x}_i = A\hat{x}_i dt + Gd\omega^i + \hat{P}_i H^\top V^{-1}(dy_i - H\hat{x}_i dt). \quad (9)$$

Here, $\omega^i$ is an independent copy of the process noise $\omega$, and $\hat{P}_i$ is the empirical covariance computed by the $i$-th node based on the values of estimates computed by node $i$ and its neighbors, that is,

$$\hat{P}_i = \frac{1}{d_i} \sum_{j \in N_i} (\hat{x}_j - \bar{x}_i)(\hat{x}_j - \bar{x}_i)^\top,$$
where, $\bar{x}_i$ denotes the mean of the estimates that are available at i-th node, that is,

$$ \bar{x}_i = \frac{1}{d_i} \sum_{j \in \mathcal{N}_i} \hat{x}_j. $$

In what follows, we will analyze how this algorithm behaves over time and whether it provides meaningful estimates.

**Remark 1.** Note that, in the conventional ensemble Kalman filter, the innovation term contains an independent copy of the measurement noise. This correction term is necessary because the particle dynamics are driven by a common observation process and hence the noise in the observation process cancels out when we look at $S^\ell - \hat{S}^M$, for $\ell \in \{1, \dots, M\}$. In (9), we do not write such an additional measurement noise because the observation process at each node is independent of each other and this noise does not cancel out in the equation for $\hat{x}_i - \bar{x}_i$, for $i \in \mathcal{V}$.

## B. Evolution of local mean and empirical covariance

To analyze the estimation algorithm proposed in (9), we start by writing the differential equation for the evolution of the local mean $\bar{x}_i$. We observe that, for each $i \in \mathcal{V}$, with its neighboring nodes denoted by $\mathcal{N}_i$, we have

$$ d\bar{x}_i = A\bar{x}_i dt + \frac{1}{\sqrt{d_i}} G d\tilde{\omega}^i + \frac{1}{d_i} \sum_{j \in \mathcal{N}_i} \hat{P}_j H^\top V^{-1} (dy_j - H\hat{x}_j dt) \quad (10) $$

where we used the notation

$$ \tilde{\omega}^i := \frac{1}{\sqrt{d_i}} \sum_{j \in \mathcal{N}_i} \omega^j \quad \text{so that} \quad E[d\tilde{\omega}_t^i d\tilde{\omega}_t^{i\top}] = I_m dt. $$

Next, for each $i \in \mathcal{V}$, and $j \in \mathcal{N}_i$, we consider the deviation of $\hat{x}_j$ from $\bar{x}_i$, which is denoted by $q_{ij} := \hat{x}_j - \bar{x}_i$. We compute its differential as follows:

$$ d q_{ij} = A q_{ij} dt + \hat{P}_j H^\top V^{-1} (dy_j - H \hat{x}_j dt) + G d \omega^j - \frac{1}{d_i} \sum_{k \in \mathcal{N}_i} \hat{P}_k H^\top V^{-1} (dy_k - H \hat{x}_k dt) - \frac{1}{\sqrt{d_i}} G d \tilde{\omega}^i. \quad (11) $$

To get the evolution equation for $\hat{P}_i$, we next compute $d(q_{ij}q_{ij}^\top)$. Using Itô's Lemma, we have

$$ d(q_{ij}q_{ij}^\top) = d q_{ij}q_{ij}^\top + q_{ij}d q_{ij}^\top + d q_{ij}d q_{ij}^\top. \quad (12) $$

Let us use the notation

$$ \bar{H} := H^\top V^{-1} H, \quad \text{and} \quad L_j := \hat{P}_j H^\top V^{-1} $$

so that $\bar{H}$ is a symmetric matrix. Using (11) and Itô rule for multiplication, the last term in (12) equates to

$$ d q_{ij} d q_{ij}^{\top} = \left(1 - \frac{1}{d_i}\right)^2 \hat{P}_j \bar{H} \hat{P}_j dt + \left(1 - \frac{1}{d_i}\right) G G^{\top} dt + \frac{1}{d_i^2} \sum_{k \in \mathcal{N}_i \setminus \{j\}} \hat{P}_k \bar{H} \hat{P}_k dt. \quad (13) $$

Plugging the expressions from (11) and (13) in (12), we get

$$ \begin{aligned} d(q_{ij}q_{ij}^{\top}) = & A(q_{ij}q_{ij}^{\top})dt + (q_{ij}q_{ij}^{\top})A^{\top}dt + \hat{P}_j\bar{H}(x-\hat{x}_j)q_{ij}^{\top}dt \\ & + q_{ij}^{\top}(x-\hat{x}_j)^{\top}\bar{H}\hat{P}_jdt - \frac{1}{d_i}\sum_{k \in \mathcal{N}_i} \hat{P}_k\bar{H}(x-\hat{x}_k)q_{ij}^{\top}dt \\ & - \frac{1}{d_i}\sum_{k \in \mathcal{N}_i} q_{ij}(x-\hat{x}_k)^{\top}\bar{H}\hat{P}_kdt \\ & + \left(1-\frac{1}{d_i}\right)^2 \hat{P}_j\bar{H}\hat{P}_jdt + \left(1-\frac{1}{d_i}\right) GG^{\top}dt \\ & + \frac{1}{d_i^2}\sum_{k \in \mathcal{N}_i \setminus \{j\}} \hat{P}_k\bar{H}\hat{P}_kdt + \mathcal{M}_{ij} + \mathcal{M}_{ij}^{\top} \end{aligned} \quad (14) $$

In the last equation, $\mathcal{M}_{ij}$ contains all the noise terms:

$$ \mathcal{M}_{ij} := \left( L_j d\nu_j - \frac{1}{d_i} \sum_{k \in \mathcal{N}_i} L_k d\nu_k \right) q_{ij}^{\top} + G \left( d\omega^j - \frac{1}{\sqrt{d_i}} d\tilde{\omega}^i \right) q_{ij}^{\top}. \quad (15) $$

We now write the equation for $d\hat{P}_i = \frac{1}{d_i} \sum_{j \in \mathcal{N}_i} d(q_{ij}q_{ij}^\top)$, while noting that,

$$ \begin{aligned} & \frac{1}{d_i} \sum_{j \in \mathcal{N}_i} \sum_{k \in \mathcal{N}_i} \hat{P}_k \bar{H}(x - \hat{x}_k) q_{ij}^{\top} = 0 \\ & \frac{1}{d_i} \sum_{j \in \mathcal{N}_i} \sum_{k \in \mathcal{N}_i} L_k d\nu_k q_{ij}^{\top} = 0 \end{aligned} $$

which yields

$$ \begin{aligned} d\hat{P}_i = & A\hat{P}_i dt + \hat{P}_i A^\top dt + \frac{1}{d_i} \sum_{j \in \mathcal{N}_i} \hat{P}_j \bar{H}(x - \hat{x}_j) q_{ij}^\top dt \\ & + \frac{1}{d_i} \sum_{j \in \mathcal{N}_i} q_{ij}(x - \hat{x}_j)^\top \bar{H} \hat{P}_j dt \\ & + \frac{1}{d_i} \left(1 - \frac{1}{d_i}\right) \sum_{j \in \mathcal{N}_i} \hat{P}_j \bar{H} \hat{P}_j dt \\ & + \left(1 - \frac{1}{d_i}\right) GG^\top dt + \bar{M}_{ij} + \bar{M}_{ij}^\top \end{aligned} \quad (16) $$

where we used the notation

$$ \bar{M}_{ij} := \frac{1}{d_i} \sum_{j \in \mathcal{N}_i} \left( L_j d\nu_j + G d\omega^j - \frac{1}{\sqrt{d_i}} G d\tilde{\omega}^i \right) q_{ij}^{\top}. $$

By writing,

$$ \hat{P}_j \bar{H}(x - \hat{x}_j) = - \hat{P}_j \bar{H}(\hat{x}_j - \bar{x}_i) + \hat{P}_j \bar{H}(x - \bar{x}_i), \quad (17) $$

we get

$$ \sum_{j \in \mathcal{N}_i} \hat{P}_j \bar{H}(x - \hat{x}_j) q_{ij}^{\top} = - \sum_{j \in \mathcal{N}_i} \hat{P}_j \bar{H} q_{ij} q_{ij}^{\top} + \sum_{j \in \mathcal{N}_i} \hat{P}_j \bar{H}(x - \bar{x}_i) q_{ij}^{\top}. \quad (18) $$

Substituting this last expression in (16), we obtain
$$
\begin{align}
d\hat{P}_i ={}& A\hat{P}_i dt + \hat{P}_i A^\top dt - \frac{1}{d_i} \sum_{j \in \mathcal{N}_i} \hat{P}_j \overline{H} q_{ij} q_{ij}^\top dt \nonumber \\
& - \frac{1}{d_i} \sum_{j \in \mathcal{N}_i} q_{ij} q_{ij}^\top \overline{H} P_j dt + \frac{1}{d_i} \sum_{j \in \mathcal{N}_i} \hat{P}_j \overline{H} (x - \bar{x}_i) q_{ij}^\top dt \nonumber \\
& + \frac{1}{d_i} \sum_{j \in \mathcal{N}_i} q_{ij} (x - \bar{x}_i)^\top \overline{H} \hat{P}_j dt \nonumber \\
& + \frac{1}{d_i} \left(1 - \frac{1}{d_i}\right) \sum_{j \in \mathcal{N}_i} \hat{P}_j \overline{H} \hat{P}_j dt \nonumber \\
& + \left(1 - \frac{1}{d_i}\right) GG^\top dt + \overline{M}_{ij} + \overline{M}_{ij}^\top. \tag{19}
\end{align}
$$

This completes the derivation for the evolution of the local covariance matrix $\hat{P}_i$ at each node $i \in \mathcal{V}$.

### C. Asymptotic invariance of $\hat{P}_i$

With equations (9) and (19) at hand, it is of interest to analyze how the processes $(\hat{x}_i, \hat{P}_i)$ evolve over time in our setup. Complete analysis of this sort requires further investigation, but at this moment, let us identify an important property of the covariance process $\hat{P}_i$ that can be obtained from (19).

For that purpose, let us consider the dynamics associated with the drift term in (19), which can be described by the following differential equation

$$
\begin{align}
d\tilde{P}_i ={}& A\tilde{P}_i dt + \tilde{P}_i A^\top dt - \frac{1}{d_i} \sum_{j \in \mathcal{N}_i} \tilde{P}_j \overline{H} q_{ij} q_{ij}^\top dt \nonumber \\
& - \frac{1}{d_i} \sum_{j \in \mathcal{N}_i} q_{ij} q_{ij}^\top \overline{H} \tilde{P}_j dt + \frac{1}{d_i} \sum_{j \in \mathcal{N}_i} \tilde{P}_j \overline{H} (x - \bar{x}_i) q_{ij}^\top dt \nonumber \\
& + \frac{1}{d_i} \sum_{j \in \mathcal{N}_i} q_{ij} (x - \bar{x}_i)^\top \overline{H} \tilde{P}_j dt \nonumber \\
& + \frac{1}{d_i} \left(1 - \frac{1}{d_i}\right) \sum_{j \in \mathcal{N}_i} \tilde{P}_j \overline{H} \tilde{P}_j dt + \left(1 - \frac{1}{d_i}\right) GG^\top dt \tag{20}
\end{align}
$$

with the initial condition $\tilde{P}_i(0) = \hat{P}_i(0)$, for each $i = 1, \dots, N$. Introduce the diagonal set $\Delta_P$ in $\mathbb{R}^{n_N \times n_N}$, such that

$$ \Delta_P := \{(P_1, \dots, P_N) \in \mathbb{R}^{n_N} \mid P_1 = P_2 = \dots = P_N\}. $$

On this set, the deterministic part of the dynamics for the local covariance matrix at each node indeed simplifies and we have the following result:

**Proposition 1.** Consider a regular graph with $N$ nodes, and for each $i = 1, \dots, N$, let the processes $\hat{P}_i$ be described by the differential equations (20). Then, the set $\Delta_P \in \mathbb{R}^{n_N}$ is forward invariant for the joint process $(\tilde{P}_1, \dots, \tilde{P}_N)$.

*Proof.* Let $(\tilde{P}_1, \dots, \tilde{P}_N) \in \Delta_P$, then for each $i = 1, \dots, N$, the dynamics (20) take the form

$$ d\tilde{P}_i = A\tilde{P}_i dt + \tilde{P}_i A^\top dt - \tilde{P}_i \overline{H} \frac{1}{d_i} \sum_{j \in \mathcal{N}_i} q_{ij} q_{ij}^\top dt $$

$$
\begin{align}
& - \frac{1}{d_i} \sum_{j \in \mathcal{N}_i} q_{ij} q_{ij}^\top \overline{H} \tilde{P}_i dt + \tilde{P}_i \overline{H} (x - \bar{x}_i) \frac{1}{d_i} \sum_{j \in \mathcal{N}_i} q_{ij}^\top dt \nonumber \\
& + \frac{1}{d_i} \sum_{j \in \mathcal{N}_i} q_{ij} (x - \bar{x}_i)^\top \overline{H} \tilde{P}_i dt \nonumber \\
& + \frac{1}{d_i} \left(1 - \frac{1}{d_i}\right) \sum_{j \in \mathcal{N}_i} \tilde{P}_i \overline{H} \tilde{P}_i dt + \left(1 - \frac{1}{d_i}\right) GG^\top dt. \tag{21}
\end{align}
$$

By using the definition of empirical covariance and noting that $\sum_{j \in \mathcal{N}_i} q_{ij} = 0$ for each $i$, we get

$$
\begin{align}
d\tilde{P}_i ={}& A\tilde{P}_i dt + \tilde{P}_i A^\top dt - \tilde{P}_i \overline{H} \tilde{P}_i dt - \tilde{P}_i \overline{H} \tilde{P}_i dt \nonumber \\
& + \left(1 - \frac{1}{d_i}\right) \tilde{P}_i \overline{H} \tilde{P}_i dt + \left(1 - \frac{1}{d_i}\right) GG^\top dt. \tag{22}
\end{align}
$$

This last equation is completely independent of the neighboring nodes, and is therefore identical for each $i = 1, \dots, N$ for regular graphs. Hence, with the same initial condition for each node $i \in \mathcal{V}$, the processes $\tilde{P}_i$ remain the same, proving that $\Delta_P$ is invariant. $\square$

In the proof of Proposition 1, equation (22) can be simply rewritten as

$$
\begin{align}
d\tilde{P}_i ={}& (A - \tilde{P}_i \overline{H}) \tilde{P}_i dt + \tilde{P}_i (A^\top - \overline{H} \tilde{P}_i) dt \nonumber \\
& + \left(1 - \frac{1}{d_i}\right) \tilde{P}_i \overline{H} \tilde{P}_i dt + \left(1 - \frac{1}{d_i}\right) GG^\top dt. \tag{23}
\end{align}
$$

In the literature on Riccati equations associated with the filtering problem, we see that the observability of the pair $(A, H)$ implies that the matrix $(A - \tilde{P}_i \overline{H})$ is Hurwitz. This property ensures that $\tilde{P}_i$ goes to a steady state asymptotically. This property is essential in understanding the global asymptotic behavior of $\hat{P}_i$ modeled by (19) with stochastic perturbations and coupled terms.

## IV. ILLUSTRATIVE EXAMPLE

As an illustration of the discussions presented in the previous section, let us consider an academic example with

$$ A = \begin{bmatrix} a & 2a & 3a & a \\ a & 0 & a & 4a \\ 0 & a & a & 2a \\ a & a & 0 & 2a \end{bmatrix}, \quad H = \begin{bmatrix} 1 & 0 & 1 & 0 \\ 0 & 1 & 0 & 1 \end{bmatrix}, $$

and we take $a = 0.1$. We choose $G = I_4$ (an identity matrix in $\mathbb{R}^{4 \times 4}$), and $V = 0.5 I_2$. The state process is described by (1). We run three different filtering algorithms: centralized Kalman filtering, ensemble Kalman filtering, and our proposed distributed filtering with empirical covariances.

The case of centralized Kalman filter is simulated with $dy = \text{col}(dy_1, \dots, dy_5)$, where $dy_i = H x dt + d\nu_i$ and $d\nu_i$ is an independent copy of the Wiener process $d\nu$ in $\mathbb{R}^2$ with infinitesimal variance $V$. The plots of the simulation are shown in Fig. 1a and Fig. 1b, where the evolution of the trace of error covariance matrix in Figure 1b is obtained from the standard Riccati equation.

For the ensemble Kalman filter (4), we take 5 particles and each of these particles uses a single observation process $dy =$
<figure><img src="image_3.png" /><figcaption>(a) Optimal state estimates with centralized Kalman filter.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>(c) Suboptimal estimates from ensemble Kalman filter.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(e) State estimates at Node 0 of (f) Empirical covariances at all the cyclic graph. (g) Nodes of the cyclic graph.</figcaption></figure>

<caption>Fig. 1: Plots for estimates (left column) and trace of error covariance matrix (right column) for different algorithms.</caption>

$Hxdt + dν$. The simulation results are plotted in Fig. 1c and Fig. 1d. In this case, Fig. 1d shows the evolution of empirical covariance which fluctuates around the optimal solution of centralized Kalman filter with varying intensity due to noise levels. The state estimates in Fig. 1c on the other hand are very close to the optimal ones.

Lastly, we consider the distributed filter where the underlying graph is undirected and cyclic with 5 nodes, so that for each node $i$, we have $d_i = 3$ and hence the regularity assumption holds. Each node has an independent observation process $dy_i = Hxdt + dν_i$, for $i = 1, \dots, 5$, with $dν_i$ being an independent copy of $dν$. The simulation results are plotted in Fig. 1e and Fig. 1f. Here, Fig. 1f shows the evolution of local empirical covariance at one node with fluctuations around a steady state, whereas the state estimates are very close to their optimal counterparts.

V. CONCLUSIONS AND PERSPECTIVES

We considered the problem of distributed filtering using a design technique that is inspired by the algorithms in ensemble filtering. The sensor nodes exchange information about their estimate, which allows them to compute the local covariance at each node. The innovation term used by the filter at each node is then multiplied by a gain term which depends on the empirical covariance computed from the estimates of the neighboring nodes. To build on this work,

one immediate question to address is the stability of process $(\hat{x}_i, \hat{P}_i)$ with complete dynamics.

REFERENCES

[1] F. Bullo, J. Cortés, and S. Martinez, *Distributed Control of Robotic Networks*. Princeton University Press, 2009.

[2] R. Kalman and R. Bucy, “New results in linear filtering and prediction theory,” *Journal of Basic Engineering*, vol. 83, no. 1, pp. 95–108, 1961.

[3] B. Anderson and J. Moore, *Optimal Filtering*. Dover Publications, 1979.

[4] J. Tsitsiklis, “Problems in decentralized decision making and computation,” Ph.D. dissertation, Massachusetts Institute of Technology, USA, 1984.

[5] V. Borkar and P. Varaiya, “Asymptotic agreement in distributed estimation,” *IEEE Transactions on Automatic Control*, vol. 27, no. 3, pp. 650–655, 1982.

[6] R. Olfati-Saber, “Distributed Kalman filtering for sensor networks,” in *Proc. of 46th IEEE Conf. on Decision and Control*, 2007.

[7] —, “Kalman-Consensus Filter : Optimality, stability, and performance,” in *Proc. of 48th IEEE Conf. on Decision and Control*, 2009.

[8] R. Carli, A. Chiuso, L. Schenato, and S. Zampieri, “Distributed Kalman filtering based on consensus strategies,” *IEEE Journal on Selected Areas in Communications*, vol. 26, no. 4, pp. 622–633, 2008.

[9] A. Speranzon, C. Fischione, K. Johansson, and A. Sangiovanni-Vincentelli, “A distributed minimum variance estimator for sensor networks,” *IEEE Journal on Selected Areas in Communications*, vol. 26, no. 4, pp. 609–621, 2008.

[10] I. Matei and J. Bara, “Consensus-based linear distributed filtering,” *Automatica*, vol. 48, no. 8, pp. 1776–1782, 2012.

[11] G. Battistelli and L. Chisci, “Kullback-Leibler average, consensus on probability densities, and distributed state estimation with guaranteed stability,” *Automatica*, vol. 50, pp. 707–718, 2014.

[12] A. Tanwani, “Suboptimal filtering over sensor networks with random communication,” *IEEE Transactions on Automatic Control*, vol. 67, no. 10, pp. 5456–5463, 2022.

[13] K. Ryu and J. Back, “Consensus optimization approach for distributed Kalman filtering: Performance recovery of centralized filtering,” *Automatica*, vol. 149, no. 110843, 2023.

[14] A. Nedic and J. Liu, “Distributed optimization for control,” *Annual Review of Control, Robotics, and Autonomous Systems*, vol. 1, pp. 77–103, 2018.

[15] A. Doucet, N. De Freitas, and N. Gordon, Eds., *Sequential Monte Carlo methods in practice*, 1st ed. Springer, 2001.

[16] G. Evensen, “Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte Carlo methods to forecast error statistics,” *Journal of Geophysical Research: Oceans*, vol. 99, pp. 10143–10162, 1994.

[17] Z. Chen, “Bayesian filtering: From Kalman filters to particle filters, and beyond,” *Statistics*, vol. 182, no. 1, pp. 1–69, 2003.

[18] E. Calvello, S. Reich, and A. Stuart, “Ensemble Kalman methods: A mean field perspective,” *Acta Numerica*, vol. 34, pp. 123–291, 2025.

[19] G. Evensen, *Data Assimilation: The Ensemble Kalman Filter*, 2nd ed. Springer Science & Business Media, 2009.

[20] S. Särkkä and L. Svensson, *Bayesian Filtering and Smoothing*, 2nd ed. Cambridge University Press, 2023.

[21] J. de Wiljes, S. Reich, and W. Stannat, “Long-time stability and accuracy of the ensemble Kalman–Bucy filter for fully observed processes and small measurement noise,” *SIAM Journal on Applied Dynamical Systems*, vol. 17, no. 2, pp. 1152–1181, 2018.

[22] P. Del Moral and J. Tugaut, “On the stability and the uniform propagation of chaos properties of ensemble Kalman–Bucy filters,” *The Annals of Applied Probability*, vol. 28, no. 2, pp. 790–850, 2018.

[23] A. N. Bishop and P. Del Moral, “On the mathematical theory of ensemble (linear-Gaussian) Kalman–Bucy filtering,” *Math. Control, Signals and Systems*, 2023.

[24] A. Taghvaei and P. Mehta, “A survey of feedback particle filter and related controlled interacting particle systems (CIPS),” *Annual Reviews in Control*, vol. 55, pp. 356–378, 2023.

[25] —, “An optimal transport formulation of the ensemble Kalman filter,” *IEEE Transactions on Automatic Control*, vol. 66, no. 7, pp. 3052–3067, 2021.

[26] B. Øksendal, *Stochastic differential equations*. Springer, 2003.