Contributed by Andrew J. Majda, March 27, 2006
Many contemporary problems in science involve making predictions based on partial observation of extremely complicated spatially extended systems with many degrees of freedom and with physical instabilities on both large and small scale. Various new ensemble filtering strategies have been developed recently for these applications, and new mathematical issues arise. Because ensembles are extremely expensive to generate, one such issue is whether it is possible under appropriate circumstances to take long time steps in an explicit difference scheme and violate the classical Courant–Friedrichs–Lewy (CFL)-stability condition yet obtain stable accurate filtering by using the observations. These issues are explored here both through elementary mathematical theory, which provides simple guidelines, and the detailed study of a prototype model. The prototype model involves an unstable finite difference scheme for a convection–diffusion equation, and it is demonstrated below that appropriate observations can result in stable accurate filtering of this strongly unstable spatially extended system.
Many contemporary problems in science ranging from protein folding in molecular dynamics to scale-up of small-scale effects in nanotechnology to making accurate predictions of the coupled atmosphere–ocean system involve partial observations of extremely complicated systems with many degrees of freedom. Novel mathematical issues arise in the attempt to quantify the behavior of such complex multiscale systems (1, 2). For example, in the coupled atmosphere–ocean system, the current practical models for prediction of both weather and climate involve general circulation models where the physical equations for these extremely complex flows are discretized in space and time, and the effects of unresolved processes are parametrized according to various recipes; the result of this process involves a model for the prediction of weather and climate from partial observations of an extremely unstable, chaotic dynamical system with several billion degrees of freedom. The nature of the physical instabilities that strongly affect the predictive properties of this system range from (i) comparatively low-dimensional large-scale instabilities involving synoptic scale weather activity to (ii) inherently statistical instabilities on shorter spatiotemporal scales such as those involving moist convection, which crucially affects the water vapor in the atmosphere and rainfall.
Bayesian hierarchical modeling (3) and reduced order filtering strategies (4⇓⇓⇓⇓⇓⇓–11) have been developed with some success in these extremely complex systems including the role of observations in tracking the instabilities of type i. The basis for such dynamic prediction strategies for the complex spatially extended systems is the classical Kalman filtering algorithm (12⇓–14). New issues arise in the practical application of these filtering strategies to complex spatially extended systems, and this topic is the focus for the present work.
One new mathematical issue that emerges is the following one: ensemble filtering requires multiple realizations of an extremely expensive dynamical system with many degrees of freedom; with these practical limitations, it is extremely interesting to see whether it is possible to use large time steps, which violate the classical Courant–Friedrichs–Lewy (CFL)-stability condition, for an explicit difference scheme (15) and, nevertheless, obtain stable and statistically accurate filtering. Such counterintuitive behavior has emerged in recent practical applications (3, 7, 8) without documentation or mathematical understanding of this potentially practically important phenomenon. The remainder of this work is devoted to elucidating this phenomenon as well as the emerging reduced order filtering strategies for tracking physical instabilities of types i and ii mentioned earlier. First, some elementary mathematical theorems are developed to supply a context and justification for the possibility for various types of stable filtering strategies for strongly unstable systems. Then, these mathematical results are used as guidelines for a detailed investigation of the key issue raised at the beginning of this paragraph for a prototype model involving an unstable finite difference approximation for a convection–diffusion equation.
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 um|m to um+1|m for m = 0, 1, …, with u0|0 a specified Gaussian distribution is given by um+1|m = (um+1|m−, um+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 um+1|m+1 is determined by a linear filtering strategy recursively from the prior distribution in um+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+1o, the Gaussian white noise for observations, which is assumed independent for each m with M × M covariance matrix, Ro, and Ro > 0. See ref. 14 for this formulation of Bayesian filtering.
With the assumption that u0|0 is a Gaussian distribution with mean ū0|0 and covariance R0|0 > 0, the dynamics in Eq. 1 and the observations in Eq. 3 through a filtering strategy determine the state um+1|m, um+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 Rm+1|m, Rm+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, Rm+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 Rm|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‖ ≤ C0, then a filtering strategy is stable provided that the Gaussian distribution, um|m, satisfy
where the constant C depends on C0. Clearly, the bound in Eqs. 5 and 1 guarantees a similar bound for ūm+1|m, Rm+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.
A necessary and sufficient condition for stable filtering in the sense of Eq. 5 for the unstable system in Eq. 1 is that there exists some integer p ≥ 0 so that with P+u+ = u+ the matrix
has full rank, i.e., there exists c > 0 so that
The algebraic test condition in Eq. 6 generalizes the classical observability condition in filtering theory (12, 13), as a requirement only on the unstable modes.
It is easy to show that if Eq. 6 is violated, then the stability requirement in Eq. 5 cannot be satisfied. If Eq. 6 is violated, then for any p, there exists a vector, up+, with ‖up+‖ = 1 and
Assume the initial mean, ū0|0 = (0, u+) and the mean of the observations vanishes, i.e., ῡm = 0 for all m. Then for 0 ≤ m ≤ p, ūm|m = (F+)mup+ and since Eq. 2 easily guarantees ‖F+u+‖ ≥ (1 + δ)‖u+‖ with some fixed δ > 0 and a suitable norm ‖ ‖, this sequence violates the bound in Eq. 5 for the mean.
The proof of the stability condition in Eq. 5 when Eq. 6 is satisfied is a consequence of the explicit estimators constructed in the next section; there, we construct unbiased estimators that are stable and satisfy Eq. 5. Because the Kalman filter has the smallest covariance matrix among all unbiased filters, this approach guarantees the uniform bound on the covariance matrix, Rm|m, for the Kalman filter. The uniform bound on the mean ūm|m in Eq. 5 is proved below through the unbiased estimators presented next. This completes the proof of Theorem 1.
It is worthwhile to note here that the matrix problem in Eqs. 6 and 7 can still be ill-conditioned so that the theoretical constant c in Eq. 7 is very small. Thus, care is needed in interpreting Theorem 1. This issue will be discussed in the applications below where such phenomena actually occur.
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(um+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|um+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(u0), we use Bayes’ theorem and write
We use the first formula in Eq. 11 and generate an estimator, um+1|m+1− = um+1|m−, in trivial fashion by solving the first stable Eqs. 1 with initial data p(u−) and completely ignoring the observations for F−. Because F− has all eigenvalues smaller than one in modulus, it is easy to see that um+1|m+1− satisfies the required stability bounds in Eq. 5. Next, let um+1|m+1+ be the Kalman-filtered solution for the second unstable dynamical equation in 1 with the observations
From the hypothesis in Eq. 1, the unstable dynamical equations for u+ in Eqs. 1 and 12 combine to form an observable Kalman filtering problem for the variable u+ alone with the initial density p0(u+|u−). This formula operates on the reduced subspace of u+. Lemma 6.1 from ref. 12 and theorem 7.4 of ref. 13 guarantee the stability of this reduced Kalman-filtering problem in Eq. 5 since ūm|m−, Rm|m−,− are uniformly bounded and R+ > 0. Furthermore, with Eq. 12 it is easy to check that the estimator (um|m−, um|m+) determined in the above fashion is unbiased and satisfies Eq. 4 by construction. This completes the proof of Theorem 1.
Although the above unbiased filter is constructed in a straightforward fashion, it establishes that there is a stable Kalman filter operating on the reduced subspace of N+ unstable modes provided that these unstable modes are determined from the observations through the precise observability condition in Eq. 6. In the present context, this filter provides rigorous justification for the reduced order Kalman-filtering strategies with N+ ≪ N developed by several authors (6, 9, 10, 16) and also provides a precise quantitative requirement for stability of such algorithms through the requirement of observability in Eq. 6 or 7.
The first unbiased estimator does not use any of the observations to estimate the mean and reduce the variance in the stable directions. It provides useful theoretical support for a lower-order reduced filtering strategy for the N+ unstable modes alone when these are the quantities of fundamental physical interest. Under appropriate hypotheses (6, 9, 10, 16), it is possible to filter observations in both u+ and u− in a stable fashion taking all of the observations into account. This setting is an appropriate one for the context of unstable difference schemes where one wants accurate filtering of the stable modes and stability through the observations for the unstable modes (3, 7). To develop these hypotheses, assume that the range of the observations restricted to the stable subspace is not all of ℝM, i.e., GP− has a nontrivial orthogonal complement, or, equivalently, assume
Let P+o denote the orthogonal projection on the subspace, Ker(P−, G⊤); this condition is automatically satisfied if M > N− and by construction
Now, assume that the reduced family of observations, P+oG, are observable for the unstable dynamical system in the second equation in 9; i.e., there exists an integer p ≥ 0 and c > 0 so that
Under these hypotheses, reverse the roles of u+ and u− and use the second formula in 11 for the initial probability density. Determine um+1|m+, um+1|m+1+ recursively from the Kalman-filtering algorithm for the dynamical system, with observations
According to lemma 6.1 in ref. 12 and theorem 7.4 from ref. 13 with Eq. 15 and R+ > 0, this filtering algorithm is stable and satisfies the bounds in Eq. 6 for um|m+. Next, given um|m+ determined in the above fashion, let um|m− be determined by the Kalman filter applied to the stable dynamical system in Eq. 1 with the observations
with initial probability distribution, p(u−|u+). Because the dynamical system in 1 for u− without observations is stable, and the observations only decrease variance (see lemma 6.2 in ref. 12), the filtering algorithm in 17 is stable, and um+1|m+1− satisfies Eq. 6. Adding Eq. 16 with Eq. 17 and using Eq. 15 shows that this algorithm generates an unbiased filter, which satisfies Eq. 4. All these facts combine to show that the present algorithm is a stable unbiased estimator under the stronger hypothesis of observability of unstable modes in Eq. 15.
These two filtering algorithms might prove useful beyond their theoretical role when combined iteratively in a fractional step algorithm.
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 = (uk)*, u0 = 0. Here each Fourier mode for k > 0 satisfies the stochastic linear differential equation
where the Wk 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 σ0 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 uh(x, t) is defined at the equidistant grid points.
and satisfies periodicity, uh(x0) = uh(xN). Here, N is even so that π/h = N/2 and the N (discrete) Fourier coefficients of uh are given through the discrete Fourier and inverse Fourier transforms
Because uh is real, so are u0h and uN/2h because exp(i(N/2)xj) = (−1)j is real for any j with N even. Below, N = 42 is used with u0h, uN/2h 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.
The amplification factor per Euler step with Δt = 1/8 of the 40 Fourier modes. For the typical observation time Tobs = 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 xl, 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 = 0.05. The observation time Tobs is an integer multiple of Δt = 1/8 in all experiments below. The standard value Tobs = 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, Tobs, 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 maximum of the Kalman-filtered solution (in physical space), max |um|m|, and the condition number, Γmax/Γmin, of the asymptotic covariance matrix, Rm|m, at t = 100 for L observations
The results of the Kalman-filtering algorithm with 10 random observations in [0, 2π] for Tobs = 1 using the strongly unstable difference method are reported in Fig. 2 and the first row of Table 1 for the noise parameter σ0 = 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 104 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 Tobs = 1 and the less-noisy physical signal with σ0 = 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−4 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.
Comparison of physical and filtered solutions. (Top) The minimal and maximal eigenvalues of the Kalman-filtered covariance matrix, Rm|m. (Middle and Bottom) The solution (solid lines) u(x, t) of Eqs. 18–20 and the filtered solution (dashed lines), um|m, are shown at t = 50 (Middle) and t = 100 (Bottom) with 10 observations, σ0 = 0.5, and Tobs = 1 (see also the first row of Table 1).
Comparison of physical and filtered solutions. (Top) The minimal and maximal eigenvalues of the Kalman filtered covariance matrix, Rm|m (Top). (Middle and Bottom) The solution (solid lines) u(x, t) of Eqs. 18–20 and the filtered solution (dashed lines), um|m, are shown at t = 50 (Middle) and t = 100 (Bottom) with 40 observations, σ0 = 0.5, and Tobs = 1 (see also the second row of Table 1).
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 Tobs = 1 and σ0 = 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 1013, and the physical space filtered amplitudes are of order 104 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.
Comparison of physical and filtered solutions. (Top) The minimal and maximal eigenvalues of the Kalman filtered covariance matrix, Rm|m (Top); (Middle and Bottom) The solution (solid lines) u(x, t) of Eqs. 18–20 and the filtered solution (dashed lines), um|m, are shown at t = 50 (Middle) and t = 100 (Bottom) with one observation, σ0 = 0.5, and Tobs = 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 Tobs = 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(106) 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.
Comparison of physical and filtered solutions. (Top) The minimal and maximal eigenvalues of the Kalman filtered covariance matrix, Rm|m (Top). (Middle and Bottom) The solution (solid lines) u(x, t) of Eqs. 18–20 and the filtered solution (dashed lines), um|m, are shown at t = 50 (Middle) and t = 100 (Bottom) with one observation, σ0 = 0.5, and Tobs = 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.
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.
The topic of this work was inspired by a lecture of Chris Wikle on ref. 3 at an interdisciplinary workshop at the National Center for Atmospheric Research (NCAR), sponsored by the Statistical and Applied Mathematical Sciences Institute, where A.J.M. noted publicly that successful Bayesian filtering was done with a strongly unstable numerical method; Jeff Anderson of NCAR then communicated to us similar examples in his own research. We warmly acknowledge these discussions. A.J.M. was supported in part by grants from the National Science Foundation and the Office of Naval Research.
Author contributions: M.J.G. and A.J.M. designed research; M.J.G. and A.J.M. performed research; and A.J.M. wrote the paper.
Conflict of interest statement: No conflicts declared.
