New Research In
Physical Sciences
Social Sciences
Featured Portals
Articles by Topic
Biological Sciences
Featured Portals
Articles by Topic
 Agricultural Sciences
 Anthropology
 Applied Biological Sciences
 Biochemistry
 Biophysics and Computational Biology
 Cell Biology
 Developmental Biology
 Ecology
 Environmental Sciences
 Evolution
 Genetics
 Immunology and Inflammation
 Medical Sciences
 Microbiology
 Neuroscience
 Pharmacology
 Physiology
 Plant Biology
 Population Biology
 Psychological and Cognitive Sciences
 Sustainability Science
 Systems Biology
Stable time filtering of strongly unstable spatially extended systems

Contributed by Andrew J. Majda, March 27, 2006
Abstract
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 scaleup of smallscale 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 lowdimensional largescale 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.
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 Ndimensional 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_{mm} to u_{m+1m} for m = 0, 1, …, with u_{00} a specified Gaussian distribution is given by u_{m+1m} = (u_{m+1m}^{−}, u_{m+1m}^{+}) 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}_{+1m}_{+1} is determined by a linear filtering strategy recursively from the prior distribution in u_{m}_{+1m} in Eq. 1 through a set of observations
The observation matrix G is M × N, so there are M observations with ῡ_{m}_{+1}, the Mdimensional 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_{00} is a Gaussian distribution with mean ū_{00} and covariance R_{00} > 0, the dynamics in Eq. 1 and the observations in Eq. 3 through a filtering strategy determine the state u_{m}_{+1m}, u_{m}_{+1m}_{+1} for m = 0, 1, 2, …, which is a Gaussian distribution at each stage with mean, ū_{m}_{+1m}, ū_{m}_{+1m}_{+1} and N × N covariance matrix R_{m}_{+1m}, R_{m}_{+1m}_{+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 crosscovariance 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}_{+1m}_{+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_{0}, then a filtering strategy is stable provided that the Gaussian distribution, u_{m}_{m}, satisfy where the constant C depends on C_{0}. Clearly, the bound in Eqs. 5 and 1 guarantees a similar bound for ū_{m}_{+1m}, R_{m}_{+1m}. 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.
Theorem 1.
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, u_{p}^{+}, with ‖u_{p}^{+}‖ = 1 and Assume the initial mean, ū_{00} = (0, u^{+}) and the mean of the observations vanishes, i.e., ῡ_{m} = 0 for all m. Then for 0 ≤ m ≤ p, ū_{m}_{m} = (F_{+})^{m}u_{p}^{+} 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, R_{m}_{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.
Remark 1:
It is worthwhile to note here that the matrix problem in Eqs. 6 and 7 can still be illconditioned 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.
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(xy) 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 lefthand side of Eq. 10 is the desired estimated probability distribution including the observations at m + 1, whereas p(u_{m}_{+1m}) 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}_{+1m}_{+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_{0}), we use Bayes’ theorem and write
First Unbiased Estimator.
We use the first formula in Eq. 11 and generate an estimator, u_{m+1m+1}^{−} = u_{m+1m}^{−}, 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 u_{m+1m+1}^{−} satisfies the required stability bounds in Eq. 5. Next, let u_{m+1m+1}^{+} be the Kalmanfiltered 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 p_{0}(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 Kalmanfiltering problem in Eq. 5 since ū_{mm}^{−}, R_{mm−,}^{−} are uniformly bounded and R^{+} > 0. Furthermore, with Eq. 12 it is easy to check that the estimator (u_{mm}^{−}, u_{mm}^{+}) determined in the above fashion is unbiased and satisfies Eq. 4 by construction. This completes the proof of Theorem 1.
Remark 2:
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 Kalmanfiltering 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.
Second Unbiased Estimator.
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 lowerorder 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_{+}^{o}G, 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 u_{m+1m}^{+}, u_{m+1m+1}^{+} recursively from the Kalmanfiltering 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 u_{mm}^{+}. Next, given u_{mm}^{+} determined in the above fashion, let u_{mm}^{−} 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 u_{m+1m+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.
Remark 3:
These two filtering algorithms might prove useful beyond their theoretical role when combined iteratively in a fractional step algorithm.
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} = 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 convectiondiffusion 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 u^{h}(x, t) is defined at the equidistant grid points. and satisfies periodicity, u^{h}(x_{0}) = 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_{0}^{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_{0}^{h}, u_{N/2}^{h} set to zero so there are 40 active discrete modes. This special setup 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 convectiondiffusion 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.
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} = 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 Kalmanfiltering 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 σ_{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 10^{4} with the largest eigenvalue of order 1. Fig. 3 and the second row of Table 1 depict the results of the Kalmanfiltering algorithm with the plentiful number of 40 random observations in [0, 2π] for T_{obs} = 1 and the lessnoisy 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.
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 Kalmanfiltering experiment with a single random observation satisfying the criterion from ref. 17 are reported with T_{obs} = 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 10^{13}, and the physical space filtered amplitudes are of order 10^{4} 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.
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^{6}) 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.
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 offline 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.
Acknowledgments
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.
Footnotes
 ↵^{§}To whom correspondence should be addressed. Email: jonjon{at}cims.nyu.edu

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.
 © 2006 by The National Academy of Sciences of the USA
References
 ↵
 Majda A.,
 Wang X.
 ↵
 Majda A.,
 Abramov R.,
 Grote M. J.
 ↵
 Berliner L. M.,
 Milliff R. E.,
 Wikle C. K.
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 Farell B.,
 Ioannou P.
 ↵
 Farell B.,
 Ioannou P.
 ↵
 Ott E.,
 Hunt B.,
 Szunyogh I.,
 Zimin A.,
 Kostelich E.,
 Corrazza M.,
 Kalnay E.,
 Yorke J.
 ↵
 Chui C. K.,
 Chen G.
 ↵
 Jazwinski A.
 ↵
 Kiapio J.,
 Somersalo E.
 ↵
 Richtmeyer R.,
 Morton K. W.
 ↵
 Chorin A.,
 Krause P.
 ↵