Sampling, feasibility, and priors in data assimilation

. Importance sampling algorithms are discussed in detail, with an emphasis on implicit sampling, and applied to data assimilation via particle ﬁl- ters. Implicit sampling makes it possible to use the data to ﬁnd high-probability samples at relatively low cost, making the assimilation more eﬃcient. A new analysis of the feasibility of data assimilation is presented, showing in detail why feasibility depends on the Frobenius norm of the covariance matrix of the noise and not on the number of variables. A discussion of the convergence of particular particle ﬁlters follows. A major open problem in numerical data assimilation is the determination of appropriate priors; a progress report on recent work on this problem is given. The analysis highlights the need for a careful attention both to the data and to the physics in data assimilation problems.


1.
Introduction. Bayesian methods for estimating parameters and states in complex systems are widely used in science and engineering; they combine a prior distribution of the quantities of interest, often generated by computation, with data from observations, to produce a posterior distribution from which reliable inferences can be made. Some recent applications of these methods, for example in geophysics, involve more unknowns, larger data sets, and problems that are more nonlinear than earlier applications, and present new challenges in their use (see, e.g. [7,50,21,48,16] for recent reviews of Bayesian methods in engineering and physics).
The emphasis in the present paper is on data assimilation via particle filters, which requires effective sampling. We give a preliminary discussion of importance sampling and explain why implicit sampling [36,3,14,13] is an effective importance sampling method. We then present implicit sampling methods for calculating a posterior distribution of the unknowns of interest, given a prior distribution and a distribution of the observation errors, first in a parameter estimation problem, then in a data assimilation problem where the prior is generated by solving stochastic differential equations with a given noise. Conditions for data assimilation to be feasible in principle and their physical interpretation are discussed (see also [12]). A linear convergence analysis for data assimilation methods shows that Monte Carlo methods converge for many physically meaningful data assimilation problems, provided that the numerical analysis is appropriate and that the size of the noise is small enough in the appropriate norm, even when the number of variables is very large. The keys to success are a careful use of data and a careful attention to correlations.
An important open problem in the practical application of data assimilation methods is the difficulty in finding suitable error models, or priors. We discuss a family of problems where priors can be found, with an example. Here too a proper use of data and a careful attention to correlations are the keys to success.
2. Importance sampling. Consider the problem of sampling a given probability density function (pdf) f on the computer, i.e., given a pdf for an m-dimensional random vector, construct a sequence of vectors X 1 , X 2 , . . . , which can pass statistical independence tests, and whose histogram converges to the graph of f , in symbols, find X i ∼ f. We discuss this problem in the context of numerical quadrature. Suppose one wishes to evaluate the integral I = g(y)f (y)dy, where g(y) is a given function and f (y) is a pdf, i.e., f (y) ≥ 0 and f dy = 1. The integral I equals E[g(x)], where x is a random variable with pdf f and E[·] denotes an expected value. This integral can be approximated via the law of large numbers as where X i ∼ f and M is a large integer (the number of samples we draw). The error is of the order of σ(g(x))/M 1/2 , where σ(·) denotes the standard deviation of the variable in the parentheses. This assumes that one knows how to find the samples X i , which in general is difficult. One way to proceed is importance sampling (see, e.g., [26,9]): introduce an auxiliary pdf f 0 , often called an importance function, which never vanishes unless f does, and which one knows how to sample, and rewrite the integral as where the pdf of x * is f 0 and w(x * ) = f (x * )/f 0 (x * ). The expected value can then be approximated by where X * i ∼ f 0 and the w(X * i ) = f (X * i )/f 0 (X * i ) are the sampling weights. The error in this approximation is of order σ(gw)/M 1/2 (here the standard deviation is computed with respect to the pdf f 0 ). The sampling weights make up for the fact that one is sampling the wrong pdf.
Importance sampling can work well if f and f 0 are fairly close to each other. Suppose they are not (see figure 1). Then many of the samples X i (drawn from the importance function) fall where f is negligible and their sampling weight is small. The computing effort used to generate such samples is wasted since the low-weight samples contribute little to the approximation of the expected value. In fact, one can define an effective number of samples [16,2,53,28] as where expected values are computed with respect to the importance function f 0 . In particular, M eff = M if all the particles have weights 1/M . The effective number of samples approximates the equivalent number of independent, unweighted samples one has after importance sampling with M samples. Importance sampling is efficient if R is much smaller than M . For example, suppose that you find that one weight is much larger than all the others. Then R ≈ M , so that one has effectively one sample and this sample does not necessarily have a high probability. In this case, the importance sampling scheme has "collapsed". To find f 0 that prevents a collapse, one has to know the shape of f quite well, in particular know the region where f is large. In some of the examples below, the whole problem is to estimate where a particular pdf is large, and it is not realistic to assume that this is known in advance.
As the number of variables increases, the value of a carefully designed importance function also increases. This can be illustrated by an example. Suppose that f = N (0, I m ), where I m is the identity matrix of order m (here and below we denote a Gaussian with mean µ and covariance matrix Q by N (µ, Q)). To sample this Gaussian we pick a Gaussian importance function importance function has a shape similar to that of the function we wish to sample and f 0 is large where f is large (for moderate σ). Nonetheless, sampling becomes increasingly costly as the number of components of x increases. We find that where the superscript T denotes a transpose, so that Here, the expected values are computed with respect to the importance function. Equation (2) shows that the number of samples required for a given M eff grows exponentially with the number of variables m: The situation is illustrated in figure 2. As the number of variables increases, the cost of sampling, measured by the number of samples required to obtain the pre-specified number of effective samples, grows exponentially with the number of variables for any σ > 1/2 except σ = 1; see also the discussion in [40].
3. Implicit sampling. Implicit sampling is a general numerical method for constructing effective importance functions [36,3,14,13]. We write the pdf to sample as f = e −F (x) and, temporarily, assume that F (x) is convex. The idea is to write x as a one-to-one and onto function of an easy-to-sample reference variable ξ with pdf g ∝ e −G(ξ) , where G is convex. To find this function, we first find the minima of F, G, call them φ F = min F, φ G = min G, and define φ = φ F − φ G . Then define a mapping ψ : ξ → x such that ξ and x satisfy the equation With our convexity assumptions a one-to-one and onto map exists, in fact there are many, because equation (3) is underdetermined (it is a single equation that connects the m elements of X i to the m elements of Ξ i , where m is the number of variables). To find samples X 1 , X 2 , . . . , first find samples Ξ 1 , Ξ 2 , . . . of ξ, and for each i = 1, 2, . . . , solve (3). A sample X i obtained this way has a high probability (with a high probability) for the following reasons. A sample Ξ i of the reference density g is likely to lie near the minimum of G, so that the difference (G(Ξ i ) − φ G ) is likely to be small. By equation (3) the difference (F (X i ) − φ F ) is also likely to be small and the sample X i is in the region where f is large. Thus, by solving (3), we map the high-probability region of the reference variable ξ to the high-probability region of x, so that one obtains a high-probability sample X i from a high-probability sample Ξ i . The mapping ψ defines the importance function implicitly by the solutions of (3): A short calculation shows that the sampling weight w(X i ) is where J is the Jacobian det(∂x/∂ξ). The factor e −φ , common to all the weights, is immaterial here, but plays an important role further below. We now give an example of a map ψ that makes this construction easy to implement. We assume here and for the rest of the paper that the reference variable ξ is the Gaussian N (0, I m ), where m is the dimension of x. With a Gaussian ξ, φ G = 0 and equation (3) becomes We can find a map that satisfies this equation by looking for solutions in a random direction η = ξ/(ξ T ξ), i.e. use a mapping ψ such that where µ = argmin F is the minimizer of F , L is a given matrix, and λ is a scalar that depends on ξ. Substitution of the above mapping into (5) gives a scalar equation in the single variable λ (regardless of the dimension of x). The matrix L can be chosen to optimize the distribution of the sampling weights. For example, if L is a square root of the inverse of the Hessian of F at the minimum, then the algorithm is affine invariant, which is important in applications with multiple scales [23]. Equation (6) can be readily solved and the resulting Jacobian is easy to calculate (see [36] for details).
Other constructions of suitable maps ψ are presented in e.g. [13,22] and some of these are analyzed in [22]. With these constructions, often the most expensive part of the calculation is finding the minimum of F . Note that this needs to be done only once for each sampling problem, and that finding the maximum of f is an unavoidable part of any useful choice of importance function, explicitly or implicitly. Addressing the need for the optimization explicitly has the advantage that effective optimization methods can be brought to bear. Connections between optimization and (implicit) sampling also have been pointed out in [3,37]. In fact, an alternate derivation of implicit sampling can be given by starting from maximum likelihood estimation, followed by a randomization that removes bias and provides error estimates.
We now relax the assumption that F is convex. If F is U -shaped, the above construction works without modification. A scalar function F is U -shaped if it is piecewise differentiable, its first derivative vanishes at a single point which is a minimum, F is strictly decreasing on one side of the minimum and strictly increasing on the other, and F (x) → ∞ as |x| → ∞; in the general case, F is U -shaped if it has a single minimum and each intersection of the graph of the function y = F (x) with a vertical plane through the minimum is U -shaped in the scalar sense. If F is not U -shaped, but has only one minimum, one can replace it by a U -shaped approximation, say F 0 , and then apply implicit sampling as above. The bias created by this approximation can be annulled by reweighting [13]. If there are multiple minima, one can represent f as a mixture of components whose logarithms are Ushaped, and then pick as a reference pdf g a cross-product of a Gaussian and a discrete random variable. However, the decomposition into a suitable mixture can be laborious (see [54]). 4. Beyond implicit sampling. Implicit sampling produces a sequence of samples that lie in the high-probability domain and are weighted by a Jacobian. It is natural to wonder whether one can absorb the Jacobian into the map ψ : ξ → x and obtain "optimal" samples that need no weights (here and below, optimal refers to minimum variance of the weights, i.e. an optimal sampler is one with a constant weight function). If a random variable x with pdf f is a one-to-one function of a variable ξ with pdf g, then where J is the Jacobian of the map ψ. One can obtain samples of x directly (i.e. without weights) by solving (7). Notable efforts in that direction can be found in [43,38].
However, optimal samplers can be expensive to use; the difficulties can be demonstrated by the following construction. Consider a one-variable pdf f with a convex F = − log f . Find the maximizer z of f , with maximum M f = f (z), and let g be the pdf of a reference variable ξ, with its maximizer also at z and maximum M g = g(z). To simplify the notations, change variables so that z = 0. Then construct a mapping ψ : ξ → x as follows. Solve the differential equation where t is the independent variable and u is a function of t, in the t-interval [0, ξ], with initial condition u = 0 when t = 0. Define the map from ξ to x by x = u(ξ); then f (x) = g(ξ)J, where J = |du/dx| evaluated at t = ξ. Under the assumptions made, this map is one-to-one and onto. We have achieved optimal sampling. The catch is that the initial value of g(t)/f (u) is M g /M f , which requires that the normalization constants of f and g be known. In general the solution of the differential equation does not depend linearly on the initial value of g(0)/f (u(0)), so that unless M f and M g are known, the samples are in the wrong place while weights to remove the resulting bias are unavailable. Analogous conclusions were reached in [17] for a different sampling scheme. In the special case where F = − log f and G = − log g are both quadratic functions the constants M f , M g can be easily calculated, but under these conditions one can find a linear map that satisfies equation (7) and implicit sampling becomes identical to optimal sampling. The elimination of all variability in the weights in the nonlinear case is rarely worth the cost of computing the normalization constants.
Note that the normalization constants are not needed for implicit sampling. The mapping (6) for solving (3) is well-defined even when the normalization constants of f and g are unknown, because if one multiplies the functions f or g by a constant, the logarithm of this constant is added both to F (G) and to φ F (φ G ) and cancels out. Implicit sampling simplifies equation (7) by assuming that the Jacobian is a constant. It produces weighted samples where the weights are known only up to a constant, and this indefiniteness can be removed by normalizing the weights so that their sum equals 1.

5.
Bayesian parameter estimation. We now apply implicit sampling to Bayesian estimation. For the sake of clarity we first explain how to use Bayesian estimation to determine physical parameters needed in the numerical modeling of physical processes, given noisy observations of these processes. This discussion explains how a prior and data combine to create a posterior probability density that is used for inference. In the next section we extend the methodology to data assimilation, which is our main focus.
Suppose one wishes to model the diffusion of some material in a given medium. The density of the diffusing material can be described by a diffusion equation, provided the diffusion coefficients θ ∈ R m can be determined. Given the coefficients θ, one can solve the equation and determine the values of the density at a set of points in space and time. The values of the density at these points can be measured experimentally, by some method with an error η whose probability density is assumed known; once the measurements d ∈ R k are made, this assigns a probability to θ. The relation between d and θ can be written as where h : R k → R m is a generally nonlinear function, and η ∼ p η (·) is a random variable with known pdf that represents the uncertainty in the measurements. The evaluation of h often requires a solution of a differential equation and can be expensive. In the Bayesian approach, one assumes that information about the parameters is available beforehand in form of a prior probability density p 0 (θ). This prior and the likelihood p(d|θ) = p η (d − h(θ)), defined by (9), are combined via Bayes' rule to give the posterior pdf, which is the probability density function for the parameters θ given the data d, where γ(d) = p 0 (θ)p(d|θ)dθ is a normalization constant (which may be hard to compute). In a Monte Carlo approach to the determination of the parameters, one finds samples of the posterior pdf. These samples can, for example, be used to approximate the expected value which is the best approximation of θ in a least squares sense (see, e.g. [9]); an error estimate can also be computed, see [48].
When sampling the posterior p(θ|d), it is natural to set the importance function equal to the prior probability p 0 (θ) and then weight the samples by the likelihood p(d|θ). However, if the data are informative, i.e., if the measurements d modify the prior significantly, then the prior and the posterior can be very different and this importance sampling scheme is ineffective. When implicit sampling is used to sample the posterior, the data are taken into account already when one generates the samples, not only in the weights of the samples.
As the number of variables increases, the prior p 0 and the posterior p(θ|d) may become nearly mutually singular. In fact, in most practical problems these pdfs are negligible outside a small spherical domain so that the odds that the prior and the posterior have a significant overlap decrease quickly with the the number of variables, making implicit sampling increasingly useful (see [37] for an application of implicit sampling to parameter estimation).
6. Data assimilation. We now apply Bayesian estimation via implicit sampling to data assimilation, in which one makes inferences from an unreliable time-dependent approximate model that defines a prior, supplemented by a stream of noisy and/or partial data. We assume here that the model is a discrete recursion of the form where n is a discrete time, n = 0, 1, 2, . . . , the set of m-dimensional vectors x 0:n = (x 0 , x 1 , . . . , x n , . . . ) are the state vectors we wish to estimate, f (·) is a given function, and w n is a Gaussian random vector N (0, Q). The model often represents a discretization of a stochastic differential equation. We assume that k-component data vectors b n , n = 1, 2, . . . are related to the states x n by where h(·) is the (generally nonlinear) observation function and η n is a N (0, R) Gaussian vector (the η n and w n are independent of η k and w k for k < n and also independent of each other). In addition, initial data x 0 , which may be random, are given at time n = 0. For simplicity, we consider in this section the problem of estimating the state x of the system rather than estimating parameters in the equations as in the previous section. Thus, we wish to sample the conditional pdf p(x 0:n |b 1:n ), which describes the probability of the sequence of states between 0 and n given the data b 1:n = (b 1 , . . . , b n ). The samples of this conditional pdf are sequences X 0 , X 1 , X 2 , . . . , usually referred to as "particles". The conditional pdf satisfies the recursion where p(x n+1 |x n ) is determined by equation (11) and p(b n+1 |x n+1 ) is determined by equation (12). (see e.g. [16]). We wish to sample the conditional pdf recursively, time step after time step, which is natural in problems where the data are obtained sequentially, and drastically reduces the computational cost. To do this we use an importance function π of the form where the π k are increments of the importance function. The recursion (13) and the factorization (14) yield a recursion for the weight W n+1 j of the j-th particle, With this recursion, the task at hand is to pick a suitable update π k (x k |x 0:k−1 , b 1:k ) for each particle, to sample p(X n+1 j |X n j )p(b n+1 |X n+1 j ), and to update the weights so that the pdf p(x 0:n+1 |b 1:n+1 ) is well-represented. The solution of the discrete model (11) plays the same role in data assimilation as the prior in the parameter estimation problem of section 5.
In the SIR (sequential importance sampling) filter one picks . Thus, the SIR filter proposes samples by the model (11) and determines their probability from the data (12).
In implicit data assimilation one samples π n+1 (x n+1 |x 0:n , b 1:n+1 ) by implicit sampling, thus taking the most recent data b n+1 into account for generating the samples. The weight of each particle is given by equation (15) W n+1 . A minimization is required for each particle at each step.
The "optimal" importance function [17,56,28] is defined by (14) with and its weights are W n+1 j = W n j p(b n+1 |X n j ). This choice of importance functions is optimal in the sense that it minimizes the variance of the weights at the n + 1 step given the data and the particle positions at the previous step. In fact, this importance function is optimal over all importance functions that factorize as in (14), and in the sense that the variance of the weights var(w n ) is minimized (with expectations computed with respect to π(x 0:n+1 |b 0:n+1 )), see [47]. In a linear problem where all the noises are Gaussian, or in a problem with nonlinear dynamics and a linear observation functional in which observations are available every time step, the implicit filter and the optimal filter coincide [14,36,32]. When a problem is nonlinear, the optimal filter may be hard to implement, as discussed above.
A variety of other constructions is available (see e.g. [54,53,52,51,1]). The SIR and optimal filters are extreme cases, in the sense that one of them makes no use of the data in finding samples while the other makes maximum use of the data. The optimal filter becomes impractical for nonlinear problems while the implicit filter can be implemented at a reasonable cost. Implicit sampling balances the use of data in sampling and the computational cost.
7. The size of a covariance matrix and the feasibility of data assimilation. Particle filters do not necessarily succeed in assimilating data. There are many possible causes-for example, the recursion (11) may be unstable, or the data may be inconsistent. The most frequent cause of failure is filter collapse, in which only one particle has a non-negligible weight, so that there is effectively only one particle, which is not sufficient for valid inference. In the next section we analyze how particle collapse happens and how to avoid it. In preparation for this discussion, we discuss here the Frobenius norm of a covariance matrix, its geometrical significance, and physical interpretation.
The effectiveness of a sampling algorithm depends on the pdf that it samples. For example, suppose that the posterior you want to sample is in one variable, but has a large variance. Then there is a large uncertainty in the state even after the data are taken into account, so that not much can be said about the state. We wish to assess the conditions under which the posterior can be used to draw reliable conclusions about the state, i.e. we make the statement "the uncertainty is not too large" quantitative. We call sampling problems with a small enough uncertainty "feasible in principle". There is no reason to attempt to solve a sampling problem that is not feasible in principle.
Our analysis is linear and Gaussian and we allow the number of variables to be large. In multivariate Gaussian problems, feasibility requires that a suitable norm of the covariance matrix be small enough. This analysis is inspired by geophysical applications where the number of variables is large, but the nonlinearity may be mild; parts of it were presented in [12].
We describe the size of the uncertainty in multivariate problems by the size the m × m covariance matrix P = (p ij ), which we measure by its Frobenius norm where the λ k , k = 1, . . . , m are the eigenvalues of P . The reason for this choice is the geometric significance of the Frobenius norm. Let x ∼ N (µ, P ) be an mdimensional Gaussian variable, and consider the distance between a sample of x and its most likely value µ

If all eigenvalues of P are O(1) and if m is large, then E[r]
, the expected value of r, is O(||P || F ) and var(r) = O(1), i.e. the samples concentrate on a thin shell, far from their most likely value. Different parts of the shell may correspond to physically different scenarios (such as "rain" or "no rain"). Moreover, the expected value and variance of the distance distance r = √ y of a sample to its most likely value are bounded above and below by multiples of ||P || F (see [12]). If one tries to estimate the state of a system with a large ||P || F , the Euclidean distance between a sample and the resulting estimate is likely to be large, making the estimate useless. For a feasible problem we thus require that ||P || F not be too large. How large "large" can be depends on the problem. In [12] and in earlier work [46,6,4] on the feasibility of particle filters, the Frobenius norm of P was called the "effective dimension". This terminology is confusing and we abandon it here. The norm ||P || F quantifies the strength of the noise, but we define the effective dimension of the noise by min ∈ Z : where is a predetermined small parameter. This definition of effective dimension is similar in spirit to the standard procedure in regression analysis in which an effective number of independent variables is determined by applying a formula essentially identical to (17) to the covariance matrix of the independent variables (e.g., [20]). The small parameter is usually chosen to retain only those components with amplitudes above some predetermined noise level. A small ||P || F can imply a small , however the reverse is not necessarily true and ||P || F can be small, even though is large (for example if x ∼ N (0, I m / √ m)) and vice versa. There are low-dimensional problems with a large noise (large ||P || F ) which are not feasible in principle, and there are large dimensional problems with a small noise which are feasible in principle (small ||P || F ). Estimation based on a posterior pdf is feasible in principle only if ||P || F is small. We now examine the relations between ||P || F , the effective dimension , and the correlations between the components of the noise (i.e., the size of the off-diagonal terms in P ). We assume that our random variables are point values x i of a stationary stochastic process on the real line, measured at the points ih of the finite interval [0, 1], where i = 1, 2, . . . , m and h = 1/m. As the mesh is refined, m increases. The condition for the discrete problem to be feasible is that the Frobenius norm of its covariance matrix be bounded. Let the covariance matrix of the continuous process be k(x, x ) = k(x − x ); the corresponding covariance operator is defined by for every function ϕ = ϕ(x). The eigenvalues of K can be approximated by the eigenvalues of P multiplied by h (see [42]). The Frobenius norm of K, , describes the distance of a sample to its most likely value in the L 2 sense, as can be seen by following the same steps as in [12], replacing the Euclidean norm by the appropriate grid-norm ||x|| 2 2 = m k=1 h x 2 k . We consider a family of Gaussian stationary processes with differing correlation structures. The discussion is simplified if we keep the energy e defined by e = ∞ −∞ k(s) 2 ds, (s = x − x ) equal to 1 for all the members of the family (so that all the members of the family are equally noisy); that e is an energy follows from Khinchin's theorem (see e.g. [9]). An example of such a family is the family of zero-mean processes with where L is the correlation length. The infinite limits in the definition of e make the calculations easier, and are harmless as long as L is less than about 1/3. The elements p ij of the discrete matrix P are p ij = k(ih, jh), i, j = 1, . . . , m.
In the left panel of figure 3 we demonstrate that for a fixed correlation length L, the Frobenius norm where the λ k are the eigenvalues of the covariance matrix P , is independent of the mesh (for small enough h) and, therefore independent of the discretization. In the calculation of the dimension in (17), we use = 5%. In the figure, we plot the Frobenius norms of the m × m matrices P as the purple dashed lines. The solid line is the Frobenius norm of the covariance operator K of the continuous process, which, for this simple example, can be computed analytically: where theλ are the eigenvalues of the operator K. We find a good agreement between the infinite dimensional results for K and and their finite dimensional approximations by ||P || F . (the dashed lines are mostly invisible because they coincide with the results calculated from the infinite dimensional problem). The right panel of the figure shows the variation of ||P || F and the dimension with the correlation length L. We observe that the Frobenius norm remains almost constant for L < 10 −2 . On the other hand, the dimension in (17) increases as the correlation length decreases. What this figure shows is that the feasibility of data assimilation depends on the level of noise, but not on the effective number of variables. One can find the same level of noise for very different effective dimensions. If data assimilation is not feasible, it is because the problem has too much noise, not too many variables.
It is interesting to consider the limit L → 0 in our family of processes. A stationary process u(x) such that for every x, u(x) is a Gaussian variable, with E[u(x)u(x )] = 0 for x = x and E[u(x) 2 ] = A, where A is a finite constant, has very little energy; white noise, the most widely used Gaussian process with independent values, has E[u(x)u(x )] = δ x,x , where the right-hand-side is a δ function, so that E[u(x) 2 ] is unbounded. The energy of white noise is infinite, as is well-known. If k(x − x ) in our family blew up like L −1 at the origin as L → 0, the process would be a multiple of Brownian motion; here it blows up like L −1/2 , allowing the energy to remain finite while not allowing the E[u 2 ] to remain bounded. The moral of this discussion is that a sampling problem where P = I m for all m is unphysical.
An apparent paradox appears when one applies our theory to determine the feasibility of a Gaussian random variable with covariance matrix P = I m , as in [46,6,4,45,47]. In this case, ||P || F = √ m, so that the problem is infeasible by our definition if the number of variables m is large. On the other hand, the problem seems to be physically plausible. For example, suppose one wishes to estimate the wind velocity at m geographical sites based on independent local models for each site (e.g. today's wind velocity is yesterday's wind velocity with some uncertainty), and independent noisy measurements at each site. The local P j at each site is onedimensional and equals 1. Why then not try to determine the m velocities in one fell swoop, by considering the whole vector x and setting P = I m ? Intuition suggest that the set of local problems is equivalent to the P = I m problem, e.g. that one can use local stochastic models to predict the velocities at nearby sites, and then mark these on a weather map and obtain a plausible global field. Thus, ||P || F seems to be large in a plausible and feasible problem. This, however, is not so. First, in reality, the uncertainties in the velocities at nearby sites are highly correlated, while the problem P = I m assumes that this is not so. The resulting velocities map from the latter will be unphysical and lead to wrong forecasts-large scale coherent flows in today's weather map will have a different impact on tomorrow's forecast than a set of uncorrelated wind patterns, and carry a much larger energy, even if the local amplitudes are the same. If one has a global problem where the covariance matrix is truly I m , then replacing the solution of the full problem by a component-wise solution merely changes the estimate of the noise from ||P || F to the numerically smaller maximum of the component variances, which is not as good a measure of the distance between the samples and their mean.
8. Linear analysis of the convergence of data assimilation. We now summarize the linear analysis of the conditions under which particular particle filters produce reliable estimates when the problem is linear (see [12]). A general nonlinear analysis is not within reach, while the linear analysis captures the main issues. The dynamical equation now takes the form: where A is an m × m matrix and, the data equation becomes where H is an k × m matrix. As before, w n are independent N (0, Q) and the η n are independent N (0, R), independent also of the w n . We assume that equation (18) is stable and the data are sufficient for assimilation (for a more technical discussion of these assumptions, see e.g. [12]). These two equations define a posterior probability density, independently of any method of estimation. The first question is, under what conditions is state estimation with the posterior feasible in principle. If one wants to estimate the state at time n, given the data up to time n, this requires that the covariance of x n |b 1:n have a small Frobenius norm. The theoretical analysis of the Kalman filter [25] can be used to estimate this covariance, even when the problem is not feasible in principle or the Kalman filter itself is too expensive for practical use. Under wide conditions, the covariance of x n |b 1:n rapidly reaches a steady state where X is the unique positive semi-definite solution of a particular nonlinear algebraic Riccati equation (see e.g. [27]) and One can calculate ||P || F and decide whether the estimation problem is feasible in principle. Thus, a condition for successful data assimilation is that ||P || F be moderate, which generaly requires that the ||Q|| F , ||R|| F in equations (18) and (19) be small enough.
A particle filter can be used to estimate the state by sampling the posterior pdf p(x 0:n |b 1:n ), defined by (18) and (19). The state at time n conditioned on the data up to time n can be computed by marginalizing samples of p(x 0:n |b 1:n ) to samples of p(x n |b 1:n+1 ) (which amounts to simply dropping the sample's past). Thus, a particle filter does not directly sample p(x n |b 1:n ), so that its samples carry weights (even in linear problems). These weights must not vary too much or else the filter "collapses" (see section 2). Thus, a particle filter can fail, even if the estimation problem is feasible in principle.
It was shown in [12,46,6,4,45] that the variance of the negative logarithm of the weight must be small or else the filter collapses. For the SIR filter, the condition that this collapse not happen is that the Frobenius norm of the matrix be small enough. For the optimal filter, which in the present linear setting coincides with the implicit filter, success requires a small Frobenius norm for see [12]. In either case, this additional condition must be satisfied as well as the condition that ||P || F be small. To understand what these formulas say, we apply them to a simple model problem which is popular as a test problem in the analysis of numerical weather prediction schemes. The problem is defined by H = A = I m and Q = qI m , R = rI m , where we vary the number of components m. Note that for a fixed m, the problem is parametrized by the variance r of the noise in the model and the variance q of the noise in the data. This problem is feasible in principle if is not too large. For a fixed m, this means that where the number 1 stands in for a sharper estimate of the acceptable variance for a given problem (and this choice gives the complete qualitative story). In figure 4, this condition is illustrated by showing feasible problems in white and all infeasible problems in grey. The analysis thus shows that for data assimilation to be feasible, either the model or the data must be accurate. If both are noisy, data assimilation is useless.
In the SIR filter, the additional condition for success becomes and for the optimal/implicit case, the additional condition is In both cases, the number 1 on the right-hand-side stands in for a sharper estimate of the acceptable variance of the weights for a given problem (which also depends on the computing resources). In both cases, the added condition is quadratic and homogeneous in the ratio q/r, and thus slices out conical regions from the region where data assimilation is feasible in principle. These conical regions are shown in figure 4. The region where estimation with a particular filter is feasible is labeled by the name of that filter. Note that we also show results for the ensemble Kalman filter (EnKF) [18,49], which are derived in [35], but not discussed in the present paper. The analysis of the EnKF relates to the situation where it is used to sample  . Conditions for feasible data assimilation and for successful sampling with the optimal particle filter (which coincides here with the implicit filter), the SIR particle filter, and the ensemble Kalman filter (from [35]). the same posterior density as the other filters quoted in the figure. The analysis in [35] also includes the situation where the EnKF is used to sample a marginal of that pdf directly, when the conclusions are different.
In practical problems one usually tries to sharpen estimates obtained from approximate dynamics with the help of accurate measurements, and the optimal/ implicit filter works well in such problems. However, there is a region in which not even the optimal filter succeeds. What fails there is the sequential approach to the estimation problem, i.e. the factorization of the importance function as in (14), see [35,47]. Non-recursive filters can succeed in this region. However, one can use filters that are not recursive, and there too implicit sampling can be helpful i.e. one can apply implicit sampling to p(x 0:n |b 1:n ).
One conclusion from our analysis is that if Bayesian estimation is feasible in principle then it is feasible in practice. However, it is assumed in the previous analysis so far that one has suitable priors, or equivalently, one knows what the noise w n in equation (11) really is. The fact is that generally one does not. We discuss this issue in the next section. 9. Estimating the prior. The previous sections described how one may use a prior and a distribution of observation errors to estimate the states or the parameters of a system. This estimation depends on knowing the prior. As we have seen, in data assimilation, the prior is generated by the solution of an equation such as the recursion (11) and depends on knowing the noise w n in that equation. It is often assumed that the noise in that equation is white. However, one can show that the noise is not white in most problems (see e.g. [55,15,10]). We now present a preliminary discussion of methods for estimating the noise.
First, one has to make some assumptions about the origin of the noise. A reasonable assumption (see e.g. [10,5]) is that there is noise in the equations of motion because a practical description of nature is necessarily incomplete. For example, one can write a solvable equation of motion for a satellite turning around the earth by assuming that the gravitational pull is due to a perfectly spherical earth with a density that is a function only of distance from the center (see [34]). Reality is different, and the difference produces noise, also known as model error. The problems to be solved are: (i) estimate this difference, (ii) identify it, i.e., find a concise approximate representation of this difference that can be effectively evaluated or sampled on a computer, and (iii) design an algorithm that imbeds the identified noise in a practical data assimilation scheme. These problems have been discussed in a number of recent papers, e.g. [10,31], in particular the review paper [24] which contains an extensive bibliography.
Assume that the full description of a physical system has the form: where t is the time, x = (x 1 , x 2 , . . . , x m ) is the vector of resolved variables, and y = (y 1 , y 2 , . . . , y ) is the vector of unresolved variables, with initial data x(0), y(0). Consider a situation where this system is too complicated to solve, but where data are available, usually as sequences of measured values of x, assumed here to be observed with negligible observation errors. Write R(x, y) in the form where R 0 is a practical approximation of R(x, y) and one is able to solve the equation However, x does not satisfy equation (25) because the true equation, the first of equations (23), is more complicated. The difference between the true equation and the approximate equation is the remainder z(x, y) = R(x, y) − R 0 (x), which is the noise in the determination of x. In general, z has to be determined from the data, i.e., from the observed values of x.
A usual approach to the problem of estimating z is to use equation (24) to obtain its values from x data, i.e. calculate z = d dt x − R 0 (x), and then identify it as a continuous stochastic process. This may be difficult, in particular, calculating z requires that one differentiate x, which is generally impractical or inaccurate because z may have high-frequency components or fail to be sufficiently smooth, and the data may not be available at sufficiently small time intervals (an illuminating analysis in a special case can be found in [41,44]). Once one has values of z, identifying it as a function of the continuous variable t often requires making unwarranted assumptions on the small-scale structure of z, which may be unknown when the data are available only at discrete times.
An alternative is supplied by a discrete approach as in [11]. Equation (25) is always solved on the computer, i.e., in discrete form, the data are always given at discrete points, and it is x one wishes to determine but in general one is not interested in determining z per se. We can therefore avoid the difficult detour through a continuous-time z followed by a discretization, as follows. We pick once and for all a particular discretization of equation (25) with a particular time step time step δ, where R δ is obtained, for example, from a Runge-Kutta method, and where n indexes the result after n steps; the differential equation has been reduced to a recursion such as equation (11). We then use the data to identify the discrepancy sequence, z n+1 δ = (x n+1 − x n )/δ − R δ (x n ), which is available from x data without approximation.
Then assume, as one does in the continuous-time case, that the system under consideration is ergodic, so that its long-time statistics are stationary. The sequence z n δ becomes a stationary time series, which can be represented by one of the representations of time series, e.g. the NARMAX (nonlinear auto-regression moving average with exogenous inputs) representation, with x as an exogenous input. The observed x of course may include observation errors, which have to be separated from the model noise via filtering. A similar approach was presented in [39,33] for the Lorenz '63 model [29]. In [8] an autoregressive process was fitted to the time series of model-data misfits in a linear model of waves in the tropical Pacific in order to construct an optimized Kalman filter.
As an example, we applied this construction to the Lorenz 96 system [30], created to serve as a metaphor for the atmosphere, which has been widely used as a test bench for various reduction and filtering schemes. It consists of a set of chaotic differential equations, which, following [19], we write as: d dt y j,k = 1 ε [y j+1,k (y j−1,k − y j+2,k ) − y j,k + h y x k ] with z k = h x J j y j,k , and k = 1, . . . , K, j = 1, . . . , J. The indices are cyclic, x k = x k+K , y j,k = y j,k+K and y j+J,k = y j,k+1 . This system is invariant under spatial translations, and the statistical properties are identical for all x k . The parameter ε measures the timescale separation between the resolved variables x k and the unresolved variables y j,k . The parameters are ε = 0.5, K = 18, J = 20, F = 10, h x = −1 and h y = 1. The ergodicity of the Lorenz 96 system has been established numerically in earlier work (see e.g. [19]). We pretend that this system is too difficult to integrate in full; we take R 0 in equation (25) as the right hand side of equation (27) with z k = 0. The noise is then z k , which in our special case can actually be calculated by solving the full system of equations; in general the noise has to be estimated from data as described above and in [55,11]. In figure 5 we plot the covariance function of the noise component z 1 determined as just described. Note that z is far from white. To the extent that the Lorenz 96 is a valid metaphor for the atmosphere, we find that the noise in a realistic problem is not white. This can of course be deduced from the construction after equations (23), where the noise appears as the difference between solutions of differential equations, which is not likely to be white. This relation between the preceding discussion of the noise (which defines the prior through equation (24)) should be compared with the discussion of implicit sampling and particle filters. Here too the data have been put to additional use (there to define samples and not only to weight them, here to provide information about the noise as well as about the signal); here too an assumption of a white noise input has been found wanting, and realism requires non-trivial correlations, (there in space, here in time).
The next steps are to extract noise estimates from noisy observations of a system, and then combine the result with the observations to produce a state estimate. Recent work on this topic is summarized in [24]. An example of this construction could not be produced in time to appear in this article. 10. Conclusions. We have presented and analyzed algorithms for data assimilation and for estimating the prior. The work presented shows that the keys to success are a better use of the data and a careful analysis of what is physically plausible and useful. We feel that significant progress has been made. One conclusion we draw from this work is that data assimilation, which is often considered as a problem in statistics, should also, or even mainly, be viewed as a problem in computational physics. This conclusion has also been reached by others, see e.g. [31,15].