Stable time filtering of strongly unstable spatially extended systems

Formulation

Motivated by the goal of finding a simple testable criterion for filtering stability in strongly unstable systems, we consider the generic discrete constant coefficient model filtering problem with these features for an N-dimensional vector u = (u⁻, u⁺)^⊤ with u⁻, u⁺ the subspaces of stable and unstable dynamical directions with dimensions, N⁻, N⁺, respectively, where N = N⁻ + N⁺. The linear dynamics that advances the system from the state u_m|m to u_m+1|m for m = 0, 1, …, with u_0|0 a specified Gaussian distribution is given by u_m+1|m = (u_m+1|m⁻, u_m+1|m⁺) with where the stable and unstable dynamics operators F₋, F₊ satisfy

We remark that any general linear system with only stable and unstable modes can be written in this block form through a change of variables. The forcing terms, σ_m+1⁻, σ_m+1⁺, are assumed to be Gaussian white noise vectors, independent for each m, with N⁻ × N⁻ and N⁺ × N⁺ covariance matrices given respectively by R⁻, R⁺, with R⁻ > 0, R⁺ > 0. The state u_m+1|m+1 is determined by a linear filtering strategy recursively from the prior distribution in u_m+1|m in Eq. 1 through a set of observations

The observation matrix G is M × N, so there are M observations with ῡ_m+1, the M-dimensional deterministic mean of the observations at step m+1 and σ_m+1^o, the Gaussian white noise for observations, which is assumed independent for each m with M × M covariance matrix, R^o, and R^o > 0. See ref. 14 for this formulation of Bayesian filtering.

With the assumption that u_0|0 is a Gaussian distribution with mean ū_0|0 and covariance R_0|0 > 0, the dynamics in Eq. 1 and the observations in Eq. 3 through a filtering strategy determine the state u_m+1|m, u_m+1|m+1 for m = 0, 1, 2, …, which is a Gaussian distribution at each stage with mean, ū_m+1|m, ū_m+1|m+1 and N × N covariance matrix R_m+1|m, R_m+1|m+1. Below, R^+,+, R^−,− denote the N⁺ × N⁺ and the N⁻ × N⁻ covariance matrices of the stable and unstable modes, u⁺ and u⁻, respectively, with R^+,− = (R^−,+)^⊤ the cross-covariance matrix between u⁻ and u⁺. A filtering strategy is unbiased if

The most important unbiased filtering strategy is the Kalman filter (12⇓–14), which is the unbiased linear filtering strategy that minimizes the covariance matrix, R_m+1|m+1, among all unbiased linear filtering strategies. This property is extensively discussed in the above references, and a systematic derivation of the explicit recursion formulas for the mean ū_m|m and covariance R_m|m of the Kalman filter from the Bayesian perspective will be used in the next section to construct two stable unbiased filters.

Here, the main interest is the stability of such a filtering process despite the presence of the strongly unstable modes u⁺ in the dynamics. Clearly, such stability can occur only if the observations interact sufficiently with the unstable dynamics in a suitable fashion to tame the growth of these modes. Assume that the mean values of the observations, ῡ_m, are uniformly bounded for all m, i.e., ‖ῡ_m‖ ≤ C₀, then a filtering strategy is stable provided that the Gaussian distribution, u_m|m, satisfy where the constant C depends on C₀. Clearly, the bound in Eqs. 5 and 1 guarantees a similar bound for ū_m+1|m, R_m+1|m. In Eq. 5, ‖ ‖ denotes any norm on the appropriate finite dimensional space.

Let P₊, P₋ denote the projections on the unstable and stable directions, respectively, so that P₊u = u⁺, P₋u = u⁻. We claim the following.

Two Stable Unbiased Estimators

As mentioned in the previous section, here unbiased stable filters are given using the Bayesian approach to filtering (14). Here and below, p(x) denotes the probability density of the random (vector) variable x, and p(x|y) is the probability of the random variable x conditioned on the variable y. Bayes’ theorem (12, 14) guarantees that In the Bayesian approach to filtering (14), the objective is to estimate the probability density where Bayes’ theorem from Eq. 9 has been used twice in Eq. 10.

The probability distribution on the left-hand side of Eq. 10 is the desired estimated probability distribution including the observations at m + 1, whereas p(u_m+1|m) is the prior distribution including all dynamics and observations before the information from the observations at m + 1 is taken into account; the probability density, p(υ_m+1|u_m+1|m+1) is evaluated explicitly by using Eq. 3. In the present setting, all probability distributions are Gaussian and are readily calculated explicitly (12, 14).

With this preliminary background, we build the two unbiased estimators. For the initial probability density, p(u₀), we use Bayes’ theorem and write

Prototype Model: Stable Filtering for a Strongly Unstable Finite Difference Scheme

First, an elementary model for the randomly fluctuating physical process is constructed. Consider a random variable u(x, t), which is 2π periodic in x and has the Fourier expansion with u_−k = (u_k)*, u₀ = 0. Here each Fourier mode for k > 0 satisfies the stochastic linear differential equation where the W_k are independent complex Wiener processes for each k, and the independent real and imaginary parts have variance σ_k. The coefficients c, μ in Eq. 19 satisfy c > 0, μ > 0 so that in physical space Eqs. 18 and 19 define the solution of a stable convection-diffusion equation with wave speed c > 0 and viscosity μ > 0 randomly forced by spatially correlated stationary white noise in time. The statistical equilibrium distribution for the kth Fourier mode in Eq. 19 is a Gaussian with zero mean and variance σ_k/ k. Here the σ_k are chosen according to the power law with σ₀ a fixed constant, either 0.5 or 5 below, and α = 2. Below, the physical signal that is filtered arises from a realization of this random field with prescribed initial data corresponding to a Gaussian pulse. The parameters c and μ are fixed to the values c = 1 and μ = 0.1 below.

The discrete model solution u^h(x, t) is defined at the equidistant grid points. and satisfies periodicity, u^h(x₀) = u^h(x_N). Here, N is even so that π/h = N/2 and the N (discrete) Fourier coefficients of u^h are given through the discrete Fourier and inverse Fourier transforms Because u^h is real, so are u₀^h and u_N/2^h because exp(i(N/2)x_j) = (−1)^j is real for any j with N even. Below, N = 42 is used with u₀^h, u_N/2^h set to zero so there are 40 active discrete modes. This special set-up is used for simplicity in rapid evaluation of the unstable rank condition in Theorem 1.

To define the approximate dynamics, the spatial derivatives in the forced convection–diffusion equation defined by Eqs. 18–20 are replaced by standard central finite differences in physical space. Their action in Fourier space (the symbol) can be computed easily by applying directly any finite difference operator to Eq. 21 yielding The time derivative is replaced by a simple forward Euler step while the form of the random forcing in each discrete Fourier component coincides with that in the physical model in Eq. 19 integrated over one time step.

The time step Δt = 1/8 is used throughout the study below. The corresponding discrete model is an unstable finite difference approximation to the stable convection-diffusion equation with random forcing. The amplification factor of the difference operator (15) at each discrete Fourier mode is plotted in Fig. 1. There are 40 active modes and 12 unstable modes corresponding to the 6 complex Fourier modes with wave numbers 15 ≤ k ≤ 20.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 1.

The amplification factor per Euler step with Δt = 1/8 of the 40 Fourier modes. For the typical observation time T_obs = 1, the effective amplification factors between observations are the 8th power of those shown here; hence, those outside the unit circle are strongly unstable.

In the experiments reported below, L observations are taken at fixed locations x_l, 1 ≤ l ≤ L, picked at random from the uniform distribution in [0, 2π] with Eq. 21 interpolating the solution in terms of the discrete Fourier coefficients. Below, the observation values are taken from the physical solution discussed in the first paragraph of this section with independent observation noise from site to site with fixed variance σ₀ = 0.05. The observation time T_obs is an integer multiple of Δt = 1/8 in all experiments below. The standard value T_obs = 1 is used in most of the numerical experiments; this is a severe test problem for stable filtering in a strongly unstable spatially extended system, because the growth rates of instability for 15 ≤ k ≤ 20 are the eighth power of those in Fig. 1.

In Table 1, we report the outcome of a variety of filtering experiments with the unstable system as the number of observations, L, and the observation time, T_obs, were varied for the time interval 0 ≤ t ≤ 100. The data reported include the largest and smallest eigenvalues of the asymptotic covariance matrix, Γ_min and Γ_max, achieved by time t = 20 in all cases, as well as their ratio, the condition number; the amplitude of the filtered solution at time t = 100 is also recorded to compare its magnitude with the physical solution at that time. The first p⁺ that satisfies the full rank test on the unstable modes also is reported there as a theoretical guideline. Although there are simple examples in the present context illustrating Theorem 1, here the focus is on the central issue of stable accurate filtering of strongly unstable systems.

The results of the Kalman-filtering algorithm with 10 random observations in [0, 2π] for T_obs = 1 using the strongly unstable difference method are reported in Fig. 2 and the first row of Table 1 for the noise parameter σ₀ = 5. Notable is the accurate filtering in Fig. 2 at times t = 50, 100, with the unstable difference scheme even though there are only 10 observations and there are 12 strongly unstable modes; in this case, the finite rank test for stability in Theorem 1 is satisfied with p⁺ = 1, and the condition number of the covariance matrix is of order 10⁴ with the largest eigenvalue of order 1. Fig. 3 and the second row of Table 1 depict the results of the Kalman-filtering algorithm with the plentiful number of 40 random observations in [0, 2π] for T_obs = 1 and the less-noisy physical signal with σ₀ = 0.5. Here the number of observations equals the total number of modes, and the largest eigenvalues of the covariance matrix are ≈1.5 × 10⁻⁴ with condition number ≈35. It is no surprise that p⁺ = 0 for the rank test in Theorem 1, and the physical solution is nearly reproduced exactly in the unstable discrete model at times t = 50, 100.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

Comparison of physical and filtered solutions. (Top) The minimal and maximal eigenvalues of the Kalman-filtered covariance matrix, R_m|m. (Middle and Bottom) The solution (solid lines) u(x, t) of Eqs. 18–20 and the filtered solution (dashed lines), u_m|m, are shown at t = 50 (Middle) and t = 100 (Bottom) with 10 observations, σ₀ = 0.5, and T_obs = 1 (see also the first row of Table 1).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 3.

Comparison of physical and filtered solutions. (Top) The minimal and maximal eigenvalues of the Kalman filtered covariance matrix, R_m|m (Top). (Middle and Bottom) The solution (solid lines) u(x, t) of Eqs. 18–20 and the filtered solution (dashed lines), u_m|m, are shown at t = 50 (Middle) and t = 100 (Bottom) with 40 observations, σ₀ = 0.5, and T_obs = 1 (see also the second row of Table 1).

Practical Failure of Standard Observability Tests in Strongly Unstable Systems

Cohn and Dee (17) interpreted the standard observability criterion for stable filtering for finite difference approximations in interesting work. In particular, they proved that if the amplification factor is distinct for different wave numbers, then a single generic observation is sufficient for stability of the discrete filter process. Although they did not apply this criterion to unstable difference schemes, there is no theoretical limitation in doing so; in fact, the result in ref. 17 is true for the strongly unstable difference scheme in the present work because the amplification factors are distinct as shown in Fig. 1. The issue arises of whether this is a practical filtering stability criterion for the strongly unstable finite difference scheme.

In Fig. 4 and the third row of Table 1, the results of a Kalman-filtering experiment with a single random observation satisfying the criterion from ref. 17 are reported with T_obs = 1 and σ₀ = 0.5. As predicted by ref. 17, the filtering process remains stable as can be seen from the level value of the largest eigenvalue of the covariance matrix in Fig. 4 Top as well as the amplitude of the filtered solution, which actually decreases at t = 100 compared with t = 50. Furthermore, p⁺ = 11 for the unstable rank criterion. However, the condition number of the covariance matrix is of order 10¹³, and the physical space filtered amplitudes are of order 10⁴ compared with the actual order 1 physical amplitude reported in Fig. 3. Thus, practical stable filtering has failed spectacularly in this example because the constant in Eq. 7 is extremely small, whereas p⁺ = 11 is very large.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 4.

Comparison of physical and filtered solutions. (Top) The minimal and maximal eigenvalues of the Kalman filtered covariance matrix, R_m|m (Top); (Middle and Bottom) The solution (solid lines) u(x, t) of Eqs. 18–20 and the filtered solution (dashed lines), u_m|m, are shown at t = 50 (Middle) and t = 100 (Bottom) with one observation, σ₀ = 0.5, and T_obs = 1 (see also the third row of Table 1).

In the experiment reported in Fig. 5 and the fourth row of Table 1, the same single random observation was used, but the frequency of observations was increased so that T_obs = 1/8, the discrete time step. Because the amplification of instability is decreased substantially between the more frequent observations and the stability criterion is satisfied with p⁺ = 11, one might anticipate better filtering properties than in the previous situation. Indeed this improvement occurs, but the condition number of the filter is O(10⁶) whereas the amplitude of the filtered estimator is ≈25 times the amplitude of the true filtered solution without any accuracy, so there is the role of the frequency and number of observations in filter stability.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 5.

Comparison of physical and filtered solutions. (Top) The minimal and maximal eigenvalues of the Kalman filtered covariance matrix, R_m|m (Top). (Middle and Bottom) The solution (solid lines) u(x, t) of Eqs. 18–20 and the filtered solution (dashed lines), u_m|m, are shown at t = 50 (Middle) and t = 100 (Bottom) with one observation, σ₀ = 0.5, and T_obs = 1/8 (see also the fourth row of Table 1).

The remaining cases listed in Table 1 indicate the variability of the stability parameters. In all cases depicted in Table 1, the largest and smallest eigenvalues of the filtered covariance matrix approach their asymptotic values by t = 20. The results strongly confirm the trends already discussed with the successful filtering cases in Figs. 2 and 3 and the stable but failed filtering cases in Figs. 4 and 5. Both increasing the number of observations and decreasing the observation time improves the stability properties of the filter.

Concluding Discussion

The above results strongly suggest the possibility that there is a practical off-line stability criterion for practical stable and accurate filtering with strongly unstable finite difference operators using low values of p⁺ in the rank test and sufficiently plentiful observations in a quantitative fashion. Lack of space prevents a detailed discussion here.