Explicit off-line criteria for stable accurate time filtering of strongly unstable spatially extended systems

Summary of the Development of the Off-Line Stability Criteria

The goal of this study was to develop an explicit off-line test criteria for stable accurate time filtering in the above context, which is akin to the classical frozen linear constant stability test for finite difference schemes for systems of nonlinear partial differential equations (PDEs) (19). In the applications for complex spatially extended systems, the actual dynamics is typically turbulent and energetic at the smallest mesh scales but the climatological spectrum of the turbulent modes is known; for example, a mesh truncation of the compressible primitive equations with a fine mesh spacing of 10–50 km still has substantial random and chaotic energy on the smallest 10-km scales because of chaotic motion of clouds, topography, and boundary layer turbulence that is not resolved. Thus, the first step in the research program is the development of an appropriate constant coefficient stochastic PDE test problem incorporating suitable observations. In classical finite difference scheme stability analysis, the next step is to use discrete Fourier series to separate variables and reduce the analysis for an s × s PDE system to the stability properties of s × s matrix amplification factors for each spatial wave number separately (19); even though time-filtering algorithms such as the Kalman filter are nonlinear, in the second step we establish that such a reduction remains possible under natural hypotheses either exactly or with upper and lower bounds that involve the properties of an additional explicit off-line matrix inversion. With the reduction from the second step, one only needs to examine the filter stability and model error characteristic for much simpler off-line test problems depending on wave number (14, 16); for a scalar field, only decoupled complex scalar test problems need to be analyzed. This is a comparatively simple task that is illustrated here for filtering a stochastically forced heat equation.

The Constant Coefficient Test Problem

The simplest models for representing complex turbulent fluctuations involve replacing nonlinear interactions by additional linear damping and stochastic white noise forcing in time that incorporates the observed climatological spectrum and turnover time for the turbulent field (1, 2). Thus, the first step in developing analogous off-line test criteria is to use the above approximations. First, as in standard test criteria for finite difference schemes (19), the complex system of s × s PDEs is linearized at a constant coefficient background, resulting in the frozen coefficient PDE, u⃗_t = P̃(∂_x)u⃗ where at each space-time location, u⃗(x, t) belongs to s-dimensional real space. In accordance with the approximations, additional damping −γ∂_xu⃗ and white noise forcing σ(x) (t) is added to the PDE to represent the small-scale unresolved turbulent motions resulting in the basic frozen coefficient.

Reduction to Simpler Filtering Problems for Each Spatial Wave Number

For a standard finite difference approximation to Eq. 1 without any random noise with a time step Δt, the solution is expressed in standard fashion (19) in terms of the s × s amplification matrix, A_h,k, for integer k with |k| ≤ N; assume that there are p time steps so that pΔt = T and set F_h,k = (A_h,k)^p. With the finite Fourier expansion: the general discrete approximation of Eq. 1 is given at the observation times mT through its Fourier coefficients û⃗ by the block diagonal operation, In Eq. 9 the zero mean complex Gaussian noises, σ⃗_k,m, are uncorrelated in time and their second moment averages satisfy: with R_k,m^h a strictly positive definite correlation matrix and δ_j the delta function. The dynamics in the truth model (Eq. 7) also has the same block diagonal structure as in Eq. 9 with F_k = exp(P(ik)T) and complex noises satisfying Eq. 10 with related positive definite correlation matrices, R_k, which we do not compute explicitly here. To make further progress, we need additional hypotheses on the L spatial observation points and their relation to the 2N + 1 difference operator grid points, x_j = jh, j = 0, 1, …, 2N. We begin with the simplest situation.

The Observation Points Coincide with the Mesh Points

The tacit assumption here is that L = 2N + 1 and x̃_j = x_j = jh, j = 0, 1, …, 2N. Expand the discrete filter solution operator in Fourier coefficients as in Eq. 8 and under the hypothesis here compute the 2N + 1 discrete Fourier coefficients of the equations for the observations in Eq. 2. The fixed matrix multiplication by G commutes with these scalar operations so that the observations in Eq. 2 are equivalent to: with σ⃗_k,m^o the complex zero mean Gaussian noise given by σ⃗_k,m^o = (σ⃗^o, exp(ikx_j))_h. Recall from Eq. 2 that σ⃗^o has the block diagonal covariance matrix, R⃗^o = (r^oI); with this additional structure it is straightforward to calculate that the observational noise covariance matrix is given by: Thus, in the special case when the observation points in Eqs. 9 and 10 coincide with the discrete mesh points, both the discrete forward operator and the observation matrix from Eqs. 11 and 12 are block diagonal including their noise matrices for each wave number separately. In particular, we immediately have the following.

Theorem 1.

If the observation points coincide with the discrete mesh points then for both the truth model and any finite difference approximation:

1. If the covariance matrix for the initial data, u⃗₀^h, has the same block diagonal structure as in Eqs.10 and 12 for the system and observation noise, i.e. different Fourier modes are uncorrelated for k ≥ 0, then the Kalman filtering test problem is equivalent to studying the independent s × s matrix Kalman filtering problems in Eqs.9–12 for each fixed spatial wave number.

2. Provided that the s × s independent Kalman filtering problems in Eqs.9–12 are observable (see refs. 13, 14, 16 and 17), then the unique steady-state limiting Kalman filter factors for the complete model into a block product of the limiting Kalman filters for each individual wave number, k.

The practical significance of this result is that off-line tests for filter stability and model error for extremely complex PDEs can be developed for the simpler s × s matrix problems; in particular, for s = 1, this only involves independent complex scalar filtering problems. For systems with s > 1 and observation matrix G with q < s observations at each fixed point, we are in the situation with filtering with fewer observations than the actual dimension of the variables; this situation readily arises for the geophysical primitive equations where only the pressure or temperature might be known at each observation point.

The Number of Observation Points Equals the Number of Discrete Mesh Points

The situation with L = 2N + 1 so that there are 2N + 1 Fourier modes and there are 2N + 1 distinct observation points {x̃_j} _j=1 ^2N+1 with random locations drawn according to some probability distribution on [0, 2π] so that x_j ≠ x̃_j is often used in practical filtering test problems (7, 8) because observation locations are not usually at the mesh points. We show that Theorem 1 generalizes immediately to this situation with observation points distinct from the discrete mesh points with upper- and lower-bound covariance matrices for the filtering problem provided by the s × s decoupled filtering problems in Eqs. 9–12. The bounding factors involve the condition numbers of an explicit (2N + 1) × (2N + 1) matrix. For 2N + 1 distinct points, {x̃_j} _j=1 ^2N+1 on [0, 2π] with f_j ≃ f(x̃_j), 1 ≤ j ≤ 2N + 1, the map solving the scalar real-valued Fourier interpolation problem: is a unique invertible linear mapping V that maps F⃗ = (f ₁, …, f _2N+1) to (f̂_k), |k| ≤ N, i.e. V(F⃗) = f̂_k. There are even well known explicit classical numerical analysis formulas for this unequally spaced Fourier interpolation problem in a single space dimension that define V (see chapter 5 of ref. 21). The fixed matrix multiplication defining the observations in Eq. 2 obviously still commutes with G, i.e. VGF⃗ = GVF⃗, so that the discrete filtering problem for the finite difference scheme or the truth model becomes the block diagonal dynamics (Eqs. 9 and 10) together with the observation conditions for each Fourier mode,

In this situation, the observation noise covariance matrix is Gaussian mean zero and generally correlated for different Fourier modes. In the present setting with the condition û_−k = û _k ^* it is convenient to view the map V as a map from R ^2N+1 to R ^2N+1 where the image is the Fourier coefficients with k ≥ 0 and û_k = a_k + ib_k is identified with (a_k, b_k). With the assumption (Eq. 2), the noise correlation matrix for Eq. 14 is given by r^oV (I, …, I)V ^T; because V is invertible this correlation matrix satisfies the bounds: with c ₊ ² > c ₋ ² > 0 and depending on h. It is well known that the covariance matrices in Kalman filtering are monotone-increasing in the sense of positive definite matrices with respect to the noise covariance argument; furthermore, the filtering problems with the upper- and lower-bound observation covariance matrices in Eqs. 9 and 10 and Eqs. 14 and 15 obviously decouple into independent block s × s filtering problems for each wave number k. Thus, we immediately have:

Theorem 2.

If there are 2N + 1 observation points {x̃_j} _j=1 ^L, which do not coincide with the grid points, then upper and lower bounds on the discrete Kalman filtering process or the truth model as in Theorem 1 are achieved through the independent decoupled s × s filtering problems for each wave number k defined in Eqs.9, 10, and 14 with the upper- and lower-bound diagonal noise covariances in Eq.15 involving c ₋ ² , c ₊ ².

As noted earlier, the upper- and lower-bound conditions and even more improved bounds can be achieved by exploiting the explicit nature of the matrix map V when the observation points do not coincide with the mesh points; this can be done off-line for general random observation points in a certain class by exploiting the explicit interpolation formulas mentioned earlier.

Finally, we remark that there are further generalizations and variants of Theorem 1 for Kalman filtering when there are sparse observations on a regular coarse grid. This is an important test problem for filtering with sparse observations. This is developed next.

Sparse Regularly Spaced Observations.

Suppose the observations in Eq. 2 are taken at 2M + 1 grid points that are regularly spaced, i.e. x̃_j = jh̃, j = 0, 1, …, 2M with (2M + 1)h̃ = 2π. When M < N where (2N + 1)h = 2π and h denotes the mesh spacing for the difference approximation, there are fewer observations than discrete mesh points. For a given smooth function, f(x), let R_h̃f denote the restriction of f to the observation points, x̃_j and let f̂_l,h̃, |l| ≤ M, denote the discrete Fourier coefficients defined in Eqs. 6 and 7 except with respect to the observation points x̃_j with spacing h̃. For any integer k there are unique integers l, q, with |l| ≤ M, so that: With these definitions and Eq. 16, the standard aliasing formula for discrete Fourier series on the observation mesh becomes:

Given the finite difference mesh, h, with (2N + 1)h = 2π and l, with |l| ≤ M, define the observation aliasing set, 𝒜(l), for any l, |l| ≤ M where: Obviously, all of the sets 𝒜(l) are mutually disjointed. As in the proof of Theorem 1, take the discrete Fourier transform on the observation lattice of the observations in Eq. 2; now, reason as in Eqs. 11 and 12 with the facts in Eqs. 17 and 18. Then the observations for the discrete filter in Eq. 2 are equivalent to: where the observational noise, σ⃗_l,m^o, is a zero mean complex Gaussian noise with covariance: The observations in Eq. 19 couple the different wave numbers in the aliasing set, 𝒜(l), for a fixed l. With Eqs. 9 and 19 we immediately have:

Theorem 3 (Regular Sparse Observations).

With the mesh size h for the difference scheme so that (2N + 1)h = 2π and sparse observations with regular spacing h̃ with (2M + 1)h̃ = 2π, assume M < N. Under these hypotheses, the basic discrete filtering problem is equivalent to the filtering problem in Eqs.9 and 19 where the different modes in each disjoint observation aliasing set 𝒜(l) are coupled through the observation map in Eq.19. In particular, if the initial covariance matrix has a block structure respecting the aliasing sets, the filter algorithm respects the same symmetries; if the filtering system is observable, the asymptotic limiting filter decomposes into a block structure along each of the disjoint observation aliasing sets 𝒜(l). Thus, under these hypotheses, the simpler filtering problems in Eqs.9 and 19 only need to be analyzed to develop off-line test criteria.

For sparse observations with values of N = 2M, 3M, only a small number of different spatial wave numbers are coupled and simple off-line strategies as developed next can be extended readily to this situation. For very sparse observations with N ≫ M, it will be interesting to see the gain or loss in filter stability in using a fine or coarser mesh for the finite difference scheme.

Application to a Stochastically Forced Diffusion Equation

To provide a simple illustration of the development of off-line filtering criteria, we consider the perfect truth model for the scalar case in Eqs. 3 and 4 with a general energy spectrum from Eq. 5 with p̃(ik) − γ̃(ik) = −vk ² for the perfect model; thus, the truth operator is a stochastically forced heat equation. Examples of the discrete model involve using standard second-order spatial finite differences for the heat equation combined with either forward Euler, backward Euler, or symmetric Crank–Nicolson differencing in time (19). With the corresponding amplification factor for each difference method (Eq. 8) is given by: with the F_h,k in Eq. 9 given by (A_h,k)^p with pΔt = T the observation time. Note that forward Euler has strongly unstable modes for 2Δtv/h ² > 1, whereas the other two methods are strongly stable and have amplification factor smaller than one in magnitude but have model error compared with the exact symbol, F_k = exp(−Tvk ²). All of the symbols, F_k and F_h,k, are real-valued so the test problems for filter performance under the hypothesis of Theorems 1 or 2 decouple into 2N + 1 one-dimensional real filtering problems for the truth model or their finite difference approximations separately for û_k = a_k + ib_k, with individual scalar equations for each a_k and b_k. We ignore the trivial k = 0 mode in the analysis.

The Off-Line Scalar Test Problem

Thus, we only need to consider the filter performance and model error for the truth model, u_t = −γu + σ̃Ẇ, with g = 1 for the observation matrix without loss of generality. The filter performance for the perfect filtering model is defined at the observation points by: where Θ_m are independent zero mean Gaussian with variance r ≡ (σ²/2γ)[1 − exp(−2γT)]. We use u _m+1|m as a prior distribution, constrained by the observations, The observation noises are independent zero mean Gaussian with variance r^o = σ_m ². The filtering model is given by Eq. 25 for u_m+1|m+1^h and the forward dynamics: where the σ_h,m are independent zero mean Gaussian random variables with variance r_h = σ _h,m ².

With the difference schemes for the heat equation, there can be either stable |F_h| ≤ 1 or unstable |F_h| > 1 amplification factors in the model. Note that the observational noise in Eq. 25 is the same for the perfect model and the approximation but the system noise covariances for the perfect model in Eq. 24 and the approximation in Eq. 26 can be completely different. Next, we develop one optimal strategy for the system filter noise covariance r_h for a given dynamical operator, F_h, with model error different from the exact dynamics, exp(−γT), based on asymptotic filter performance.

Asymptotic Model Error and Information Criteria

Both the perfect filtering model in Eqs. 24 and 25 and the model filter in Eqs. 25 and 26 are observable and controllable, thus asymptotically stable (14, 16). Therefore, they have large time-limiting filter behavior for the covariance r_∞ and Kalman gain K_∞ as functions of the observation noise covariance, r^o, the scalar operator F = exp(−γT) or F_h, and the system covariance r or r_h. In this one-dimensional case we have explicit formulas: The perfect model limiting filter has variance and Kalman gain: while the approximate model has the limiting filter values: Next, we develop a simple objective information theoretic criterion (1, 2) to choose the approximate system noise covariance r_h in Eq. 29. First, note that the explicit function, K_∞(ỹ, z̃), satisfies: Thus, the map r_h ↦ K _∞(F_h, r^o, r_h) is a one-to-one increasing map from [0, +∞) onto its range; this range is [0, 1) for the stable case, |F_h| ≤ 1, and is [1 − F̃ ⁻², 1) for the unstable case, |F_h| > 1. Recall that the relative entropy, measures the lack of information in the probability measure p_h compared with p (1, 2). In the present application, p is the asymptotic filtering limit mean zero Gaussian measure with variance r_∞, while p_h(r_h) is the asymptotic filtering limit mean zero Gaussian measure with variance r_∞(r_h) for the fixed model parameters F_h, r^o. Clearly the best choice of the system noise r _h ^* is the one for which p has the least additional information, i.e. minimizes Eq. 31 so that:

With Eq. 30 and the discussion immediately afterward it is easy to see that r _h ^* is determined uniquely by Eq. 32: Note that this simple information theoretic criterion makes the best choice of a stable filter for the model that avoids filter divergence (14, 16), i.e., the asymptotic model filter covariance must exceed the asymptotic perfect model covariance; furthermore, for the unstable modes with |F_h| > 1, it is advantageous to let the instability operate alone and use no model system noise or greatly reduced noise under the appropriate circumstances. A rational route toward generalizations for more complex systems is clear.

Off-Line Criteria for Stable/Unstable Finite Difference Filtering and Model Error

The proposed off-line criterion for filtering the unstable/stable difference schemes in Eqs. 25 and 26 under the hypotheses of Theorems 1 and 2 is simply to use the objective criteria in Eq. 33 off-line for each different spatial wave number, k; this is readily achieved for each Fourier mode separately and an arbitrary climatological spectrum. Clearly, suboptimal choices of the model system noise, σ_h,k, rather than exactly satisfying Eq. 33 can be used in practice and assessed off-line; note that the optimal model system noise is stationary and spatially correlated when viewed in physical space.

Theorems 1 and 2 also allow for an explicit dynamic assessment of model error for the finite difference filtering model for the stochastically forced PDE. In the present simplest case of a real scalar symbol, we only need to assess the model error for decoupled scalar test problems as in Eqs. 23, 25, and 26 for each fixed wave number k and sum the results. The most interesting and simplest model error statistic for the basic scalar off-line test model is the expected value: where u_m is defined in Eqs. 23 and 25 from the perfect model. By using the augmented system for model error in chapter 6 of ref. 16, we obtain the exact recursion formulas, where K_m is the Kalman gain matrix at the mth observation time. Note that the coefficient of the first term in Eq. 35 represents the system model error; this term vanishes if there is no model error and trivially yields the fact that the Kalman filter is unbiased on the perfect model provided there is no bias at time t = 0. To gain insight into the behavior of model error in the interesting setting of an unstable mode, |F_h| > 1, replace K_m by the asymptotic Kalman gain, K_∞,h(F_h, r ^o, r) from Eq. 29 and use the lower bound, K_∞,h ≥ 1 − F _h ⁻² to obtain |y_m| < ȳ_m, where: From Eq. 23, E[u_m] = exp(−γmT)|E[u₀]| so that: assuming for simplicity that ȳ ₀ = 0. Thus, in the worst case with no model system noise in the unstable mode, there is a potentially large preconstant |exp(−γT) − F_h| |E[u₀]| but compensated by a rapidly decreasing convergent sum as m → ∞: This above worst-case analysis clearly helps to explain why stable accurate filtering can be achieved with an unstable difference approximation (17). Eqs. 35 and 36 also provide a second off-line filtering criterion beyond Eq. 33 where it can be advantageous to increase the model system noise beyond the criterion in Eq. 33 to increase the Kalman gain and thereby decrease the factor, |F_h(1 − K_∞,h)| at the expense of suboptimal filter performance.