Stochastic Dynamics of Noisy Average Consensus: Analysis and Optimization

The following notation will be used throughout this paper. The symbols $\mathbb{R}$ and $\mathbb{R}_{+}$ represent the set of real numbers and the set of positive real numbers, respectively. The one dimensional Gaussian distribution with mean $\mu$ and variance $\sigma^{2}$ is denoted by ${\cal N}(\mu,\sigma^{2})$. The multivariate Gaussian distribution with mean vector $\bm{\mu}$ and covariance matrix $\bm{\Sigma}$ is represented by ${\cal N}(\bm{\mu},\bm{\Sigma})$. The expectation operator is denoted by ${\sf E}[\cdot]$. The notation $\mbox{diag}(\bm{x})$ is the diagonal matrix whose diagonal elements are given by $\bm{x}\in\mathbb{R}^{n}$. The matrix exponential $\exp(\bm{X})(\bm{X}\in\mathbb{R}^{n\times n})$ is defined by

$\displaystyle\exp(\bm{X})\equiv\sum_{k=0}^{\infty}\frac{1}{k!}\bm{X}^{k}.$

(3)

The Frobenius norm of $\bm{X}\in\mathbb{R}^{n\times n}$ is denoted by $\|\bm{X}\|_{F}$. The notation $[n]$ denotes the set of consecutive integers from $1$ to $n$.

Let $G\equiv(V,E)$ be a connected undirected graph where $V=[n]$. Suppose that a node $i\in V$ can be regarded as an agent communicating over the graph $G$. Namely, a node $i$ and a node $j$ can communicate with each other if $(i,j)\in E$. We will not distinguish $(i,j)$ and $(j,i)$ because the graph $G$ is undirected.

Each node $i$ has a state value $x_{i}(t)\in\mathbb{R}$ where $t\in\mathbb{R}$ represents continuous-time variable. The neighborhood of a node $i\in V$ is represented by

$\displaystyle{\cal N}_{i}\equiv\{j\in V:(j,i)\in E,i\neq j\}.$

(4)

Note that the node $i$ is excluded from ${\cal N}_{i}$. For any time $t$, a node $i\in V$ can access the self-state $x_{i}(t)$ and the state values of its neighborhood, i.e., $x_{j}(t),j\in{\cal N}_{i}$ but cannot access to the other state values.

In this section, we briefly review the basic properties of the average consensus processes [3]. We now assume that the set of state values $\bm{x}(t)\equiv(x_{1}(t),x_{2}(t),\ldots,x_{n}(t))^{T}$ are evolved according to the simultaneous differential equations

$\displaystyle\frac{dx_{i}(t)}{dt}=-\sum_{j\in{\cal N}_{i}}\mu_{ij}(x_{i}(t)-x_{j}(t)),\quad i\in[n],$

(5)

where the initial condition is

$\displaystyle\bm{x}(0)=\bm{c}\equiv(c_{1},c_{2},\ldots,c_{n})^{T}\in\mathbb{R}^{n}.$

(6)

The edge weight $\mu_{ij}$ follows the symmetric condition

$\displaystyle\mu_{ij}=\mu_{ji},\quad i\in[n],j\in[n].$

(7)

Let $\bm{\Delta}\equiv(\Delta_{1},\Delta_{2},\ldots,\Delta_{n})^{T}\in\mathbb{R}_{+}^{n}$ be a degree sequence where $\Delta_{i}$ is defined by

$\displaystyle\Delta_{i}\equiv\sum_{j\in{\cal N}_{i}}\mu_{ij},\quad i\in[n].$

(8)

The continuous-time dynamical system (5) is called an average consensus system because a state value converges to the average of the initial state values at the limit of $t\rightarrow\infty$, i.e,

$\displaystyle\lim_{t\rightarrow\infty}\bm{x}(t)=\frac{1}{n}\left(\sum_{i=1}^{n}c_{i}\right)\bm{1}{=\gamma\bm{1}},$

(9)

where the vector $\bm{1}$ represents $(1,1,\ldots,1)^{T}$ and $\gamma$ is defined by

$\displaystyle\gamma\equiv\frac{1}{n}\sum_{i=1}^{n}c_{i}.$

(10)

We define the Laplacian matrix $\bm{L}\equiv\{L_{ij}\}\in\mathbb{R}^{n\times n}$ of this consensus system as follows:

	$\displaystyle L_{ij}$	$\displaystyle=\Delta_{i},\quad i=j,\ i\in[n],$		(11)
	$\displaystyle L_{ij}$	$\displaystyle=-\mu_{ij},\quad i\neq j\mbox{ and }(i,j)\in E,$		(12)
	$\displaystyle L_{ij}$	$\displaystyle=0,\quad i\neq j\mbox{ and }(i,j)\notin E.$		(13)

From this definition, a Laplacian matrix satisfies

	$\displaystyle\bm{L}\bm{1}$	$\displaystyle=\bm{0},$		(14)
	$\displaystyle\mbox{diag}(\bm{L})$	$\displaystyle=\bm{\Delta},$		(15)
	$\displaystyle\bm{L}$	$\displaystyle=\bm{L}^{T}.$		(16)

Note that the eigenvalues of the Laplacian matrix $\bm{L}$ are nonnegative real because $\bm{L}$ is a positive semi-definite symmetric matrix. Let $\lambda_{1}=0<\lambda_{2}\leq\ldots\leq\lambda_{n}$ be the eigenvalues of $\bm{L}$ and $\bm{\xi}_{1},\bm{\xi}_{2},\ldots,\bm{\xi}_{n}$ be the corresponding orthonormal eigenvectors. The first eigenvector $\bm{\xi}_{1}\equiv(1/\sqrt{n})\bm{1}$ is corresponding to the eigenvalue $\lambda_{1}=0$, which results in $\bm{L}\bm{\xi}_{1}=0$.

By using the notion of the Laplacian matrix, the dynamical system (5) can be compactly rewritten as

$\displaystyle\frac{d\bm{x}(t)}{dt}=-\bm{L}\bm{x}(t),$

(17)

where the initial condition is $\bm{x}(0)=\bm{c}$. The dynamical behaviors of the average consensus system (17) are thus characterized by the Laplacian matrix $\bm{L}$. Since the ODE (17) is a linear ODE, it can be easily solved. The solution of the ODE (17) is given by

$\displaystyle\bm{x}(t)=\exp(-\bm{L}t)\bm{x}(0),\quad t\geq 0.$

(18)

Let $\bm{U}\equiv(\bm{\xi}_{1},\bm{\xi}_{2},\ldots,\bm{\xi}_{n})\in\mathbb{R}^{n\times n}$ where $\bm{U}$ is an orthogonal matrix. The Laplacian matrix $\bm{L}$ can be diagonalized by using $\bm{U}$, i.e.,

$\displaystyle\bm{L}=\bm{U}\mbox{diag}(\lambda_{1},\ldots,\lambda_{n})\bm{U}^{T}.$

(19)

On the basis of the diagonalization, we have the spectral expansion of the matrix exponential:

	$\displaystyle\exp(-\bm{L}t)$	$\displaystyle=\exp(-\bm{U}\mbox{diag}(\lambda_{1},\ldots,\lambda_{n})\bm{U}^{T}t)$
		$\displaystyle=\bm{U}\exp(-\mbox{diag}(\lambda_{1},\ldots,\lambda_{n})t)\bm{U}^{T}$
		$\displaystyle=\sum_{i=1}^{n}\exp(-\lambda_{i}t)\bm{\xi}_{i}\bm{\xi}_{i}^{T}.$		(20)

Substituting this to $\bm{x}(t)=\exp(-\bm{L}t)\bm{x}(0)$, we immediately have

$\displaystyle\bm{x}(t)=\frac{1}{n}\bm{1}(\bm{1}^{T})\bm{c}+\sum_{i=2}^{n}\exp(-\lambda_{i}t)\bm{\xi}_{i}\bm{\xi}_{i}^{T}\bm{c}.$

(21)

The second term of the right-hand side converges to zero since $\lambda_{k}>0$ for $k=2,3,\ldots,n$. This explains why average consensus happens, i.e., the convergence to the average of the initial state values (9). The second smallest eigenvalue $\lambda_{2}$, called algebraic connectivity [10], determines the convergence speed because $\exp(-\lambda_{2}t)\bm{\xi}_{2}\bm{\xi}_{2}^{\sf T}$ shows the slowest convergence in the second term.

The dynamical model (2) containing a white Gaussian noise process is mathematically challenging to handle. We will use a common approach of approximating the white Gaussian process by using the standard Wiener process. Instead of model (2), we will focus on the following stochastic differential equation (SDE)  [8]

$$d\bm{x}(t)=-\bm{L}\bm{x}(t)dt+\alpha d\bm{b}(t)$$

(22)

to study the noisy average consensus system. The parameter $\alpha$ is a positive real number, and it represents the intensity of the noises. The stochastic term $\bm{b}(t)$ represents the $n$-dimensional standard Wiener process. The elements of $\bm{b}(t)=(b_{1}(t),b_{2}(t),\ldots,b_{n}(t))^{T}$ are independent one dimensional-standard Wiener processes. For the Wiener process $b(t)$, we have $b(0)=0$, $E[b(t)]=0$, and it satisfies

$\displaystyle b(t)-b(s)\sim{\cal N}(0,t-s),\ 0\leq s\leq t.$

(23)

Our primary objective in the following analysis is to investigate the stochastic dynamics of the noisy average consensus system, focusing on deriving the mean and covariance of the solution $\bm{x}(t)$ for the SDE (22).

There are two approaches to analyze the system. The first approach relies on the established theory of Ito calculus [8], which is used to handle stochastic integrals directly (see Fig. 1). Ito calculus can be applied to derive the first and second moments of the solution of (22).

Alternatively, the second approach employs the Euler-Maruyama (EM) method [7] and utilizes the weak convergence property [7] of the EM method. We will adopt the latter approach in our analysis, as it does not require knowledge of advanced stochastic calculus if we accept the weak convergence property. Additionally, this approach can be naturally extended to the analysis on the discrete-time noisy average consensus system. Furthermore, the EM method plays a key role in the optimization method to be presented in Section V. Our analysis motivates the use EM method for optimizing the covariance.

We use the Euler-Maruyama method corresponding to this SDE so as to study the stochastic behavior of the solution of the SDE (22) defined above. The EM method is well-known numerical method for solving SDEs [7].

Assume that we need numerical solutions of a SDE in the time interval $0\leq t\leq T$. We divide this interval into $N$ bins and let $t_{k}\equiv k\eta,\ k=0,1,\ldots,N$ where the interval $\eta$ is given by $\eta\equiv{T}/{N}.$ Let us define a discretized sample $\bm{x}^{(k)}$ be $\bm{x}^{(k)}\equiv\bm{x}(t_{k}).$ It should be noted that, the choice of the width $\eta$ is crucial in order to ensure the stability and the accuracy of the EM method. A small width leads to a more accurate solution, but requires more computational time. A large width may be computationally efficient but may lead to instability in the solution.

The recursive equation of the EM method corresponding to SDE (22) is given by

$\displaystyle\bm{x}^{(k+1)}=\bm{x}^{(k)}-\eta\bm{L}\bm{x}^{(k)}+\alpha\bm{w}^{(k)},\ k=0,1,2,\ldots,N,$

(24)

where each element of $\bm{w}^{(k)}\equiv(w_{1}^{(k)},w_{2}^{(k)},\ldots,w_{n}^{(k)})^{T}$ follows $w_{i}^{(k)}\sim{\cal N}(0,\eta).$ In the following discussion, we will use the equivalent expression [7]:

$\displaystyle\bm{x}^{(k+1)}=\bm{x}^{(k)}-\eta\bm{L}\bm{x}^{(k)}+\alpha\sqrt{\eta}\bm{z}^{(k)},\ k=0,1,2,\ldots,N,$

(25)

where $\bm{z}^{(k)}$ is a random vector following the multivariate Gaussian distribution ${\cal N}(\bm{0},\bm{I})$. The initial vector $\bm{x}^{(0)}$ is set to be $\bm{c}$. This recursive equation will be referred to as the Euler-Maruyama recursive equation.

Figure 2 presents a solution evaluated with the EM method. The cycle graph with 10 nodes with the degree sequence $\bm{d}=(2,2,\ldots,2)$ is assumed. The initial value is randomly initialized as $\bm{x}(0)\sim{\cal N}(0,\bm{I})$. We can confirm that the state values are certainly converging to the average value $\gamma$ in the case of noiseless case (left). On the other hand, the state vector fluctuates around the average in the noisy case (right).

In the following, we will analyze the stochastic behavior of the residual error. This will be the basis for the MSE formula to be presented.

Recall that the initial state vector is $\bm{c}=(c_{1},c_{2},\ldots,c_{n})^{T}$ and that the average of the initial values is denoted by $\gamma$. Since the set of eigenvectors $\{\bm{\xi}_{1},\ldots,\bm{\xi}_{n}\}$ of $\bm{L}$ is an orthonormal base, we can expand the initial state vector $\bm{c}$ as

$\displaystyle\bm{c}=\zeta_{1}\bm{\xi}_{1}+\zeta_{2}\bm{\xi}_{2}+\cdots+\zeta_{n}\bm{\xi}_{n},$

(26)

where the coefficient is obtained by $\zeta_{i}=\bm{c}^{T}\bm{\xi}_{i}(i\in[n])$. Note that $\zeta_{1}\bm{\xi}_{1}=\gamma\bm{1}$ holds.

At the initial index $k=0$, the Euler-Maruyama recursive equation becomes

$\displaystyle\bm{x}^{(1)}=\bm{x}^{(0)}-\eta\bm{L}\bm{x}^{(0)}+\alpha\sqrt{\eta}\bm{z}^{(0)}.$

(27)

Substituting (26) into the above equation, we have

	$\displaystyle\bm{x}^{(1)}$	$\displaystyle=\bm{x}^{(0)}-\eta\bm{L}(\zeta_{1}\bm{\xi}_{1}+\zeta_{2}\bm{\xi}_{2}+\cdots+\zeta_{n}\bm{\xi}_{n})+\alpha\sqrt{\eta}\bm{z}^{(0)}$
		$\displaystyle=\bm{x}^{(0)}-\eta\zeta_{1}\bm{L}\bm{\xi}_{1}-\eta L(\bm{x}^{(0)}-\zeta_{1}\bm{\xi}_{1})+\alpha\sqrt{\eta}\bm{z}^{(0)}$
		$\displaystyle={\bm{x}^{(0)}-\eta\bm{L}(\bm{x}^{(0)}-\gamma\bm{1})+\alpha\sqrt{\eta}\bm{z}^{(0)}},$		(28)

where the equations $L\bm{\xi}_{1}=\bm{0}$ and $\zeta_{1}\bm{\xi}_{1}=\gamma\bm{1}$ are used in the last equality. Subtracting $\gamma\bm{1}$ from the both sides, we get

$\displaystyle\bm{x}^{(1)}-\gamma\bm{1}=(\bm{I}-\eta\bm{L})(\bm{x}^{(0)}-\gamma\bm{1})+\alpha\sqrt{\eta}\bm{z}^{(0)}.$

(29)

For the index $k\geq 1$, the Euler-Maruyama recursive equation can be written as

$\displaystyle\bm{x}^{(k+1)}=(\bm{I}-\eta\bm{L})\bm{x}^{(k)}+\alpha\sqrt{\eta}\bm{z}^{(k)}.$

(30)

Subtracting $\gamma\bm{1}$ from the both sides, we have

$\displaystyle\bm{x}^{(k+1)}-\gamma\bm{1}=(\bm{I}-\eta\bm{L})\bm{x}^{(k)}-\gamma\bm{1}+\alpha\sqrt{\eta}\bm{z}^{(k)}.$

(31)

By using the relation $(\bm{I}-\eta\bm{L})\gamma\bm{1}=\gamma\bm{1},$ we can rewrite the above equation as

	$\displaystyle\bm{x}^{(k+1)}-\gamma\bm{1}$	$\displaystyle=(\bm{I}-\eta\bm{L})\bm{x}^{(k)}-(\bm{I}-\eta\bm{L})\gamma\bm{1}+\alpha\sqrt{\eta}\bm{z}^{(k)}$
		$\displaystyle=(\bm{I}-\eta\bm{L})(\bm{x}^{(k)}-\gamma\bm{1})+\alpha\sqrt{\eta}\bm{z}^{(k)}.$		(32)

It can be confirmed the above recursion (32) is consistent with the initial equation (29). We here summarize the above argument as the following lemma.

The residual error $\bm{e}^{(k)}$ denotes the error between the average vector $\gamma\bm{1}$ and the state vector $\bm{x}^{(k)}$ at time index $k$. By analyzing the statistical behavior of $\bm{e}^{(k)}$, we can gain insight into the stochastic properties of the dynamics of the noisy consensus system.

Let a vector $\bm{x}\sim{\cal N}(\bm{\mu},\bm{\Sigma})$. Recall that the vector obtained by a linear map $\bm{y}=\bm{A}\bm{x}$ also follows the Gaussian distribution, i.e.,

$\displaystyle\bm{y}\sim{\cal N}(\bm{A}\bm{\mu},\bm{A}\bm{\Sigma}\bm{A}^{T}),$

(34)

where $\bm{A}\in\mathbb{R}^{n\times n}$. If two Gaussian vectors $\bm{a}\sim{\cal N}(\bm{\mu}_{a},\bm{\Sigma}_{a})$ and $\bm{b}\sim{\cal N}(\bm{\mu}_{b},\bm{\Sigma}_{b})$ are independent, the sum $\bm{z}=\bm{a}+\bm{b}$ becomes also Gaussian, i.e,

$\displaystyle\bm{z}\sim{\cal N}(\bm{\mu}_{a}+\bm{\mu}_{b},\bm{\Sigma}_{a}+\bm{\Sigma}_{b}).$

(35)

In the recursive equation (33), it is evident that $\bm{e}^{(1)}$ follows a multivariate Gaussian distribution because

$\displaystyle\bm{e}^{(1)}=(\bm{I}-\eta\bm{L})(\bm{c}-\gamma\bm{1})+\alpha\sqrt{\eta}\bm{z}^{(0)}$

(36)

is the sum of a constant vector and a Gaussian random vector. From the above properties of Gaussian random vectors, the residual error vector $\bm{e}^{(k)}$ follows the multivariate Gaussian distribution ${\cal N}(\bm{\mu}^{(k)},\bm{\Sigma}^{(k)})$ where the mean vector $\bm{\mu}^{(k)}$ and the covariance matrix $\bm{\Sigma}^{(k)}$ are recursively determined by

	$\displaystyle\bm{\mu}^{(k+1)}$	$\displaystyle=(\bm{I}-\eta\bm{L})\bm{\mu}^{(k)},$		(37)
	$\displaystyle\bm{\Sigma}^{(k+1)}$	$\displaystyle=(\bm{I}-\eta\bm{L})\bm{\Sigma}^{(k)}(\bm{I}-\eta\bm{L})^{T}+\alpha^{2}\eta\bm{I}$		(38)

for $k\geq 0$ where the initial values are formally given by

	$\displaystyle\bm{\mu}^{(0)}$	$\displaystyle=\bm{c}-\gamma\bm{1},$		(39)
	$\displaystyle\bm{\Sigma}^{(0)}$	$\displaystyle=\bm{O}.$		(40)

Solving the recursive equation, we can get the asymptotic mean formula as follows.

(Proof) The mean recursion is given as $\bm{\mu}^{(k)}=(\bm{I}-\eta\bm{L})^{k}(\bm{c}-\gamma\bm{1})$ for $k\geq 1$. Recall that the eigenvalue decomposition of $\bm{L}$ is given by $\bm{L}=\bm{U}\mbox{diag}(\lambda_{1},\ldots,\lambda_{n})\bm{U}^{T}$. From

$\displaystyle\bm{I}-\eta\bm{L}=\bm{U}(\bm{I}-\eta\mbox{diag}(\lambda_{1},\ldots,\lambda_{n}))\bm{U}^{T},$

(42)

we have

$\displaystyle(\bm{I}-\eta\bm{L})^{k}=\bm{U}\mbox{diag}((1-\eta\lambda_{1})^{k},\ldots,(1-\eta\lambda_{n})^{k})\bm{U}^{T}.$

(43)

This implies, from the definition of exponential function,

$\displaystyle\lim_{N\rightarrow\infty}\left(\bm{I}-\frac{T}{N}\bm{L}\right)^{N}=\exp(-\bm{L}T),$

(44)

where $\eta=T/N$.

It is easy to confirm that the claim of this lemma is consistent with the continuous solution of noiseless case (18). Namely, at the limit of $\alpha\rightarrow 0$, the state evolution of the noisy system converges to that of the noiseless system.

We here discuss the asymptotic behavior of the covariance matrix $\bm{\Sigma}^{(N)}$ at the limit of $N\rightarrow\infty$.

(Proof) Recall that

$\displaystyle\bm{I}-\eta\bm{L}=\bm{U}\mbox{diag}(1,1-\eta\lambda_{2}\ldots,1-\eta\lambda_{n})\bm{U}^{T}.$

(47)

Let $\bm{\Sigma}^{(k)}=\bm{U}\mbox{diag}(s_{1}^{(k)},\ldots,s_{n}^{(k)})\bm{U}^{T}$. A spectral representation of the covariance evolution (38) is thus given by

	$\displaystyle\mbox{diag}(s_{1}^{(k+1)},\ldots,s_{n}^{(k+1)})$
	$\displaystyle=\mbox{diag}(s_{1}^{(k)},s_{2}^{(k)}(1-\eta\lambda_{2})^{2}\ldots,s_{n}^{(k)}(1-\eta\lambda_{n})^{2})+\alpha^{2}\eta\bm{I},$		(48)

where $s_{i}^{(0)}=0$. The first component follows a recursion $s_{1}^{(k+1)}=s_{1}^{(k)}+\alpha^{2}\eta$ and thus we have $s_{1}^{(N)}=\alpha^{2}\eta N=\alpha^{2}T.$ Another component follows

$\displaystyle s_{i}^{(k+1)}=s_{i}^{(k)}(1-\eta\lambda_{i})^{2}+\alpha^{2}\eta.$

(49)

Let us consider the characteristic equation of (49) which is given by

$\displaystyle s=s(1-\eta\lambda_{i})^{2}+\alpha^{2}\eta.$

(50)

The solution of the equation is given by

$\displaystyle s=\frac{\alpha^{2}\eta}{1-(1-\eta\lambda_{i})^{2}}.$

(51)

The above recursive equation (49) thus can be transformed as

$\displaystyle s_{i}^{(k+1)}-s=(s_{i}^{(k)}-s)(1-\eta\lambda_{i})^{2}.$

(52)

From the above equation, $s_{i}^{(N)}$ can be solved as

$\displaystyle s_{i}^{(N)}=s+(s_{i}^{(0)}-s)(1-\eta\lambda_{i})^{2N}.$

(53)

Taking the limit $N\rightarrow\infty$, we have

$\displaystyle\lim_{N\rightarrow\infty}s_{i}^{(N)}=\frac{\alpha^{2}}{2\lambda_{i}}\left(1-e^{-2\lambda_{i}T}\right).$

(54)

We thus have the claim of this lemma.

As previously noted, the asymptotic mean (41) is consistent with the continuous solution. The weak convergence property of the EM method [7] allows us to obtain the moments of the error at time $t$.

We will briefly explain the weak convergence property. Suppose a SDE with the form:

$\displaystyle d\bm{x}(t)=\phi(\bm{x}(t))dt+\psi(\bm{x}(t))d\bm{b}(t).$

(55)

If $\phi$ and $\psi$ are bounded and Lipschitz continuous, then the finite order moment estimated by the EM method converges to the exact moment of the solution $\bm{x}(t)$ at the limit $N\rightarrow\infty$ [7]. This property is called the weak convergence property. In our case, the SDE (22) has bounded and Lipschitz continuous coefficient functions, i.e, $\phi(\bm{x})=-\bm{L}\bm{x}$ and $\psi(\bm{x})=\alpha$. Hence, we can employ the weak convergence property in our analysis.

Suppose $\bm{x}(t)$ is a solution of SDE (22) with the initial condition $\bm{x}(0)=\bm{c}$. Let $\bm{\mu}(t)$ be the mean vector of the residual error $\bm{e}(t)=\bm{x}(t)-\gamma\bm{1}$ and $\bm{\Sigma}(t)$ is the covariance matrix of the residual error $\bm{e}(t)$.

(Proof) Due to the weak convergence property of the EM method, the first and second moments of the error are converged to the asymptotic mean and covariance of the EM method [7], i.e.,

	$\displaystyle\bm{\mu}(T)$	$\displaystyle=\lim_{N\rightarrow\infty}\bm{\mu}^{(N)}$		(58)
	$\displaystyle\bm{\Sigma}(T)$	$\displaystyle=\lim_{N\rightarrow\infty}\bm{\Sigma}^{(N)},$		(59)

where $N$ and $T$ are related by $T=\eta N$. Applying Lemmas 2 and 3 and replacing the variable $T$ by $t$ provide the claim of the theorem.

In the following, we assume that the initial state vector $\bm{c}$ follows Gaussian distribution ${\cal N}(\bm{0},\bm{I})$.

In this setting, $\bm{\mu}(t)$ also follows multivariate Gaussian distribution with the mean vector $\bm{0}$ and the covariance matrix $\bm{Q}(t)\bm{Q}(t)^{T}$ where

$\displaystyle\bm{Q}(t)\equiv\exp(-\bm{L}t)\left(\bm{I}-\frac{1}{n}\bm{1}(\bm{1}^{T})\right)$

(60)

because $\bm{\mu}(t)$ can be rewritten as

$\displaystyle\bm{\mu}(t)$

$\displaystyle=\exp(-\bm{L}t)(\bm{c}-\gamma\bm{1})=\exp(-\bm{L}t)\left(\bm{I}-\frac{1}{n}\bm{1}(\bm{1}^{T})\right)\bm{c}.$

(61)

By using the result of Theorem 1, we immediately have the following corollary indicating the MSE formula.

(Proof) We can rewrite $\bm{x}(t)$ as:

$\displaystyle\bm{x}(t)=\gamma\bm{1}+\bm{Q}(t)\bm{c}+\bm{w},$

(63)

where $\bm{w}\sim{\cal N}(\bm{0},\bm{\Sigma}(t))$, and $\bm{w}$ and $\bm{c}$ are independent. We thus have

	$\displaystyle{\sf MSE}(t)$	$\displaystyle=\mbox{tr}(\bm{\Sigma}(t))+\mbox{tr}(\bm{Q}(t)\bm{Q}(t)^{T})$
		$\displaystyle=\alpha^{2}t+\frac{\alpha^{2}}{2}\sum_{i=2}^{n}\frac{1-e^{-2\lambda_{i}t}}{\lambda_{i}}+\mbox{tr}(\bm{Q}(t)\bm{Q}(t)^{T})$		(64)

due to Theorem 1.

Since the value of the term $\mbox{tr}(\bm{Q}(t)\bm{Q}(t)^{T})$ is exponentially decreasing with $t$, $\mbox{tr}(\bm{\Sigma}(t))$ is dominant in ${\sf MSE}(t)$ for sufficiently large $t$. For sufficiently large $t$, the MSE is well approximated by the asymptotic MSE (AMSE) as

$\displaystyle{\sf MSE}(t)\simeq{\sf AMSE}(t)\equiv\alpha^{2}t+\frac{\alpha^{2}}{2}\sum_{i=2}^{n}\frac{1}{\lambda_{i}}$

(65)

because $\mbox{tr}(\bm{Q}(t)\bm{Q}(t)^{T})$ is negligible, and $1-e^{-2\lambda_{i}t}$ can be well approximated to $1$. We can observe that the sum of inverse eigenvalue $\sum_{i=2}^{n}({1}/{\lambda_{i}})$ of the Laplacian matrix determines the intercept of the ${\sf AMSE}(t)$. In other words, the graph topology influences the stochastic error behavior through the sum of inverse eigenvalues of the Laplacian matrix.

Figure 3 presents a comparison of ${\sf MSE}(t)$ evaluated by the EM method (25) and the formula in (64). In this experiment, the cycle graph with 10 nodes is used. The values of ${\sf AMSE}(t)$ are also included in Fig. 3. We can see that the theoretical values of ${\sf MSE}(t)$ and estimated values by the EM method are quite close.

In the previous section, we demonstrated that the MSE can be expressed in closed-form. It is natural to optimize the edge weights $\{\mu_{ij}\}$ in order to decrease the value of the MSE. The optimization of the edge weights is equivalent to the optimization of the Laplacian matrix $\bm{L}$. There exist several related works that aim to achieve a similar goal for noise-free systems. For example, Xiao and Boyd [5] proposed a method to minimize the second eigenvalue to achieve the fastest convergence to the average. They formulated the optimization problem as a convex optimization problem, which can be solved efficiently. Kishida et al. [13] presented a deep unfolding-based method for optimizing time-dependent edge weights, yet these methods are not applicable to systems with noise. Optimizing the MSE may be a non-trivial task as it involves the sum of the inverse eigenvalues of the Laplacian matrix.

In this subsection, we will present two optimization problems of edge weights.

Advances in deep neural networks have had a strong impact on the design of algorithms for communications and signal processing [14, 15, 16]. Deep unfolding can be seen as a very effective way to improve the convergence of iterative algorithms. Gregor and LeCun introduced the Learned ISTA (LISTA) [21]. Borgerding et al. also proposed variants of AMP and VAMP with trainable capability [19][20]. Trainable ISTA(TISTA) [23] is another trainable sparse signal recovery algorithm with fast convergence. TISTA requires only a small number of trainable parameters, which provides a fast and stable training process. Another advantage of deep unfolding is that it has a relatively high interpretability of learning results.

The concept behind deep unfolding is rather simple. We can embed trainable parameters into the original iterative algorithm and then unfold the signal-flow graph of the original algorithm. The standard supervised training techniques used in deep learning, such as Stochastic Gradient Descent (SGD) and back propagation, can then be applied to the unfolded signal-flow graph to optimize the trainable parameters.

The combination of deep unfolding and the differential equation solvers [24] is a current area of active research in scientific machine learning. It should be noted, however, that the technique is not limited to applications within scientific machine learning. In this subsection, we introduce an optimization algorithm that is based on the deep-unfolded EM method. The central idea is to use a loss function that approximates ${\sf MSE}(t^{*})$. By using a stochastic gradient descent approach with this loss function, we can obtain a near-optimal solution for both optimization problems A and B. The proposed method can be easily implemented using any modern neural network framework that includes an automatic differentiation mechanism. The following subsections will provide a more detailed explanation of the proposed method.