Spatial localization for nonlinear dynamical stochastic models for excitable media

TL;DR

This paper develops a mathematical framework to analyze spatial localization in high-dimensional nonlinear stochastic models of excitable media, with applications in data assimilation and covariance sampling.

Contribution

It introduces a rigorous theory for covariance decay in coupled nonlinear stochastic systems, validated through models like FitzHugh-Nagumo.

Findings

Covariance decay behavior is characterized by local and nonlocal effects.

Theoretical results enable efficient covariance sampling strategies.

Validation with FitzHugh-Nagumo model confirms the framework's applicability.

Abstract

Nonlinear dynamical stochastic models are ubiquitous in different areas. Excitable media models are typical examples with large state dimensions. Their statistical properties are often of great interest but are also very challenging to compute. In this article, a theoretical framework to understand the spatial localization for a large class of stochastically coupled nonlinear systems in high dimensions is developed. Rigorous mathematical theories show the covariance decay behavior due to both local and nonlocal effects, which result from the diffusion and the mean field interaction, respectively. The analysis is based on a comparison with an appropriate linear surrogate model, of which the covariance propagation can be computed explicitly. Two important applications of these theoretical results are discussed. They are the spatial averaging strategy for efficiently sampling the covariance matrix and the localization technique in data assimilation. Test examples of a surrogate linear model and a stochastically coupled FitzHugh-Nagumo model for excitable media are adopted to validate the theoretical results. The latter is also used for a systematical study of the spatial averaging strategy in efficiently sampling the covariance matrix in different dynamical regimes.

Figures

Equations

$$\frac{d}{dt}{u}_t=\nu \Delta u_t+F(u_t)+\int h(u_t(x))dx+\dot{W}(t,x),$$

$$d{u}_i(t)=\frac{\nu }{2}{N}^2(u_{i+1}(t)+u_{i-1}(t)-2u_i(t))dt+F(u_i(t))dt+\frac{1}{N}\sum _{j=1}^{N}h(u_j(t))+d{w}_i(t),$$

$$d{\mathbf{x}}_i=\mathbf{f}(t,{\mathbf{x}}_{i-1},{\mathbf{x}}_i,{\mathbf{x}}_{i+1})dt+\frac{1}{N}\sum _{j=1}^N\mathbf{h}(t,{\mathbf{x}}_j)dt+\Sigma d{\mathbf{w}}_i(t),i=1,\dots ,N.$$

$$\mathbf{d}(i,j)=\min \{|i-j|,|i+N-j|,|j-i+N|\}.$$

$$\begin{array}{c}{\lambda }_0:=\underset{i,\mathbf{x},s\le t}{\sup }\{{\lambda }_{max}\left({\nabla }_{{\mathbf{x}}_i}\mathbf{f}(s,{\mathbf{x}}_{i-1},{\mathbf{x}}_i,{\mathbf{x}}_{i+1})+{\nabla }_{{\mathbf{x}}_i}\mathbf{f}(s,{\mathbf{x}}_{i-1},{\mathbf{x}}_i,{\mathbf{x}}_{i+1}{)}^T\right)\},\\ {\lambda }_F:=\underset{i,\mathbf{x},s\le t}{\sup }\{\Vert {\nabla }_{{\mathbf{x}}_{i-1}}\mathbf{f}(s,{\mathbf{x}}_{i-1},{\mathbf{x}}_i,{\mathbf{x}}_{i+1})\Vert ,\Vert {\nabla }_{{\mathbf{x}}_{i+1}}\mathbf{f}(s,{\mathbf{x}}_{i+1},{\mathbf{x}}_i,{\mathbf{x}}_{i+1})\Vert \},\\ {\lambda }_H:=\underset{i,{\mathbf{x}}_i,s\le t}{\sup }\Vert {\nabla }_{{\mathbf{x}}_i}\mathbf{h}(s,{\mathbf{x}}_i)\Vert .\end{array}$$

$$\text{cov}(g({\mathbf{x}}_i(t)),g({\mathbf{x}}_j(t)))\le C_{\beta ,t}\Vert \nabla g{\Vert }_{\infty }^2\left(e^{-\beta \mathbf{d}(i,j)}+\frac{1}{N}\right).$$

$${\hat{g}}_{K,N}=\frac{1}{K}\sum _{j=1}^K\frac{1}{N}\sum _{i=1}^{N}g({\mathbf{x}}_i^{(j)}(t)).$$

$${\hat{g}}_N=\frac{1}{N}\sum _{i=1}^{N}g({\mathbf{x}}_i(t)).$$

$$\mathbb{E}{\hat{g}}_{K,N}=\mathbb{E}{\hat{g}}_N=\frac{1}{N}\sum _{i=1}^N\mathbb{E}g({\mathbf{x}}_i(t))=\mathbb{E}g({\mathbf{x}}_1(t)).$$

$$\text{var}{\hat{g}}_N=\frac{1}{N^2}\sum _{i,j=1}^N\text{cov}(g({\mathbf{x}}_i),g({\mathbf{x}}_j))=\frac{1}{N}\sum _{j=1}^N\text{cov}(g({\mathbf{x}}_1),g({\mathbf{x}}_j)).$$

$$\text{var}{\hat{g}}_N\le \frac{1}{N}\sum _{j=1}^N{C}_{\beta ,t}\Vert \nabla g{\Vert }_{\infty }^2\left(e^{-\beta \mathbf{d}(1,j)}+\frac{1}{N}\right)\le \frac{1}{N}\frac{C_{\beta ,t}\Vert \nabla g{\Vert }_{\infty }^2}{1-e^{-\beta }}.$$

$$\text{var}{\hat{g}}_{K,N}\le \frac{1}{NK}\frac{C_{\beta ,t}\Vert \nabla g{\Vert }_{\infty }^2}{1-e^{-\beta }}.$$

$$\hat{C}:=\frac{1}{K-1}\sum _{j=1}^K({\mathbf{x}}^{(j)}(t)-\bar{\mathbf{x}})\cdot ({\mathbf{x}}^{(j)}(t)-\bar{\mathbf{x}}(t){)}^T,\bar{\mathbf{x}}(t)=\frac{1}{K}\sum _{j=1}^K{\mathbf{x}}_t^{(j)}.$$

$$[C^L{]}_{i,j}=[C{]}_{i,j}{1}_{\mathbf{d}(i,j)\le L}.$$

$$\mathbb{P}(\Vert {\hat{C}}^L-C^L\Vert >ϵ)\le 8\exp \left(2\log N-cK\min \left\{\frac{ϵ}{2L\Vert C\Vert },\frac{{ϵ}^2}{4L^2\Vert C{\Vert }^2}\right\}\right).$$

$$\Vert C-C^L\Vert \le \underset{i}{\max }\sum _{\mathbf{d}(i,j)>L}\left|[C{]}_{i,j}\right|\le \frac{2C_{\beta ,t}{e}^{-\beta L}}{1-e^{-\beta }}.$$

$$[A_{i,j}{]}_{m,n}=[\{A{\}}_{i,j}{]}_{m,n}=[A{]}_{(i-1)q+m,(j-1)q+n},i,j=1,\cdots ,N,m,n=1,\cdots ,q.$$

$$d{\mathbf{x}}_i(t)={\mathbf{f}}_i(t,{\mathbf{x}}_{i-1},{\mathbf{x}}_i,{\mathbf{x}}_{i+1})dt+\frac{1}{N}\sum _{j=1}^N{\mathbf{h}}_i(t,{\mathbf{x}}_j)dt+\Sigma d{\mathbf{w}}_i(t),i=1,\dots ,N.$$

$$Q(s)=\exp \left({\int }_{t-s}^t{G}_{\ast }(r)dr\right)\exp {\left({\int }_{t-s}^t{G}_{\ast }(r)dr\right)}^T.$$

$$|\text{cov}(g({\mathbf{x}}_i(t)),g({\mathbf{x}}_j(t)))|\le \sqrt{q}\Vert \nabla g{\Vert }_{\infty }^2\left(\Vert {\Sigma }_0^2{\Vert }_F[Q(t){]}_{i,j}+\Vert {\Sigma }^2{\Vert }_F{\int }_0^t[Q(s){]}_{i,j}ds\right).$$

$$d{\mathbf{x}}_i=[G_{\ast }{]}_{i,i}{\mathbf{x}}_{i}dt+[G_{\ast }{]}_{i,i-1}{\mathbf{x}}_{i-1}dt+[G_{\ast }{]}_{i,i+1}{\mathbf{x}}_{i+1}dt+\Sigma d{\mathbf{w}}_i(t),i=1,\dots ,N.$$

$$d\mathbf{x}(t)=D_{G_{\ast }}\mathbf{x}dt+D_{\Sigma }d\mathbf{w}(t),\mathbf{x}(0)\sim \mathcal{N}(0,D_{{\Sigma }_0}).$$

$$\{D_{G_{\ast }}{\}}_{i,j}=[G_{\ast }{]}_{i,j}{\mathbf{I}}_q,\{D_{\Sigma }{\}}_{i,j}=\text{mathds}{1}_{i=j}\Sigma ,\{D_{{\Sigma }_0}{\}}_{i,j}=\text{mathds}{1}_{i=j}{\Sigma }_0.$$

$$\mathbf{x}(t)=q(t)\mathbf{x}(0)+{\int }_0^{t}q(s)D_{\Sigma }d\mathbf{w}(s),q(s):=\exp \left({\int }_{t-s}^t{D}_{G_{\ast }}(r)dr\right).$$

$$\text{cov}\mathbf{x}(t)=q(t)D_{{\Sigma }_0}{D}_{{\Sigma }_0}^{T}q(t{)}^T+{\int }_0^{t}q(s)D_{\Sigma }{D}_{\Sigma }^{T}q(s{)}^{T}ds.$$

$$\text{cov}(f(X),g(X))=\mathbb{E}f(X)g(X)-f(X)g(Y)={\int }_0^{\frac{\pi }{2}}\sin \theta \mathbb{E}⟨\nabla f(X),\nabla g(X^{\theta })⟩d\theta .$$

$$d{\mathbf{x}}_i(t)={\mathbf{f}}_i(t,{\mathbf{x}}_{i-1},{\mathbf{x}}_i,{\mathbf{x}}_{i+1})dt+\frac{1}{N}\sum _{j=1}^N{\mathbf{h}}_i(t,{\mathbf{x}}_j)dt,i=1,\dots ,N.$$

$$\mathbf{X}(0)=m_0+{\Sigma }_0{\mathbf{W}}_0,\mathbf{X}(n+1)={\Psi }_h(\mathbf{X}(n))+{\Sigma }_h{\mathbf{W}}_n,{\Sigma }_h=\sqrt{h}{D}_{\Sigma }.$$

$$\begin{array}{c}{\lambda }_0:=\underset{i,\mathbf{x},s\le t}{\sup }\{{\lambda }_{max}\left({\nabla }_{{\mathbf{x}}_i}{\mathbf{f}}_i(s,{\mathbf{x}}_{i-1},{\mathbf{x}}_i,{\mathbf{x}}_{i+1})+{\nabla }_{{\mathbf{x}}_i}{\mathbf{f}}_i(s,{\mathbf{x}}_{i-1},{\mathbf{x}}_i,{\mathbf{x}}_{i+1}{)}^T\right)\},\\ {\lambda }_F:=\underset{i,\mathbf{x},s\le t}{\sup }\{\Vert {\nabla }_{{\mathbf{x}}_{i-1}}{\mathbf{f}}_i(s,{\mathbf{x}}_{i-1},{\mathbf{x}}_i,{\mathbf{x}}_{i+1})\Vert ,\Vert {\nabla }_{{\mathbf{x}}_{i+1}}{\mathbf{f}}_i(s,{\mathbf{x}}_{i-1},{\mathbf{x}}_i,{\mathbf{x}}_{i+1})\Vert \},\\ {\lambda }_H:=\underset{i,j,\mathbf{x},s\le t}{\sup }\Vert {\nabla }_{{\mathbf{x}}_j}{\mathbf{h}}_i(s,{\mathbf{x}}_j)\Vert .\end{array}$$

$$[F_{\ast }{]}_{i,i}={\lambda }_0,[F_{\ast }{]}_{i,i\pm 1}={\lambda }_F,[H_{\ast }{]}_{i,j}=\frac{1}{N}{\lambda }_H.$$

Full text

(1)Department of Mathematics, University of Wisconsin - Madison, WI, USA

(2)Department of Mathematics and Center for Atmosphere Ocean Science, Courant Institute of Mathematical Sciences, New York University, New York, NY, USA

Center for Prototype Climate Modeling, New York University Abu Dhabi, Saadiyat Island, Abu Dhabi, UAE

(3)Department of Mathematics, National University of Singapore, Singapore 119076

1 Introduction

Nonlinear dynamical stochastic models are ubiquitous in different areas [35, 36, 52, 50]. Key features of these complex systems are multiscale dynamics, high-dimensional phase space, nonlinear energy transfers, highly non-Gaussian probability density functions (PDFs), intermittent instability, random internal and external forcing as well as extreme events. Understanding and predicting these complex turbulent dynamical systems have both scientific and practical significance. Many excitable media models belong to such types of nonlinear dynamical stochastic systems [33, 56, 45, 23]. These excitable media models are ubiquitous in nature and can be found in physical, chemical, and biological systems that are far from thermodynamic equilibrium. One important feature of these excitable medium is that they are susceptible to finite perturbations, where a variety of wave patterns can be triggered. Examples of excitable medium include solitary waves, target like patterns, and spiral waves.

For nonlinear dynamical stochastic models, the high dimensionality imposes great challenges to the understanding of the underlying mechanisms, the design of effective prediction schemes and various uncertainty quantification problems [4, 47]. In particular, the high dimensionality leads to both a large computational cost and an enormous storage even for one single model simulation. Moreover, it is often of interest to calculate the time evolution of the statistics, especially the covariance matrix that includes the information of the interactions between different components [1, 29]. Since analytical solutions for complex nonlinear systems are typically unavailable, Monte Carlo methods are often applied. However, the affordable sample size is very limited due to the computational cost while the size of the covariance matrix, which is the square of the state variables, can be many magnitudes larger in many real applications that involve millions or billions of state variables. Furthermore, non-Gaussian features due to the nonlinear nature of the underlying dynamics such as extreme events can greatly affect the calculation of the covariance matrix and introduce extra difficulties. Therefore, understanding the covariance structures in nonlinear dynamical stochastic models and then developing efficient computing strategies that can exploit these structures become extremely important.

To facilitate the calculation of the covariance matrix, various approximations are applied in practice. One widely used strategy is to adopt Gaussian approximations for the underlying dynamics since there are many well-developed analytical and numerical methods for sampling high-dimensional Gaussian distributions. However, completely ignoring the non-Gaussian information in the underlying dynamics often leads to large model errors. One effective strategy is to decompose the entire phase space into a high dimensional Gaussian part and a low dimensional non-Gaussian part using the conditional Gaussian framework [11, 12, 9, 10, 38, 32]. Another strategy is to exploit the fact that the majority of the problem uncertainty is distributed within a low dimensional subspace of the Karhuen-Loève expansion, then uncertainty quantification (UQ) can be operated just in this subspace for sampling and computational simplicity [51, 15, 13, 7, 19, 46, 14].

This paper focuses on a different structure-exploiting strategy known as localization, which is widely used in many areas including the numerical weather prediction [26, 22, 25] and efficiently quantifying the uncertainty in excitable media [6, 44]. In many applications, the system state variables represent the information at different spatial locations, and the dependence of certain information between different grid points decays quickly when their distance becomes large. In other words, two grid points can be regarded as nearly independent if they are far from each other. The localization strategy exploits this feature and adopts a band matrix as the approximation of the covariance [48, 2]. Thus, the localization allows the development of much simpler approaches in analyzing the covariance matrix and thus facilitates the study of many important UQ problems.

While the localization strategy has a wide range of real-world applications, it remains unclear about the role of different dynamical components on affecting the covariance localization. In addition, despite many empirical evidences, there is no rigorous justifiction of the covariance decay as a function of the spatial distance. Nevertheless, these are crucial topics in practice for at least two reasons. First, it is known that applying localization by approximating the full covariance matrix with a band matrix may introduce some biases in quantifying the uncertainty of the underlying system [53, 42] and the current bias estimates [53, 41] require the knowledge of the underlying covariance matrix as well as a quantitative description of the decay in the covariance. Second, in order to apply localization in practice, one often needs to specify a localization radius $L$, which represents the typical distance for two components to be roughly independent. Unfortunately, in the absence of a rigorous theory to understand the covariance decay behavior, determining $L$ often relies on exhaustive tuning and data fitting, which can be expensive and inaccurate.

This paper aims at developing a theoretical framework to understand the covariance behavior for a large class of stochastically coupled nonlinear systems in high dimensions that include many excitable media models and other nonlinear complex models. Rigorous mathematical analysis shows that the local effect from the diffusion results in an exponential decay of the components in the covariance matrix as a function of the distance while the global effect due to the mean field interaction synchronizes different components and contributes to a global covariance. To achieve this, a covariance comparison principle is established. Then by explicitly computing the covariance propagation in a linear surrogate model, an upper bound on the covariance entries can be obtained.

Two important applications of these theoretical results are illustrated. First, a spatial averaging sampling strategy is developed. It makes use of the statistical symmetry to replace the samples from repeated experiments using ensemble methods by samples at different spatial grid points [8, 49], which greatly reduces the number of repeated experiments while the accuracy remains the same. In other words, combining the spatial averaging strategy with a small number of Monte Carlo simulation of the underlying nonlinear complex systems, the covariance with a localization structure can be sampled in an efficient and accurate fashion, at least for a finite time. Secondly, the same structure can be exploited to improve the accuracy of data assimilation [30, 37, 18, 31]. The theory developed here provides a practical guideline to choose the localization radius $L$ in a systematical way that facilitates an effective data assimilation for large dimensional systems.

The rest of the article is organized as follows. Section 2 states the main theoretical results for the spatial localization for nonlinear dynamical stochastic models and the two applications. Section 3 includes the detailed analysis of the covariance matrix with both local and nonlocal effects. Numerical simulations of a linear model and a stochastically coupled FitzHugh-Nagumo (FHN) model for excitable media in different dynamical regimes are shown in Section 4. The article is concluded in Section 5. Details of proofs are included in Appendix.

2 Problem setup and main results

To describe the behavior of nonlinear dynamical stochastic models for an excitable media, we consider a general stochastic diffusion reaction model with mean field interaction on a one dimensional torus $[0,1]$:

[TABLE]

where $W(t,x)$ is a spatial white noise. We can discretize (2.1) in space, and write down the dynamics at the $i$-th grid point $u_{i}(t)=u_{t}(\frac{i}{N})$ as

[TABLE]

where $w_{i}(t)$ are independent standard Wiener processes.

The deterministic forcings in (2.2) can be classified into two groups based on their range of action. Here $\frac{\nu}{2}N^{2}(u_{i+1}(t)+u_{i-1}(t)-2u_{i}(t))+F(u_{i}(t))$ describes how does a component $u_{i}$ force its neighbors and itself while $\frac{1}{N}\sum_{j=1}^{N}h(u_{j}(t))$ models the McKean Vlasov type of mean field interaction [39]. We will adapt this classification to simplify our notation. Moreover, we will consider the general setting where $u_{i}$ is not univariate. This will be useful when we discuss the FHN model, as we need both the membrane voltage and the recovery forcing.

Specifically, we let $\mathbf{x}(t)=(\mathbf{x}_{1}(t),\cdots,\mathbf{x}_{N}(t))$ be a stochastic process in $\mathbb{R}^{d}$. For simplicity, we assume all block share the same dimension, that is $\mathbf{x}_{i}(t)\in\mathbb{R}^{q}$ and $d=qN$. Suppose each block follows an SDE

[TABLE]

In (2.3), $\mathbf{w}_{i}(t)$ are independent standard Wiener processes of dimension $q$, and $\Sigma$ is symmetric. The index $i$ should be interpreted in a cyclic fashion, that is $\mathbf{x}_{0}=\mathbf{x}_{N},\mathbf{x}_{1}=\mathbf{x}_{N+1}$. The natural distance between indices is given by

[TABLE]

For simplicity, we assume $\{\mathbf{x}_{i}(0)\}_{i=1,\ldots,N}$ are i.i.d. samples from $\mathcal{N}(m_{0},\Sigma_{0}\Sigma^{T}_{0})$.

Many high dimensional stochastic models can be formulated as (2.3). Examples include the Lorenz 96 model [34] and the FHN model [33].

Our main results can be summarized as the follows.

Theorem 2.1**.**

Suppose the following constants are independent of $N$

[TABLE]

Here $\lambda_{max}$ denotes the maximum eigenvalue. Then for any test function $g$ and $\beta>0$, there is a constant $C_{\beta,t}$ independent on $N$, such that

[TABLE]

Here $\|\nabla g\|_{\infty}:=\sup_{\mathbf{x}}\|\nabla g(\mathbf{x})\|$. And if $\mathbf{h}\equiv 0$, the bound can be improved to $C_{\beta,t}\|\nabla g\|^{2}_{\infty}e^{-\beta\mathbf{d}(i,j)}$.

The exact value of $C_{\beta,t}$ can be found below in Theorem 3.3. The main conclusion here is that when the dimension $N$ is large, for components that are far from each other, their covariance is close to zero. In other words, the covariance between model components is significant only for nearby components. Following [41], we call such a covariance structure to be local. Below, we illustrate two important applications that show such a property is crucial for uncertainty quantification (UQ).

2.1 Sampling with spatial averaging

One fundamental UQ question regarding SDEs such as (2.3) is how to efficiently compute $\mathbb{E}g(\mathbf{x}_{i}(t))$ for a given test function and ensure that the estimation error is small. The standard Monte Carlo approach involves simulating (2.3) multiple times, and use the sample average $\hat{g}_{mc}=\frac{1}{K}\sum_{j=1}^{K}g(\mathbf{x}^{(j)}_{i}(t))$ as an estimator. However, with the increase of the dimension $N$, the computational costs to generate each independent sample $\mathbf{x}^{(j)}(t)$ increases dramatically. Therefore, only a small number of samples $K$ can be adopted in practice.

One important observation is that (2.3) is invariant under index shifting. Therefore, by symmetry, $\mathbf{x}_{i}(t)$ shares a common distribution for all $i$. As a consequence, each of them can be viewed as an individual, despite not independent, sample. Exploiting this perspective, it has been proposed in [8] the following spatial average as an estimator of $\mathbb{E}g(\mathbf{x}_{i}(t))$

[TABLE]

Since the samples are independent, the statistical properties of $\hat{g}_{M,N}$ is determined by the ones of

[TABLE]

It is straightforward to verify that $\hat{g}_{N}$ and $\hat{g}_{K,N}$ are unbiased:

[TABLE]

The accuracy of these estimators can be measured by their variance. By independence, we have $\text{var}\hat{g}_{K,N}=\frac{1}{K}\hat{g}_{N}$, which is given by

[TABLE]

Then plugging in the estimate from Theorem 2.1, the result in (2.6) can be bounded by

[TABLE]

Therefore we have

[TABLE]

In contrast, if we consider the standard estimator $\hat{g}_{mc}$, its variance is $\frac{1}{K}\text{var}(g(\mathbf{x}_{i}(t)))$, which is of $O(N)$ multiple of the one of $\text{var}\hat{g}_{M,N}$. In other words, the spatial averaging has increase the effective sample size by a factor of $N$.

2.2 Localization in data assimilation

In data assimilation (DA), sequential partial observations are incorporated to improve predictions of a stochastic dynamical system [30, 37, 18, 31]. One major application of DA is numerical weather prediction, where forecasts are obtained by combining simulations of the dynamical models with satellite data. The associated dimensions of such data assimilation problems are extremely high, often in orders of millions in order to describe variables in different temporal and spatial scales of the underlying nonlinear complex systems. Moreover, for these refined models, the current computational power can only provide around one hundred simulation samples. DA algorithms, for example the ensemble Kalman filter (EnKF), need to use these 100 samples to represent the covariance matrix of model components.

The classic sample covariance,

[TABLE]

is not accurate in this setting. A classic random matrix results [3] indicates that if $\mathbf{x}^{(i)}_{t}$ are i.i.d. samples from $\mathcal{N}(0,C)$, the covariance estimation error $\|C-\widehat{C}\|=O(\sqrt{N/K})$. Therefore, the estimator is very inaccurate in the situations with the small samples and high dimensions, i.e., $K\ll N$. The intuition behind this mathematical phenomenon is that, spurious correlation may exist between each pair of components with a small probability because of the sampling effect, but there are $O(N^{2})$ pairs of components. Thus, the resulting errors are very likely to be significant.

Nevertheless, if we are aware that covariance is local as described in Theorem 2.1, then the erroneous estimation can be made accurate. Note that the far-off-diagonal entries are close to zero. Thus, a natural way is to truncate these entries. In particular, with a given bandwidth $L$ and covariance matrix $C$, we define the localization of $C$ as

[TABLE]

It has been shown in [5] that such a localization technique can significantly improve the sampling accuracy. A simplified version in [53] shows that if $\mathbf{x}(t)$ is Gaussian distributed, then there exists a uniform constant $c$, so that for any $\epsilon>0$,

[TABLE]

In other words, if $L$ is a fixed bandwidth, it requires only $L^{2}O(\log N)$ many samples to find an accurate estimate of $C^{L}$.

On the other hand, one needs to notice that (2.7) indicates $\widehat{C}^{L}$ is a good estimator of $C^{L}$, but not necessarily $C$. In order for it to also be a good estimator of $C$, one needs $\|C-C^{L}\|$ to be small as well. Theorem 2.1 can guarantee this if there is only local interaction, i.e. $\mathbf{h}\equiv 0$. Since a covariance matrix’s $l_{2}$-norm is bounded by its $l_{1}$-norm, we can plug in the estimate in Theorem 2.1 with $g(x)=x$ and find

[TABLE]

In practice, given an error threshold $\epsilon>0$, one can combine (2.8) and (2.7) together to find a bandwidth $L$ and sample size $K=O(L^{2}\log N)$, such that the covariance estimation error $\|\widehat{C}^{L}-C\|$ is below $2\epsilon$ with high probability.

Since localization involves using very simple numerical operation to significantly improve the estimation of the covariance matrix, it has been widely applied in all DA algorithms, in particular, the ensemble Kalman filter [17, 24, 21, 57, 40, 27, 55, 54]. Moreover, adaptation of the localization strategy can be applied to Bayesian inverse problems where the associated covariance matrices follow the description of Theorem 2.1. In particular, it has been shown in [41] that blocking the high dimensional components, and then applying Gibbs sampler, the resulting Markov chain Monte Carlo algorithms will converge to the desired distribution with a rate independent of the dimension. One fundamental assumption of all these practical computational methods is that the underlying system is spatially localized, yet there has been no rigorous justification of this assumption. This article closes this gap for systems of type of (2.3).

2.3 Preliminaries

To facilitate our discussion, we will use the following notations: $\|v\|$ denotes the $l_{2}$ norm of $v$, while $\langle u,v\rangle=u^{T}v$ is the inner product; $\mathbf{I}_{q}$ denotes the identity matrix of dimension $q$. Given a matrix $A$, the norm $\|A\|$ is the $l_{2}$ operator norm of $A$; $\lambda_{max}(A)$ denotes the maximum eigenvalue of $A$ and $\|A\|_{F}$ is its Frobenius norm. We use $[A]_{i,j}$ to denote the $(i,j)$-entry of $A$.

When we have a vector $\mathbf{x}$ of $qN$ dimension, we often view it as the concatenation of $N$ sub-blocks of dimension $q$. The decomposition can be written as $\mathbf{x}=(\mathbf{x}_{1},\cdots,\mathbf{x}_{N})$. Given a $qN\times qN$ matrix $A$, we write its $(i,j)$-th $q\times q$ sub-block as $A_{i,j}$ or $\{A\}_{i,j}$. The readers should not confuse these notations with the matrix entries $[A]_{i,j}$. In particular, we can write down the entries of $A_{i,j}$:

[TABLE]

With a function $f$, $\nabla f$ denotes its Jacobian matrix. $\nabla_{\mathbf{x}_{i}}f(\mathbf{x})$ denotes the partial derivative with respect to $\mathbf{x}_{i}$, which is the $i$-th block-column of $\nabla f$.

For two symmetric matrices $A\preceq B$ indicates that $B-A$ is positive semidefinite.

3 Analysis of covariance matrices

We provide the explicit statements of our result in this section. The detailed proofs are allocated in Appendix.

3.1 Linear covariance propagation

Consider the following multivariate SDE:

[TABLE]

Our main result is established through a covariance comparison principle between (3.1) and a linear surrogate.

Theorem 3.1**.**

Suppose the covariance propagation of (3.1) is dominated by $(F_{*}(t),H_{*}(t))\in(\mathbb{R}^{N\times N},\mathbb{R}^{N\times N})$ in the following sense:

1. (1)

$\nabla_{\mathbf{x}_{i}}\mathbf{f}_{i}(t,\mathbf{x}_{i-1},\mathbf{x}_{i},\mathbf{x}_{i+1})+\nabla_{\mathbf{x}_{i}}\mathbf{f}_{i}(t,\mathbf{x}_{i-1},\mathbf{x}_{i},\mathbf{x}_{i+1})^{T}\preceq 2[F_{*}(t)]_{i,i}\mathbf{I}$. 2. (2)

$\|\nabla_{\mathbf{x}_{i-1}}\mathbf{f}_{i}(t,\mathbf{x}_{i-1},\mathbf{x}_{i},\mathbf{x}_{i+1})\|\leq[F_{*}(t)]_{i,i-1},\|\nabla_{\mathbf{x}_{i+1}}\mathbf{f}_{i}(t,\mathbf{x}_{i-1},\mathbf{x}_{i},\mathbf{x}_{i+1})\|\leq[F_{*}(t)]_{i,i+1}$. 3. (3)

$\frac{1}{N}\|\nabla_{\mathbf{x}_{j}}\mathbf{h}_{i}(t,\mathbf{x}_{j})\|\leq[H_{*}(t)]_{i,j}$.

Let $G_{*}(s)=F_{*}(s)+H_{*}(s)$, and

[TABLE]

Then for any test function $g$ with bounded gradient,

[TABLE]

This theorem shows that if the covariance propagation of (3.1) is dominated by $G_{*}=F_{*}+H_{*}$, then the spatial covariance of (3.1) is also dominated by the linear surrogate model associated with $G_{*}$:

[TABLE]

To see that (3.3) has the spatial covariance structure as in (3.2), we concatenate (3.3) into one equation for $\mathbf{x}(t)$

[TABLE]

Here, $D_{G_{*}}$,$D_{\Sigma}$ and $D_{\Sigma_{0}}$ are $qN\times qN$ dimensional matrices. Their $q\times q$ sub-blocks are given by

[TABLE]

Since this equation is linear, we can apply the Duhamel’s formula and write the solution as

[TABLE]

By Ito’s isometry,

[TABLE]

In particular, the covariance between $\mathbf{x}_{i}(t)$ and $\mathbf{x}_{j}(t)$ is given by $\{\text{cov}\,\mathbf{x}(t)\}_{i,j}$. Since $D_{\Sigma_{0}}$ and $D_{\Sigma}$ are block diagonal, it is not difficult to establish (3.2) for system (3.3).

3.2 A representation of the covariance matrix

One difficulty in reaching Theorem 3.1 is how to compute the covariance matrix for a multivariate SDE. We will establish a framework for covariance matrix computation, starting from a simple Stein’s lemma type of formula:

Lemma 3.2**.**

Suppose $X,Y$ are independent standard Gaussian random variables of dimension $d$, $f$ and $g$ are smooth functions with bounded first order derivatives. Let $X^{\theta}=\cos\theta\,X+\sin\theta\,Y$, then

[TABLE]

In order to apply Lemma 3.2, we need to represent $\mathbf{x}(t)$ as a function of certain Gaussian random variables. These Gaussian random variables are the increments of the Wiener processes driving $\mathbf{x}(t)$. To setup this connection, recall that one way to simulate (3.1) is to perturb its ordinary differential equation (ODE) version with random noises. In particular, consider the ODE system

[TABLE]

Let $\Psi_{h}:\mathbb{R}^{qN}\mapsto\mathbb{R}^{qN}$ be the mapping from $\mathbf{x}(t)$ to $\mathbf{x}(t+h)$. We suppress the dependence of $\Psi_{h}$ on $t$, for notational simplicity. Then one way to simulate (3.1) is generate the following iterates

[TABLE]

Then $\mathbf{X}(T)$ is a numerical approximation of $\mathbf{x}(t)$ if $Th=t$ when $T$ is large. A standard verification Lemma of this statement is provided in Appendix. We write the $i$-th block of $\mathbf{X}(n)$ as $\mathbf{X}_{i}(n)\in\mathbb{R}^{q}$. $\mathbf{W}_{n}$ are $qN$ dimensional standard Gaussian noise for $n=0,1,\cdots,T$.

The strategy here is to consider $\mathbf{X}(T)$ instead of $\mathbf{x}(t)$, because its dependence on the noise realization $\mathbf{W}_{i}$ is easier to handle using Lemma 3.2. We can show a similar result for $\mathbf{X}(T)$ as in Theorem 3.1. Then in the limit $h\to 0$, we have our desired results. The details are shown in Appendix. As a remark, it might be possible and more elegant to represent the covariance of $\mathbf{x}_{t}$ directly using tools from Malliavin calculus [43], and avoid the procedures of discretization and taking continuous limit. Yet the current proof strategy can be easier to understand for an applied math audience.

3.3 Diffusion and mean field interaction

While (3.2) provides a general upper bound, from its formulation is not easy to see the covariance decay we expect. To make the covariance decay explicit, we try to use instead a homogenous linear surrogate. We update the constant definition in (2.4) to the spatial inhomogeneous setting:

[TABLE]

They are finite under the conditions of Theorem 3.1. We can build a homogeneous linear surrogate by letting

[TABLE]

Then the linear surrogate (3.3) corresponds

[TABLE]

where the discrete Laplacian is given by

[TABLE]

which is closely related to the heat equation.

The spatial covariance of $\mathbf{x}(t)$ is built up by the diffusion effect of $\Delta\mathbf{x}_{i}$ and the mixing effect of mean field interaction $\frac{1}{N}\sum_{j=1}^{N}\mathbf{x}_{j}$. These two effects operate in different ways. The diffusion propagates the diagonal entries in the covariance matrix to the off-diagonal ones, it contributes to the exponential decay term $e^{-\beta\mathbf{d}(i,j)}$ in Theorem 2.1. The mean field interaction synchronizes the components, and contributes to a global covariance of scale $\frac{1}{N}$.

We can establish the following explicit bounds in terms of $\lambda_{G},\lambda_{F},\lambda_{H}$, instead of the matrix $Q(t)$:

Theorem 3.3**.**

Under the same conditions of Theorem 3.1, for any $\beta>0$, let

[TABLE]

with $\lambda_{0},\lambda_{F}$ and $\lambda_{H}$ defined by (3.6). Then with any function $g$,

[TABLE]

When there is only mean field interaction, i.e. $\lambda_{F}=0$, the result can be simplified by choosing $\beta\to\infty$,

[TABLE]

When there is only diffusion effect, that is $\mathbf{h}_{i}\equiv 0$, the result can be simplified as

[TABLE]

It is worth noticing the result holds for all $\beta>0$. This indicates the covariance structure actually decays faster than exponential. In practice, one can try to minimize the right hand side with respect to $\beta$, or directly compute the covariance of (3.7) to find the exact local covariance structure. While sharper theoretical upper bound maybe obtainable by applying formulas as in [28], they are complicated and work only asymptotically. Our current upper bound is simpler, and sufficient for applications discussed in Sections 2.1 and 2.2.

With Theorem 3.3, it is straightforward to show the spatial averaging estimator is consistent:

Corollary 3.4**.**

Under the same settings of Theorem 3.3, the following bound holds for the block estimator (2.5):

[TABLE]

Proof.

Since in (2.6), $\sum^{N}_{j=1}e^{-\beta\mathbf{d}(1,j)}\leq 2\sum^{\infty}_{j=0}e^{-\beta j}=\frac{2}{1-e^{-\beta}}$. ∎

3.4 Long time stability

In many applications, the long time $t\to\infty$ scenario is of interest. For example, if we want to find the equilibrium measure of $\mathbf{x}(t)$, we often simulate $\mathbf{x}(t)$ for a long time, and use its distribution as an approximate. Likewise, if we are interested in data assimilation with models as (2.3), the algorithms are usually iterated for an extended period. In order to apply the spatial averaging strategy in Section 2.1, and the localization strategy in Section 2.2 for these operations, we need the local covariance structure to be stable in time.

In the view of Theorem 3.3, we can guarantee the local covariance structure is stable if $\lambda_{0}+2\lambda_{F}+\lambda_{H}<0$. In particular we have the following corollary by letting $t\to\infty$ in the estimate.

Corollary 3.5**.**

*If $\lambda_{0}+\lambda_{H}+2\lambda_{F}<0$, then there is a $\beta_{0}>0$, so that if $0<\beta<\beta_{0}$, $\lambda_{\beta}\leq\eta_{\beta}<0$ we have *

[TABLE]

4 Test Examples

In this section, we demonstrate our theoretical findings with several numerical test examples.

4.1 A linear model

We start with a linear model, where the analytic solutions associated with this linear model can be written down explicitly (See Appendix B). Therefore, this model is able to provide insights and validation of the theories developed in the previous sections.

The linear model is as follows:

[TABLE]

with periodic boundary conditions. Here $\bar{u}$ is the averaged value of all the state variables,

[TABLE]

In the linear model (4.1), the coefficients of damping $a$, diffusion $d_{u}$, mean field interaction $w$ and the stochastic noise $\sigma_{u}$ are all constants. Thus, the flow field is statistically homogeneous. Below, we fix

[TABLE]

while the diffusion and mean field interaction coefficients $d_{u}$ and $w$ vary in different cases.

In the linear model (4.1), it is easy to write down the constants in (3.6),

[TABLE]

In addition, only one-layer model is utilized here and we consider the covariance of different $u_{i}$ themselves. Therefore,

[TABLE]

Furthermore, $\|\Sigma^{2}\|_{F}=\sigma_{u}^{2}$. Below, we study the role of the mean field interaction and the diffusion. All the results shown in this subsection are computed using the exact solution given in Appendix B, which allows us to understand Theorem 3.3 without being interfered by numerical errors.

1. Linear model with only mean field interaction. First, we study the case with only mean field interaction. In other words, we set the diffusion coefficient $d_{u}=0$. The constants in (4.3) reduce to

[TABLE]

According to (3.8) in Theorem 3.3, the covariance bound is proportional to $1/N$ and the coefficient is provided by letting $\beta\to\infty$. Therefore, assuming the initial covariance is zero, we have the following result

[TABLE]

In Figure 4.1, we show the covariance between $u_{1}$ and $u_{i}$ for $i=1,\ldots,N$ at $t=5$. Here the initial conditions are zero everywhere. Different columns show the covariance with different $N$. The spatiotemporal simulations are also included to provide intuitions. Without the diffusion and the initial covariance, the covariance between $u_{1}$ and all other $u_{i}$ with $i\neq 1$ has the same non-zero value due to the global effect from the mean field interaction.

In Figure 4.2, the dependence of the covariance between $u_{1}$ and $u_{2}$ as a function of $N$ is shown. It is clear that the covariance decays as a function of $1/N$, which validates the results in Theorem 3.3 (and that in (4.6)). The black dashed line in panel (b) (log-log scale) shows the bound in (4.6), where the constant $C_{1}$ is given by the coefficient on the right hand side of (4.6) in front of $1/N$.

2. Linear model with only diffusion. Next, we study the case with only diffusion. In other words, we set the mean field interaction coefficient $w=0$. The constants in (4.3) reduce to

[TABLE]

According to (3.9) in Theorem 3.3, the covariance bound is proportional to $e^{-\beta\mathbf{d}(i,j)}$. Again, assuming the initial covariance is zero, we have the following result

[TABLE]

where

[TABLE]

In Figure 4.3, we show the covariance between $u_{1}$ and $u_{i}$ for $i=1,\ldots,N$ at $t=5$ with $N=64$ fixed. Here the initial conditions are zero everywhere. Different columns show the covariance with different strengths of the diffusion coefficient $d_{u}$. The spatiotemporal simulations are also included to provide intuitions. Without the mean field interaction and the initial covariance, the covariance between $u_{1}$ and $u_{i}$ for $u_{i}$ being far from $u_{1}$ decays to zero, reflecting the local contribution due to only the diffusion (without mean field interaction). On the other hand, the increase of $d_{u}$ allows the increase of the covariance between $u_{1}$ and $u_{i}$ for $i$ far from $1$. This is clearly illustrated in the spatiotemporal patterns and is consistent with our intuition as well.

Figure 4.4 shows the covariance between $u_{1}$ and $u_{i}$ for $i=1,\ldots,N/2$ with $d_{u}=20$ and $N=64$ at $t=5$. In panel (b), the logarithm scale is shown, which indicates a linear dependence of the log covariance on the distance $\mathbf{d}(1,i)$. This validates the exponential decay of the covariance as a function of the distance as shown in (4.8). The black dashed line shows the bound in (4.8), where the constant $C_{2}$ is given by the coefficient on the right hand side of (4.8) in front of $e^{-\beta\mathbf{d}(i,j)}$ by taking $\beta=1/5$. We have tested other $\beta$, with which the bounds all above the actual covariance. Therefore, we validate the theoretical results with only diffusion, described by (3.9).

3. Linear model with both the mean field interaction and diffusion. Finally, we study the case with both mean field interaction and diffusion. Figure 4.5 is similar to Figure 4.3 but with a non-zero mean field interaction $w=5$. Therefore, unlike the results in Figure 4.5, the covariance between $u_{1}$ and $u_{i}$ where $i$ is far from $1$ is non-zero due to the mean field interaction that has a global impact. In Figure 4.6, it is also clear that when $i$ increases the covariance between $u_{i}$ and $u_{1}$ first experience an exponential decay and then the covariance remains as a constant. Thus, the exponential decay is the main contribution to the covariance behavior for grid points that are close to $i=1$ while the mean field interaction plays the dominant role at the location that is far away. These simulations are consistent with the theoretical results shown in Theorem 3.3.

4.2 A stochastically coupled FHN model

In this subsection, we explore the covariance structure in a stochastically coupled FitzHugh-Nagumo (FHN) model, which is a nonlinear model that describes activation and deactivation dynamics of spiking neurons [33, 8]. The model reads,

[TABLE]

where $u_{i}$ and $v_{i}$ are activator and inhibitor variables. Periodic boundary conditions are imposed on $u_{i}$ variables. In (4.9), the timescale ratio $\epsilon=0.01\ll 1$ leads to a slow–fast structure of the model. The parameter $a=1.05>1$ such that the system has a global attractor in the absence of noise and diffusion [20]. The random noise is able to drive the system above the threshold level of global stability and triggers limit cycles intermittently. Note that with $N=1$, the model reduces to the classical FHN model with a single neuron and it contains the model families with both coherence resonance and self-induced stochastic resonance [16]. With different choices of the noise strength $\delta_{1},\delta_{2}$ and the diffusion coefficient $d_{u}$, the system in (4.9) exhibits rich dynamical behaviors. Finally, in (4.9), the variable $\bar{u}=\frac{1}{N}\sum_{i=1}^{N}u_{i}$ is the mean field interaction.

The goal of this subsection is to study the spatial structure of the covariance due to both the local effect (i.e., diffusion) and the global effect (i.e., mean field interaction).

4.2.1 Basic properties

For the convenience of stating the theoretical results, we slightly change of notation of (4.9) in this paragraph,

[TABLE]

To apply the framework we developed, it is natural to let $\mathbf{x}_{i}=(u_{i},v_{i}^{\epsilon})^{T}$ such that $q=2$, and

[TABLE]

and

[TABLE]

Since the $-\epsilon^{-1}u_{i}^{2}$ is always negative, we can take

[TABLE]

Therefore, using the spatial averaging sampling strategy for high dimensional systems, we have

Corollary 4.1**.**

The block average estimator is consistent for the FHN system (4.10). In particular, estimate (3.10) holds with $\|\Sigma^{2}\|_{F}=\sqrt{\delta_{1}^{4}+\delta^{4}_{2}}/\epsilon$, and

[TABLE]

In particular, from the bound (3.10), we learn that the estimator consistency deteriorates with larger $t$, $d_{u}$ and $\tilde{w}$, which correspond to longer simulation time, stronger local interaction, and stronger mean field interaction regimes.

4.2.2 The FHN model with different diffusion and mean field interaction coefficients

1. FHN model with only diffusion. First, we focus on the situation in the FHN model where the mean field interaction is zero ($w=0$). Therefore, the covariance comes purely from the local effect, namely the diffusion. We set $\delta_{1}=\delta_{2}=0.4$ but choose different values of the diffusion parameter $d_{u}$ which provides different dynamical behavior:

[TABLE]

The simulation of the model in these three dynamical regimes is shown in Figure 4.7. Figure 4.8 shows the covariance between $u_{1}$ and $u_{i}$ for $i=1,\ldots,N$ at $t=5$. Here we use a Monte Carlo simulation with $M=8192$ samples at each grid points (namely, running the Monte Carlo for $M$ times). With the increase of $d_{u}$ from Regime (a) to Regime (c), the spatial covariance of $u$ increases. In particular, in Regime (c) with a strong diffusion, the covariance between $u_{1}$ and all other $u_{i}$ for $i=1,\ldots,N$ remains significant, and is above $0.5$. Note that there is almost no mixing of the underlying stochastic system with such a large diffusion coefficient in Regime (c) (See Figure 4.7), which may cause some issues of the controllability of the system. In fact, with a large diffusion coefficient, the role of the stochastic forcing becomes relatively weaker and the control of the system via the external stochastic forcing becomes more difficult. This also leads to some potential issues in applying the spatial averaging strategy for computing the covariance, especially for long time. Note that although Corollary 4.1 is always valid, the constant $\eta_{\beta}$ becomes very large with a large diffusion coefficient (due to the second term) and thus the consistency becomes nearly degenerated very quickly. See the remarks below Corollary 4.1 and the statement of the bound in (3.10). The numerical results using the spatial averaging strategy will be illustrated in the next subsection. On the other hand, the spatial covariance of $v$ changes in the same fashion as $u$ but is weaker. This is because the diffusion appears in the $u$ equations and it only has an indirect impact on $v$. In all three regimes, the spatial covariance decays in an exponential rate (see the subpanels). All these are consistent with the result in (3.9) of Theorem 3.3.

2. FHN model with only mean field interaction. Now we study the situation in the FHN model where the diffusion is zero ($d_{u}=0$). Therefore, the covariance comes purely from the global effect (i.e., the mean field interaction). We set again $\delta_{1}=\delta_{2}=0.4$ but choose different values of the mean field interaction parameter $w$ which provides different dynamical behavior:

[TABLE]

The simulations of the three regimes are shown between Figure 4.9. Again, these results are computed using a Monte Carlo simulation with $M=8192$ samples at each grid points. In the moderate mean field interaction regime (Regime (b)), some global effects can already been seen. In the strong mean field interaction regime (Regime (c)), it is clear that the variables at different grid points are highly synchronized with each other, which is consistent with the theoretic analysis described in the previous section.

Note the difference in Figure 4.7 and 4.9. In Figure 4.7 the diffusion leads to the development of local interaction and extends it to the global scale. In Figure 4.9 the mean field interaction starts directly with a change in the global scale and such a change becomes significant with the increase of the mean field interaction strength $w$. Again, the constant $\eta_{\beta}$ in Corollary 4.1 with a large mean field interaction coefficient can become very large (due to the first term) and thus the consistency becomes nearly degenerated.

Figure 4.10 shows the covariance $\mbox{Cov}(u_{1},u_{2})$ and $\mbox{Cov}(v_{1},v_{2})$ as a function of $N$ in two different regimes with either moderate or strong mean field interactions. With the increase of $N$, the covariance decreases and the decrease rate is proportion to $1/N$ (see the log-log plot), which is consistent with the theoretical conclusion in (3.8) of Theorem 3.3.

4.2.3 Recovery of the covariance using the spatial averaging strategy

In this subsection, we compare the covariance computed from the spatial averaging strategy with the truth in order to study the consistency of the spatial averaging strategy. Here the truth is formed by using direct Monte Carlo simulations.

Recall in Section 2.1, the spatial averaging strategy makes use of the statistical symmetry of the system. Comparing with the covariance computed by using the direct Monte Carlo simulations with a large sample size, the spatial averaging strategy requires only a small number of repeated experiments to reach the same accuracy (providing that the consistency is valid). In fact, due to the spatial homogeneity, the covariance between $u_{1}$ and $u_{2}$ is the same as that between $u_{i}$ and $u_{i+1}$ for all $i=1,\ldots,N$. Therefore, regarding $u_{1},u_{2},\ldots,u_{N}$ as a sample set, the covariance between $u_{1}$ and $u_{2}$ can be computed by taking the covariance between $u_{1},u_{2},\ldots,u_{N}$ and its shifted sample set $u_{2},u_{3},\ldots,u_{N},u_{1}$. Note that if $M$ repeat experiments are combined (namely, $M$ Monte Carlo simulations) with the spatial averaging strategy, then the effective sample size is $M_{eff}=MN$. One potential issue of the spatial averaging strategy is that the samples $u_{1},u_{2},\ldots,u_{N}$ are not independent with each other. Therefore, when these samples are strongly correlated, then the spatial averaging strategy in computing the covariance may lead to biases. In this subsection, we compare the covariance computed from the spatial averaging strategy with the truth that is generated from the Monte Carlo simulations in different dynamical regimes. The goal is to explore the consistency of this spatial averaging strategy in computing the covariance in different dynamical regimes.

Below, we always take $N=512$ spatial grid points. For a fair comparison, we generate the truth using the direct Monte Carlo simulation with $M_{MC}=N=512$ samples, which equals to the effective sample number using the spatial averaging strategy with simply one simulation. In all the tests shown below, we compare the covariance computed using the above two strategies at time $t=0.5,t=1,t=2$ and $t=5$. The regimes to be tested contain both local and non-local effects with various strengths. See Figures 4.11–4.12 for the spatiotemporal simulations of different dynamical regimes.

In Figure 4.11, strong spatial correlations are seen in all three dynamical regimes and the underlying stochastic processes is not mixing. Yet, there are some fundamental difference in these dynamical regimes. In Regime (a), there is a strong diffusion but there is no mean field interaction. Thus, the correlation comes from only the local effect but the diffusion is so strong that each grid point affects the other points that are far from itself. In Regime (b), only the mean field interaction plays an role since the diffusion is zero but the mean field interaction coefficient is quite large. Thus, a strong synchronization is seen in the spatiotemporal patterns. Regime (c) involves both a moderate diffusion and a moderate mean field interaction. The combined local and non-local effects also contribute to the strong correlation in the spatiotemporal patterns. Figure 4.13 shows the calculated covariance based on the spatial averaging strategy (red) and the truth generated from the direct Monte Carlo simulation (green). For a short time up to $t=1$, the consistency is found in Regimes (a) and (b). But the covariance computed from the spatial averaging strategy becomes significantly different from the truth after $t=2$, which indicates the loss of the long-term consistency using the spatial averaging strategy. These results are related to the loss of practical controllability as was discussed in the previous subsection. In fact, the constant $\eta_{\beta}$ in Corollary 4.1 becomes very large in both regimes due to the strong diffusion and/or the strong mean field interaction and thus the consistency of the spatial averaging strategy becomes nearly degenerated very quickly, which has been pointed out in the remarks below Corollary 4.1. For Regime (b), due to the strong global effect (via mean field interaction), the consistency is again lost very quickly at $t=1$. Note in all these three regimes the covariance is not localized after a certain time due to either strong diffusion or strong mean field interaction (or both). Thus, the spatial averaging strategy does not work well for long time.

On the other hand, no strong spatial synchronization is observed in the three dynamical regimes shown in Figure 4.12. Some weak spatial synchronization is observed in Regimes (d) and (e), where the former involves a moderate diffusion and a weak mean field interaction while the latter is given by a moderate mean field interaction but no diffusion. Regarding the spatial covariance, Figure 4.14 shows that the spatial averaging strategy is able to provide accurate solutions up to $t=2$. At $t=5$, a slight inconsistency in the solution computed from the spatial averaging strategy is found, which is due to the degeneration of the consistency for long time as indicated in Corollary 4.1. Thus, the spatial averaging strategy is skillful at short and moderate times but the long-term consistency cannot be guaranteed in these regimes.

Finally, in Regime (f), which involves only a moderate diffusion but no mean field interaction, the consistency of the solution using the spatial averaging strategy lasts longer and is at least up to $t=5$. See Figure 4.15.

To summarize, both a strong global effect (mean field interaction) and a strong local effect (diffusion) will lead to the loss of the long-term consistency in the solutions using the spatial averaging strategy, which has been clearly indicated in Corollary 4.1 and the remarks below it. The spatial averaging strategy is skillful and can be applied for a relatively long time when the FHN model has a weak global or a moderate local effect.

5 Conclusions

In this article, a theoretical framework to understand the spatial localization for a large class of stochastically coupled nonlinear systems in high dimensions is developed and is applied to nonlinear complex models for excitable media. Rigorous mathematical theories (Theorem 3.3) show covariance decay behavior due to both local and nonlocal effects, which result from the diffusion and the mean field interaction, respectively. The analysis is based on a comparison with an appropriate linear surrogate model, of which the covariance propagation can be computed explicitly. Two important applications of these theoretical results are illustrated. They are the spatial averaging strategy for efficiently sampling the state variables and the localization technique in data assimilation. Test examples of a linear model and a stochastically coupled FitzHugh-Nagumo model for excitable media are adopted to validate the theoretical results. The later is also used for a systematical study of the spatial averaging strategy in efficiently sampling the covariance matrix in different dynamical regimes. The spatial averaging strategy is always skillful for short time. The long-term consistency of the spatial averaging strategy requires that the mean field interaction and the diffusion cannot be too strong and the covariance matrix has localized structures.

Acknowledgement

The research of N.C. is supported by the Office of Vice Chancellor for Research and Graduate Education (VCRGE) at University of Wisconsin-Madison. The research of A.J.M. is partially supported by the Office of Naval Research Grant ONR MURI N00014-16-1-2161 and the Center for Prototype Climate Modeling (CPCM) at New York University Abu Dhabi Research Institute. X.T.T is supported by NUS grant R-146-000-226-133.

Appendix A Proof

We provide the proofs of the theoretical results in this appendix. In particular for Theorem 3.1, its arguments are broken into subsections

A.1 Consistency of the numerical formulation

Lemma A.1**.**

Under the conditions of Theorem 3.3, let $\mathbf{X}(0)=\mathbf{x}(0)$, and $\mathbf{W}_{n}=\mathbf{w}_{hn}-\mathbf{w}_{h(n-1)}$, and consider the corresponding $X(n)$ given by (3.5). Then if we write $\Delta(t)=\mathbf{X}(T)-\mathbf{x}(t)$, there is a constant $D_{t}$ such that for any sufficiently small $h$,

[TABLE]

Proof.

Let $\mathbf{z}(nh)=\mathbf{X}(n)$, and $\mathbf{z}((n+s)h)=\Psi_{sh}(\mathbf{X}(n))$ for $0\leq s<1$. Clearly $\mathbf{z}(0)=\mathbf{x}(0)$. Note that

[TABLE]

So

[TABLE]

Note that under our conditions

[TABLE]

But

[TABLE]

Therefore

[TABLE]

So by Gronwall’s inequality

[TABLE]

Next notice that

[TABLE]

By Ito’s isometry

[TABLE]

Therefore

[TABLE]

Therefore, when $h$ is small enough, by Taylor expansion,

[TABLE]

Then by the discrete Gronwall’s inequality, if $\lambda_{0}+2\lambda_{H}>0$, at $Th=t$, since $\exp(a)-1>a$ for all $a>0$,

[TABLE]

∎

A.2 Decomposition of covariance

Clearly, we can view $\mathbf{X}(T)$ as a function of the realization of $\mathbf{W}_{i}$,

[TABLE]

For the simplicity of notation, we write $\mathbf{W}_{0},\mathbf{W}_{1},\mathbf{W}_{2},\cdots,\mathbf{W}_{T}$ as $\mathbf{W}_{0:T}$ in below. And the $i$-th block of $\Phi\in\mathbb{R}^{qN}$ as $\Phi_{i}\in\mathbb{R}^{q}$. Then we can write $\mathbf{X}_{i}(T)=\Phi_{i}(\mathbf{W}_{0:T})$.

We can study the dependence of $\mathbf{X}(T)$ on the $\mathbf{W}_{i}$ through the gradients, written as $\nabla_{i}\Phi:=\nabla_{\mathbf{W}_{i}}\Phi$. By the chain rule, it has this back propagation formulation,

[TABLE]

To investigate the covariance, we consider another independent sample

[TABLE]

Then clearly the covariance can be written as

[TABLE]

We can decompose it by considering

[TABLE]

Because $\mathbf{X}(T)=\Phi(\mathbf{W}_{0:T}),\mathbf{X}^{\prime}(T)=\Phi(\mathbf{W}^{\prime}_{0:T})$, so

[TABLE]

We will bound $\mathbb{E}\Delta_{k}$ using a Stein’s lemma type of formula, Lemma 3.2. Its proof is given below.

A.3 Stein’s lemma

Proof of Lemma 3.2.

The first part comes simply by

[TABLE]

The note that

[TABLE]

Note that when $X\sim\mathcal{N}(0,\mathbf{I}_{d})$, and smooth vector field $h:\mathbb{R}^{d}\to\mathbb{R}^{d}$, let $c=(2\pi)^{-\frac{d}{2}}$, by the divergence theorem of calculus, $\int\nabla\cdot(h(\mathbf{x})e^{-\frac{1}{2}|\mathbf{x}|^{2}})d\mathbf{x}=0$, therefore

[TABLE]

Here, $\nabla\cdot$ denotes the divergence operator. Therefore, if we condition on the value of $Y$ and use law of total expectation

[TABLE]

Likewise, we have the following

[TABLE]

Therefore, we have

[TABLE]

and this finishes our proof. ∎

To apply Lemma 3.2, we view $\mathbf{W}_{k}$ as $X$, and $\mathbf{W}^{\prime}_{k}$ as $Y$ in (1.2). Let

[TABLE]

Then

[TABLE]

Note that by the chain rule,

[TABLE]

So by Lemma 3.2, we have

[TABLE]

Note that $\nabla_{k}\Phi_{i}(\mathbf{W}_{0:T})\nabla_{k}\Phi_{j}(\mathbf{W}^{\theta,k})^{T}\in\mathbb{R}^{q\times q}$ is the $(i,j)$-th sub-block of

[TABLE]

So essentially we want bound (1.5) in each of its block. To do that, we need to consider the derivative of the flow $\Psi_{t}$.

A.4 Derivative flow

Lemma A.2**.**

Fix any realization of $\mathbf{W}_{0:T}$ and $\mathbf{W}_{0:T}^{\prime}$. Let $Q_{h}(s)=\exp(hG_{*}s)\exp(hG_{*}s)^{T}$. The $(i,j)$-th sub-block of $R^{k}$ is bounded by

[TABLE]

Proof.

We will discuss only the $k\geq 1$ case, the analysis $k=0$ case is identical. For any $s\leq T-k$, we can extend the definition of $\mathbf{X}(t)$ to noninteger $t$, by letting

[TABLE]

Here $\lfloor s\rfloor$ denotes the largest integer that is less than $s$. The derivative flow of $\Psi_{t}$ can also be extended. We define the following $\mathbb{R}^{qN\times qN}$ matrix:

[TABLE]

Then following the back propagation formula (1.1), $\nabla_{k}\Phi=P(T-k)$. $P(s)$ also follows the following differential equation with initial condition $P(0)=\Sigma_{h}$,

[TABLE]

$G(s,\mathbf{x})\in\mathbb{R}^{qN\times qN}$ is the Jacobian created by the ODE flow (3.4). So its $(i,j)$-th sub-block is given by

[TABLE]

Clearly, it is a block-tridiagonal matrix.

Because of (1.5), we are also interested in $\nabla_{k}\Phi(\mathbf{W}^{\theta,k})$. We consider a numerical sequence generated by (3.5), but the driving noise is $\mathbf{W}^{\theta,k}$. Fix a $\theta\in[0,\frac{\pi}{2}]$. Let $\mathbf{Y}(n)=\mathbf{X}(n)$ for $n=0,\ldots,k$,

[TABLE]

and $\mathbf{Y}(n+1)=\Psi_{h}(\mathbf{Y}(n))+\Sigma_{h}\mathbf{W}^{\prime}_{n}$ for $n\geq k+1$. We also define $\mathbf{Y}(k+s)=\Psi_{sh}(\mathbf{Y}(k+\lfloor s\rfloor))$ for noninteger time $s\neq\lfloor s\rfloor$. Define also the following matrix

[TABLE]

It follows the ODE

[TABLE]

By the back propogation formula (1.1), $P^{\prime}(T-k)=\nabla_{k}\Phi(\mathbf{W}^{\theta,k})$.

Let $R(s)=P(s)(P^{\prime}(s))^{T}$. The quantity we are interested is $R^{k}=R(T-k)$. $R(s)$ follows the following ODE:

[TABLE]

Here and below, we write $G(k+s,\mathbf{X}_{k+s})$ as $G(\mathbf{X}_{k+s})$ for simplicity. Consider the $(i,j)$-th sub block of $R_{i,j}(s)$. In below, we use $\{AB\}_{i,j}$ to denote the $(i,j)$-th sub block of matrix $AB$. We also define the inner product of matrices as $\langle A,B\rangle=\text{tr}(A^{T}B)$. Then $\|R_{i,j}(s)\|^{2}_{F}=\langle R_{i,j}(s),R_{i,j}(s)\rangle$. We take time derivative, and recall that $G$ is block tridiagonal, so

[TABLE]

Because of the covariance propagation assumption in Theorem 3.1, and an elementary A.3 for matrix inner product, we can bound each term as below

[TABLE]

[TABLE]

[TABLE]

[TABLE]

In conclusion, we have the following with the dependence of $G_{*}$ on $s+k$ suppressed,

[TABLE]

Since $\frac{d}{ds}\|R_{i,j}(s)\|^{2}_{F}=2\|R_{i,j}(s)\|_{F}\frac{d}{ds}\|R_{i,j}(s)\|_{F}$, we have

[TABLE]

Next we consider $Q_{h}(s)$. Note that $\|h\Sigma^{2}\|_{F}[Q_{h}(0)]_{i,j}=\|R_{i,j}(0)\|_{F}$. Moreover, $Q_{h}$ follows the ODE

[TABLE]

Use the fact that $G_{*}$ is tridiagonal,

[TABLE]

If we view $[Q_{h}(s)]_{i,j}$ as a linear ODE, the non diagonal term of this ODE system , $[G_{*}]_{i,i+m}$ and $[G_{*}]_{j,j+m}$, are all nonnegative. So by Lemma A.4, we have that

[TABLE]

This leads to our claims. ∎

Lemma A.3**.**

Suppose $A+A^{T}\preceq 2\lambda_{A}I$, then $\langle AX,X\rangle\leq\lambda_{A}\langle X,X\rangle$, and $\langle XA,X\rangle\leq\lambda_{A}\langle X,X\rangle$.

Proof.

Let $v_{i}$ be the right eigenvectors of $X$, so $Xv_{i}=\lambda_{i}v_{i}$, then

[TABLE]

For the second claim, let $v_{i}$ be the right eigenvectors of $X^{T}$, so $X^{T}v_{i}=\lambda_{i}v_{i}$

[TABLE]

∎

Lemma A.4**.**

If $A\in\mathbb{R}^{d\times d}$ has all non-diagonal entries being nonnegative, then $\exp(sA)$ has all entries being nonnegative for any $s$. In particular, if $X(s)\in\mathbb{R}^{d}$ and $Y(s)\in\mathbb{R}^{d}$ satisfies

[TABLE]

where $\leq$ is interpreted entry-wise. Then $X(0)\leq Y(0)$ leads to $X(s)\leq Y(s)$ entry-wise.

Proof.

The statement is elementary if all of $A$ diagonal entries are also nonnegative, since

[TABLE]

On the other hand, we can always find a $\lambda$ so that $\lambda\mathbf{I}+sA$ has all diagonal entries being nonnegative. Because $\mathbf{I}$ commute with $sA$, so

[TABLE]

and all entries of $\exp(sA)$ are nonnegative.

For the ODE system, we first consider the case when all entries of $A$ are positive. Define $\Delta(s)=Y(s)-X(s)$. $\Delta(0)\geq\mathbf{0}$ entry-wise. Also when $\Delta(s)\geq 0$,

[TABLE]

meaning $\Delta(s)$ is increasing entry-wise. Therefore $\Delta(s)\geq\mathbf{0}$ entry-wise for $s\geq 0$.

If the diagonal terms of $A$ are negative, we can consider $X_{\lambda}(t)=e^{\lambda t}X(t)$ and $Y_{\lambda}(t)=e^{\lambda t}Y(t)$, then

[TABLE]

By our previous analysis, $X_{\lambda}(t)\leq Y_{\lambda}(t)$ entry-wise, which leads to $X(t)\leq Y(t)$ entry-wise. ∎

A.5 Summarizing argument

Proof of Theorem 3.1.

Recall in (1.4), we have

[TABLE]

Note that $R^{k}_{i,j}$ is $q\times q$ dimensional, so by Cauchy-Schwartz, $|\text{tr}(R^{k}_{i,j})|\leq\sqrt{q}\|R^{k}_{i,j}\|_{F}.$ Lemma A.2 shows that

[TABLE]

Therefore for $k\geq 1$,

[TABLE]

and likewise $|\mathbb{E}\Delta_{0}|\leq\sqrt{q}\|\nabla g\|^{2}_{\infty}\|\Sigma^{2}_{0}\|_{F}[Q_{h}(T)]_{i,j}$. Consequentially, we have

[TABLE]

When we let $h\to 0$, $\mathbf{X}(T)\to\mathbf{x}(t)$ by Lemma A.1. Then

[TABLE]

By the Lipschitzness of SDE coefficient, it is standard to show that $\mathbb{E}|g(\mathbf{x}_{i}(t))|^{2}\leq\|\nabla g\|_{\infty}^{2}\mathbb{E}\|\mathbf{x}_{i}(t))\|^{2}$ is bounded. By Lemma A.1 we know

[TABLE]

Also note that $Q_{h}(n)=Q(hn)$, by letting $h\to 0$, we have our claim. ∎

A.6 Continuous time random walk

Lemma A.5**.**

The $Q(s)$ in Theorem 3.3 has the following upper bound

[TABLE]

and

[TABLE]

Proof.

Let $F_{H}=F_{*}-\lambda_{G}\mathbf{I}_{N}$, then its entries are

[TABLE]

Then $\exp(F_{H}s)=e^{s\lambda_{G}}\exp(F_{*}s)$. Moreover, because $H_{*}F_{H}=F_{H}H_{*}=\mathbf{0}$, so $\exp(G_{H}s)=\exp(F_{H}s)\exp(H_{*}s)$. We first study the behavior of $\exp(F_{H}s)$.

Let $Z_{t}$ be a continuous time random walk on $\mathbb{Z}$, with rate $\lambda_{z,z\pm 1}=\lambda_{H}.$ Denote its projection on the modular space $\mathbb{Z}/N\mathbb{Z}$ as $Y_{t}$, which becomes a continuous time random walk on $\mathbb{Z}/N\mathbb{Z}$. Clearly $F_{H}$ is the transition rate matrix of $Y_{t}$. Therefore,

[TABLE]

Then by the Markov property and the translation invariance of random walk,

[TABLE]

Next we investigate the Laplace transform $e^{\beta z}$. Applying the generator of $Z_{t}$

[TABLE]

By Dynkin’s formula,

[TABLE]

therefore, $\mathbb{E}(e^{\beta Z_{t}}|Z_{0}=0)=\exp(\lambda_{F}t(e^{\beta}+e^{-\beta}-2))$. By Markov inequality,

[TABLE]

This leads to our first claim on $[\exp(F_{H}s)]_{i,j}\leq 2\exp\left(s\lambda_{G}+\lambda_{F}s(e^{\beta}+e^{-\beta}-2)-\beta\mathbf{d}(i,j)\right)$.

Next we try to compute $\exp(H_{*}s)$. Because $(H_{*}/\lambda_{H})^{n}=H_{*}/\lambda_{H}=\frac{1}{N}\mathbf{1}\cdot\mathbf{1}^{T}$ for any $n\geq 0$, so

[TABLE]

Then

[TABLE]

Its integral can be further bounded by

[TABLE]

∎

Appendix B Analytic solution of the linear model

For the linear model (4.1), collecting all the state variables $\mathbf{u}=(u_{1},u_{2},\ldots,u_{N})$, the abstract form the the system is given by

[TABLE]

According to (4.1), the coefficient matrices $\mathbf{A}$ and $\boldsymbol{\Sigma}$ are given by

[TABLE]

It is easy to verify that $\mathbf{A}$ is a negative definite matrix. Therefore, eigenvalue decomposition allows the following

[TABLE]

where $\boldsymbol{\Lambda}$ is a diagonal matrix with all diagonal entries $\lambda_{1},\ldots,\lambda_{N}$ being negative.

The solution of $\mathbf{u}$ is given by

[TABLE]

The time evolution of the covariance matrix is given by

[TABLE]

where

[TABLE]