# Large-scale Kalman filtering using the limited memory BFGS method.

1. Introduction. The Kalman filter (KF) for linear dynamical systems and the extended Kalman filter (EKF) for nonlinear but smoothly evolving dynamical systems are popular methods for use on state space estimation problems. As the dimension of the state space becomes very large, as is the case, for example, in numerical weather forecasting, the standard formulations of KF and EKF become computationally intractable due to matrix storage and inversion requirements.Computationally efficient variants of KF and EKF have been proposed for use on such large-scale problems. The Reduced Rank Kalman Filter or Reduced Order extended Kalman filter (see, e.g., [4, 7, 30]) project the dynamical state vector of the model onto a lower dimensional subspace. The success of the approach depends upon a judicious choice of the reduction operator. Moreover, since the reduction operator is typically fixed in time, the dynamics of the system may not be correctly captured; see [9] for more details.

In the context of numerical weather forecasting, a great deal of attention has been given to the filtering problem. The current state of the art is 4D-Var (see, e.g., [9, 23]), which utilizes a variational formulation of an initial value estimation problem [11, 14, 17]. 4D-Var has been shown to be identical to a Kalman smoother when the model is assumed to be perfect [16]. The resulting quadratic minimization problem is very large-scale ([10.sup.4]-[10.sup.7] unknowns) and so efficient numerical optimization methods are needed. Similar to the methods in the previous paragraph, the partial orthogonal decomposition is used in [5] to reduce the dimensionality of the 4D-Var minimization problem. A more standard approach is to implement a preconditioned conjugate gradient method [8, 12, 22, 26]. In this context, a number of different preconditioners have been tested.

In this paper, we take a different approach. In particular, we focus our attention on the Kalman filter itself, using the limited memory BFGS (LBFGS) [22] iterative method for the required large-scale matrix storage and inversion within KF and EKF. More specifically, suppose Ax = b is a system, with symmetric positive definite matrix A, that requires solution, and or a low storage approximation of [A.sup.-]1 is needed, within the KF or the EKF algorithms. In such cases, LBFGS is a natural choice, since it generates both a sequence of approximations of [A.sup.-]1b, and a sequence of symmetric positive definite matrices {[B.sup.-1.sub.k]} approximating [A.sup.-1]. This choice is further supported by the result that [A.sup.-1]b is reached within a finite number of LBFGS iterations (assuming exact arithmetic) [20, 21], and that if k and full memory BFGS is used [B.sup.-1.sub.n+1] = [A.sup.-1] [22, Chapter 8]. A result regarding the accuracy of the LBFGS Hessian approximations to [A.sup.-1] is also given in [21]; it depends upon the eigenvalues of [A.sup.1/2] [B.sup.k] [A.sup.1/2]; see Appendix A. The use of LBFGS and its success in the examples that we consider can be motivated by the fact that the covariance matrices being approximated are approximately low rank. In applications of interest [6], covariance information is contained in slowly varying, low dimensional subspaces, making accurate low-rank approximations possible.

The idea of using the LBFGS method in variational data assimilation is not new; see, e.g., [10, 13, 17, 25, 27, 28, 29, 31]. In many of these references, the LBFGS Hessian or inverse Hessian is used as a preconditioner for conjugate gradient iterations, and even as an approximate error covariance matrix for the background term in 3D- and 4D-Var variational data assimilation. However, in the method presented here, the LBFGS method is further used for matrix inversion, in order to propagate effectively the state estimate covariance information forward in time. Moreover, we apply our methodology to the Kalman filter itself, not to the variation formulation used by the 3D- and 4D-Var methods [17]. LBFGS can also be incorporated in a fully variational formulation of the Kalman filter; see [2].

As has been stated, the approach presented here uses the LBFGS algorithm directly within the context of the Kalman filter. The equivalence of LBFGS and a certain preconditioned conjugate gradient method (see [21] and Appendix A) suggests that our approach and those cited above are similar. One advantage of the cited approaches, however, is that they can be incorporated into existing 3D- and 4D-Var codes used in practice.

An application of a similar methodology that could be used in conjunction with 3Dand 4D-Var is presented in [1, 2]. The aim of the current paper is to demonstrate the use of LBFGS within the standard (non-variational) formulation of the linear or extended Kalman filter.

The paper is organized as follows. We present KF and EKF in Section 2, and then in Section 3 we present LBFGS-KF and LBFGS-EKF. We test these methods with two numerical experiments in Section 4. Conclusions are then given in Section 5, and implementation details of the LBFGS algorithm are contained in Appendix A.

2. The Kalman filter. We consider the coupled system of discrete, linear stochastic difference equations given by

(2.1) [x.sub.k] = [M.sub.k][x.sub.k-1] + [[epsilon].sup.p.sub.k],

(2.2) [x.sub.k] = [K.sub.k][x.sub.k-1] + [[epsilon].sup.o.sub.k]

In the first equation, [x.sub.k] denotes the n x 1 state of the system at time k; [M.sub.k] is the n x n linear evolution operator; and [[epsilon].sup.o.sub.k] is a n x 1 random vector known as the prediction error and is assumed to characterize errors in the model and corresponding numerical approximations. In the second equation, [y.sub.k] denotes the m x 1 observed data; [K.sub.k] is the m x n linear observation operator; and [[epsilon].sup.o.sub.k] is an m x 1 random vector known as the observation error. The prediction error [[epsilon].sup.p.sub.k] and observation error [[epsilon].sup.o.sub.k] are assumed to be independent and normally distributed, with zero means and symmetric positive definite covariance matrices [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], respectively.

We assume, in addition, that we have in hand estimates of both the state [x.sup.est.sub.k-1] and its positive definite covariance matrix [C.sup.est.sub.k-1] at time k - 1. Moreover, we assume that [x.sup.est.sub.k-1], [[epsilon].sup.p.sub.k], and [[epsilon].sup.p.sub.k] are independent random vectors. The goal is to then estimate [x.sup.est.sub.k] and its covariance [C.sup.est.sub.k-1]. In Bayesian terms, equation (2.2) provides the likelihood function in the estimation step, while (2.1) gives the prior. By the assumptions on [[epsilon].sup.p.sub.k] and [[epsilon].sup.o.sub.k], the prior mean and the prior covariance matrix are directly obtained from ( 2.1),

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The full negative-log posterior density given (up to an additive constant) by Bayes' theorem, takes the form

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII],

and hence, we have

(2.3) [x.sup.est.sub.k] = arg min l(x|[y.sub.k]),

(2.4) [C.sup.est.sub.k] = [[nabla].sup.2] l [(x|[y.sub.k]).sup.-1].

Equations (2.3) and (2.4) motivate the variational Kalman filter, which is the subject of [1]. However, they can also be used to derive the Kalman filter. In particular, noting that (2.3) and

(2.4) can be alternatively written (see [24] for detail)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII],

we have the following standard formulation of the Kalman filter.

The Kalman Filter

Step 0: Select initial guess [x.sup.est.sub.0] and covariance [C.sup.est.sub.0], and set k = 0.

Step 1: Compute the evolution model estimate and covariance:

(i) Compute [x.sup.p.sub.k] = [M.sub.k][x.sup.est.sub.k-1];

(ii) Compute [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Step 2: Compute the Kalman filter estimate and covariance:

(i) Compute the Kalman Gain [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

(ii) Compute the Kalman filter estimate [x.sup.est.sub.k] = [x.sup.p.sub.k] + [G.sub.k] ([y.sub.k] - [K.sub.k] [x.sup.p.sub.k);

(iii) Compute the estimate covariance [C.sup.est.sub.0] = [C.sup.p.sub.k] - [G.sub.k] [K.sub.k] [C.sup.p.sub.k].

Step 3: Update k := k + 1 and return to Step 1.

Note that it is typical to take the initial covariance [C.sup.est.sub.k] to be diagonal.

The extended Kalman filter (EKF) is the extension of KF when (2.1), (2.2) are replaced by

(2.5) [x.sub.k] = M([x.sub.k-1]) + [[epsilon].sup.p.sub.k],

(2.6) [x.sub.k] = K([x.sub.k]) + [[epsilon].sup.0.sub.k],

where M and K are (possibly) nonlinear functions. EKF is obtained by the following simple modification of the above algorithm: in Step 1 (i), use instead [x.sub.p] = M ([x.sup.est.sub.k]), and define

(2.7) [M.sub.k] [partial derivative]M([x.sup.est.sub.k-1]/[partial derivative]x, and [K.sub.k] [partial derivative]K([x.sup.p])/[partial derivative]x.

We note that [M.sub.k] and [K.sub.k] can be computed or estimated in a number of ways. For example, the numerical scheme that is used in the solution of either the evolution or the observation model defines a tangent linear code, which can be used to compute (2.7); see, e.g., [12, 15]. However, a more common, but also more computationally expensive, approach is to use finite differences to approximate (2.7).

3. Using LBFGS for large-scale Kalman filtering. When the model size n is large, the Kalman filter is known to be prohibitively expensive to implement. This motivates several alternative approaches-most notably the 4D-Var method-used in large-scale applications such as numerical weather forecasting [8, 9, 11, 12, 17, 23, 26] and oceanography [4]. We instead focus our attention on the Kalman filter itself, using the limited memory BFGS (LBFGS) [22] iterative method for the required large-scale matrix storage and inversion within KF and EKF.

First, we give a general description of the LBFGS method for minimizing

(3.1) q(u) = 1/2(Au,u) - (b,u);

where A is an n x n symmetric positive definite matrix and b is an n x 1 vector. It is given by

The LBFGS method for quadratic minimization

v := 0;

[u.sub.0] := initial guess;

[B.sup.-1.sub.0] := initial inverse Hessian approximation; begin quasi-Newton iterations

[g.sub.v]:= [nabla]q([u.sub.v]) = A[u.sub.v] - b;

[v.sub.v] = [B.sup.-1.sub.v][g.sub.v];

[[tau].sub.v] = ([g.sub.v], [v.sub.v]}/([v.sub.v], [Av.sub.v]};

[u.sub.v+1] := [u.sub.v] - [[tau].sub.v] [v.sub.v];

[B.sup.-1.sub.v] := LBFGS approximation to [A.sup.-1]; end quasi-Newton iterations

In all of the examples considered in this paper, [B.sub.0] was taken to be the identity matrix. The limited memory formulations for [B.sup.-1.sub.v], and corresponding formulas for [B.sub.v], are found in

Appendix append1. The stopping criteria for the LBFGS iterations is discussed in Section 4.

Some insight into the convergence properties of the LBFGS method can be obtained by an appeal to its connection with the well-known conjugate gradient (CG) method, which described in detail in [20, 21, 22, Section 9.1]. In particular, CG can be formulated as a memoryless BFGS method. Moreover, in the presence of exact arithmetic, LBFGS and iterates from a certain preconditioned CG iteration are identical [21], and hence finite convergence is guaranteed. Thus it seems reasonable to suspect that LBFGS will have convergence properties similar to that of CG, which are well-known and have been extensively studied. In particular, the early convergence of CG iterates within the dominant subspaces corresponding to the largest eigenvalues of the coefficient matrix is likely shared by LBFGS iterates.

Next, we describe how LBFGS was used to make the Kalman filter more efficient. We make the reasonable assumption that multiplication by the evolution and observation matrices [M.sub.k] and [K.sub.k], and by the covariance matrices [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], is efficient, both in terms of storage and CPU time. Additional computational challenges arise for sufficiently large n due to the storage requirements for [C.sup.est.sub.k], which becomes a full matrix as the iterations proceed. The same is also true for [C.sup.p.sub.k]. However, given that

(3.2) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

storage issues are restricted to those for [C.sup.est.sub.k]; typically the matrix [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is assumed to be diagonal.

A low storage approximation of [C.sup.est.sub.k] can be obtained by applying the LBFGS algorithm to the problem of minimizing (3.1) with A = [C.sup.est.sub.k] and b = 0. The LBFGS matrix [B.sup.-1.sub.v] is then a low storage approximation of [([C.sup.est.sub.k]).sup.-1] and formulas for [C.sup.est.sub.k] from [3]-and also found in Appendix A.2-can be used.

Additionally, when m is sufficiently large the computation of ([MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII])) that is required in Step 2, (ii) of the Kalman filter iteration will be prohibitively expensive.

For the approximation of [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (3.1) and apply LBFGS to the problem of minimizing (3.1).

The LBFGS Kalman filter method can now be presented.

The LBFGS Kalman Filter (LBFGS-KF)

Step 0: Select initial guess ([C.sup.est.sub.k] and covariance [B.sub.#] = [C.sup.est.sub.k], and set k = 0.

Step 1: Compute the evolution model estimate and covariance:

(i) Compute ([x.sup.p.sub.k] = [M.sub.k]([e.sup.est.sub.k]

(ii) Define [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Step 2: Compute the Kalman filter estimate and covariance:

(i) Define A = [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] in (3.1) and compute the LBFGS approximations [B.sub.*] of [A.sup.-1] and[ .sup.*] of [A.sup.-1]b.

(ii) Compute the LBFGS-KF estimate [e.sup.est.sub.k+1] = [x.sup.p.sub.k] [C.sup.p.sub.k][K.sup.T.sub.k][u.sub.*];

(iii) Define A = [C.sup.p.sub.k] - [C.sup.p.sub.k] [K.sup.T.sub.k][B.sub.*][K.sub.k][C.sup.P.sub.k] ([approximately equal to] [C.sup.est.sub.k+1]) and b = 0 in (3.1) and compute the LBFGS approximation B# of [C.sup.est.sub.k+1] using (A.2).

Step 3: Update k := k + 1 and return to Step 1.

All operations with the [C.sup.est.sub.k] and [A.sup.-1] are done using the LBFGS formulas; see Appendix A. As a result, LBFGS-KF is much less memory and computationally intensive than KF making its use on large-scale problems more feasible. Specifically, the storage requirements for the LBFGS estimate of [C.sup.est.sub.k] are on the order of 2nl + 4n, where l is the number of stored vectors in LBFGS (typically 10-20), rather than n2 + 4n [22, Section 9.1], and the computational cost for both obtaining and using this estimate is order n. Furthermore, the inversion of the m x m matrix [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]is carried out in order m operations and its storage requirements are on the order of 2ml + 4m rather than [m.sup.2] + 4m [22].

The accuracy of the LBFGS covariance approximations is an important question. An analysis addressing this question in the similar variational setting is performed in [2]. We believe that the results of that analysis should be similar for LBFGS-KF. Thus, we choose not to repeat it here.

In the first example considered in the numerical experiments, LBFGS-KF and KF are compared and it is noted that LBFGS-KF is roughly 10 times faster, in terms of CPU time, than KF when applied to the same problem. Moreover, using our MATLAB implementation, LBFGS-KF can be used on significantly larger-scale problems.

As we have mentioned, in our implementations of KF and LBFGS-KF, the covariance matrices [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] are taken to be diagonal. This is not a necessary requirement. More structured covariances can be used, containing important prior information [17], however in order to maintain the computational efficiency and low storage requirements of LBFGS-KF, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] must be comparable to [M.sub.k], [B.sub.#] and [K.sub.k], [B.sub.*], respectively, in terms storage requirements and the computational cost required for their multiplication.

In the next section, we test the algorithm on two examples. The first is large-scale and linear, while the second is small-scale and nonlinear. The purpose of these experiments is to demonstrate that LBFGS-KF and LBFGS-EKF are effective algorithms. We leave their comparison with other state of the art methods for approximate Kalman filtering [5,8,9, 11, 12, 14, 17, 23, 26] for a later paper.

4. Numerical experiments. In this section, we present numerical results that justify the use of LBFGS-KF. In particular, we apply the method to two examples. The first is sufficiently large-scale that the use of LBFGS-KF is justified. In particular, we assume the following forced heat equation evolution model

(4.1) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where x is a function of u and v over the domain [OMEGA] = {(u, v) [less than or equal to] 0 < u, v [less than or equal to] 1} and [alpha] [greater than or equal to] 0. In our experiment, we will generate synthetic data using (4.1) with [alpha] > 0 and assume that the evolution model is given by (4.1) with [alpha] = 0, which gives a model bias. The problem can be made as large-scale as one wants via the choice of a sufficiently fine discretization of the domain [OMEGA].

However, the well-behaved nature of solutions of (4.1)-in particular the fact that its solutions tend to a steady state-makes further experiments with a different test case a necessity. For this reason, we also test our method on a second example, which contains chaotic solutions, and hence has unpredictable behavior. In particular, we consider the simple non-linear model introduced and analyzed in [18, 19] and which is given by

(4.2)[partial derivative][x.sup.i]/[partial derivative]t = ([x.sup.i+1] - [x.sup.i-2]) [x.sup.i-1]- [x.sup.i] + 8, i = 1,2,..., 40,

with periodic state space variables, i.e., [x.sup.-1] = [x.sup.n-1], [x.sup.0] = [x.sup.n] and [x.sup.n+1] = [x.sup.1], n = 40. Then (4.2) is a chaotic dynamical system (cf. [19]), which is desirable for testing purposes. As the model is computationally light and shares many characteristics with realistic atmospheric models (cf. [19]), it is commonly used for testing different data analysis schemes for weather forecasting.

4.1. An example with a large-scale linear evolution model. We perform our first experiments using model (4.1) using a uniform N x N computational grid and the standard finite difference discretization of both the time and spatial derivatives, which yields the following time stepping equation [x.sub.k+1] = [Mx.sub.k] + f, where M = I - [DELTA]tL. Here L is given by the standard finite difference discretization of the two-dimensional Laplacian operator with homogeneous Dirichlet boundary conditions, [DELTA]t is chosen to guarantee stability, and f is the constant vector determined by the evaluation of the forcing term in (4.1) at each of the points of the computational grid.

We define [K.sub.k] = K for all k in (2.2), where K is a matrix modeling an array of square sensors on the computational grid. Assuming that each sensor collects a weighted average of the state values in a 3 x 3 pixel region centered at every 8th pixel in both the x and y directions, K will have dimension (n/64) x n. We assume, further, that the weighted average in the 3 x 3 region is defined by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

In our first test, we generate synthetic data using the linear stochastic equations

[x.sub.k+1] = M[x.sub.k] = f + [N(0.5[[sigma].sub.ev]).sup.2]I),

[y.sub.k+1] = K[x.sub.k-1] = + [N(0.[(0.8[[sigma].sub.ev]).sup.2]I)

with [alpha] = 3/4 in (4.1) and where [[sigma].sup.2.sub.ev] and [[sigma].sup.2.sub.obs], are chosen so that the signal to noise ratios, defined by [parallel][x.sub.0][[parallel].sup.2]/[n.sup.2][[sigma].sup.2.sub.ev] and [parallel]K[x.sub.0][[parallel].sup.2]/[n.sup.2][[sigma].sup.2.sub.obs] respectively, are both 50. The initial condition used for the data generation was

[[[x.sub.o]].sub.ij] = exp[-[(([u.sub.i] - 1/2).sup.2] + [([v.sub.j] - 1/2).sup.2])],

where ([u.sub.i]; [v.sub.j]) is the ijth grid point.

In our implementation of KF, we used the biased models

(4.3) [x.sub.k+1] = M[x.sub.k] = + [N(0,[[sigma].sup.2.sub.ev]).sup.2]I),

(4.4) [y.sub.k+1] = K[x.sub.k+1] = + [N(0,[[sigma]sup.2.sub.obs]I),

with initial conditions [x.sub.0] = 0 and [C.sup.est.sub.0] = 0.001I in Step 0 of the filter. We compare the results obtained with the LBFGS-KF and KF, where N = [2.sup.j] with j taken to be the largest positive integer so that memory issues do not arise in the MATLAB implementation for the standard KF. For the computer on which the simulations were done (a laptop with 2G RAM memory and a 1.8 GHz Core 2 Duo processor) the largest such j was 5, making N = 32 and n = 1024. We note that in our implementation of the LBFGS method within LBFGS-KF, we have chosen to take only 10 LBFGS iterations with 9 saved vectors. These choices may seem crude at first, however, more stringent stopping tolerances and/or a larger number of stored vectors did not appreciably affect the results for the examples that we considered.

The purpose of this test is to show that the results obtained with LBFGS-KF are comparable results with those obtained with KF. To do this, we present a plot in Figure 4.1 of the relative error vector, which has kth component

[[relative_error].sub.k]]: = [parallel][e.sup.est.sub.k] -[x.sub.k][parallel]/[parallel][x.sub.k][parallel],

for both the LBFGS Kalman Filter and for the standard Kalman Filter. We see that results obtained using the two approaches yield similar, though not identical, relative error curves. Both curves eventually begin to increase once the forcing term, which is not used in the state space model in KF, has a prominent effect on the data; in early iterations, it is overwhelmed by the diffused initial temperature. We also mention that in the large number of test runs that we did using this large-scale model, our implementation of the LBFGS-KF was on average about 10 times faster than was the standard KF.

Additionally, in Figure 4.2, we present the filter estimates obtained from both KF and LBFGS-KF together with the true state values at time points 35 and 70. Note that in the early iterations of the filter, represented by time point 35, the filter does not detect the source because it is overwhelmed by the initial temperature and is not contained in the model (4.3), (4.4). However, the source is detected once the initial temperature has sufficiently dissipated.

For a thorough comparison, we perform the same test using values for [[sigma].sup.2.sub.ev] and [[sigma].sup.2.sub.obs] that yield signal-to-noise ratios of 10. The relative error curves in Figure 4.3 result. Interestingly, LBFGS-KF provides better results at the beginning of the filtering period than does KF. This can be explained, we believe, by the fact that a regularization of sorts is implicitly implemented via the use of a truncated LBFGS algorithm.

Finally, we choose [[sigma].sup.2.sub.ev] and [[sigma].sup.2.sub.obs] as in the original experiment, but take [alpha] = 2, which has the effect of making the state space model that is used within LBFGS-KF and KF less accurate. When this is done, we obtain the solution curves appearing in Figure 4.4. Thus it seems that as the underlying evolution model becomes less accurate and the noise level remains moderately low KF provides better results

[FIGURE 4.1 OMITTED]

In order to show that satisfactory results can also be obtained for much larger scale problems, we take j = 8 which gives N = 256 and n = 65536. We take all other parameter values to be those of the original experiment. Note, however, that the stability condition of the time stepping scheme requires a much smaller time step for this problem. A relative error plot similar to those in the previous example is given in Figure 4.5. We do not include an error curve for the Kalman filter because memory issues prevent its implementation on our computer for either N = 128 (n = 16384) or N = 256 (n = 65536).

The previous large-scale example remains orders of magnitude smaller than the typical size of systems considered in practical weather models. We stopped at N = 256 because our experiments were performed on a laptop that could not handle a larger-scale problem. However, the discussion of computational cost and storage in the paragraph following the description of the LBFGS-KF algorithm suggests that it scales well with problem size. Thus the use of LBFGS-KF on much larger-scale problems should be feasible. Efficiency can be further improved if several time steps are allowed in the forward model for each Kalman filter iteration, much as is done in 4D-Var implementations. In addition, to the degree that LBFGS is parallelizable, LBFGS-KF will also be parallelizable.

4.2. An example with a small-scale, nonlinear evolution model. In our next example, we apply EKF and LBFGS-EKF to the problem of estimating the state variables from data generated using the nonlinear, chaotic evolution model (4.2). To generate the data, a time integration of the model was first performed using a fourth order Runge-Kutta (RK4) method with time-step [DELTA]t = 0.025. Analysis in [19] suggests that when (4.2) is used as a test example for weather forecasting data assimilation algorithms, the characteristic time scale is such that the above [DELTA]t corresponds to 3 hours, which we will use in what follows. It is also noted in [19] that for [DELTA]t [less than or equal to] 0.5, the RK4 method is stable. The "true data" was generated by taking 42920 time steps of the RK4 method, which corresponds to 5365 days. The initial state at the beginning of the data generation was [x.sup.20] = 8 + 0.008 and [x.sub.i] = 8 for all i [not equal to] 20.

[FIGURE 4.2 OMITTED]

The observed data is then computed using this true data. In particular, after a 365 day long initial period, the true data is observed at every other time step and at the last 3 grid points in each set of 5; that is, the observation matrix is m x n with nonzero entries

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII],

The observation error is simulated using the Gaussian random vector N(0, (0.15 [[[sigma].sub.clim]).sup.2]I) where [[sigma].sub.clim] is a standard deviation of the model state used in climatological simulations, [[sigma].sub.clim] := 3.6414723. The data generation codes were written in MATLAB and were transcribed by us from the scilab codes written by the author of [15].

[FIGURE 4.3 OMITTED]

[FIGURE 4.4 OMITTED]

[FIGURE 4.5 OMITTED]

In our application of EKF and LFBGS-EKF, we assume the coupled stochastic system

(4.5) [x.sub.k+1] = M[x.sub.k] + [N(0, (0, (0.05[[sigma].sub.clim]).sup.2]I),

(4.6) [y.sub.k+1] = K[x.sub.k+1] + [N(0, (0, (0.05[[sigma].sub.clim]).sup.2]I)

where M ([x.sub.k]) is obtained by taking two steps of the RK4 method applied to (4.2) with initial condition [x.sub.k] with time-step 0.025. We note that if the noise term is removed from (4.5) and the above initial condition is used, our data generation scheme results.

Due to the fact that M is a nonlinear function, EKF must be used; see (2.5) and (2.6). Since K := K in (2.6) is linear, [K.sub.k] = K for all k in (2.7). However, a linearization of the nonlinear evolution function M is required. Fortunately, the computation of [M.sub.k] in (2.7) is performed by a routine in one of the scilab codes mentioned above and that we have adapted for our use in MATLAB.

The initial condition used in implementation of both the EKF and LBFGS-EKF is defined by [[[x.sub.t0]].sub.i] = [[x.sup.true.sub.t0]].sub.i] + N(0, (0.3 [[[sigma].sub.clim]).sup.2]) for all i, and the initial covariance was taken to be [C.sup.est.sub.0] = (0.13 [[sigma].sub.clim]).sup.2]I. In our implementation of the LBFGS method within LBFGS-EKF, we computed 10 iterations with 9 saved vectors.

In order to analyze the accuracy of the state estimates [x.sup.est.sub.k] obtained by both EKF and LBFGS-EKF we plot the vector with components

(4.7) [[rms].sub.k] = [square root of 1/40[parallel][x.sup.est.sub.k] - [x.sup.true.sub.k][parallel].sup.2]]

in Figure 4.6. We can see that the two methods yield comparable results.

[FIGURE 4.6 OMITTED]

In order to compare the forecasting abilities of the two approaches, we compute the following forecast statistics at every 8th observation. Take [member of] I := {8i | i = 1, 2,..., 100} and define

(4.8) [[forcast.error.sub.j].sub.i] = 1/40 [parallel][M.sub.4i]([x.sup.est.sub.j]) - [x.sup.true.sub.j+4i][parallel].sup.2]], i = 1,..., 20,

where [M.sub.n] denotes a forward integration of the model by n time steps with the RK4 method. Thus this vector gives a measure of forecast accuracy given by the respective filter estimate up to 80 time steps, or 10 days out. This allows us to define the forecast skill vector

(4.9) [[forecast_skill].sub.i] = 1//[[sigma].sub.clim][square root of 1/00[summation (j[member of]I][forecast_error.sub.j].sub.i]], i = 1,..., 20,

which is plotted in Figure 4.7. The results show that the forecasting skill of the two methods is very similar, which suggests that on the whole, the quality of the LBFGS-EKF estimates is as high as those obtained using EKF. Figure 4.7 also illustrates the fact that the Lorenz 95 model (4.2) is truly chaotic.

In the test cases considered here, a linear or linearized model matrix [M.sub.k] has been available. This is not true in important examples such as in numerical weather forecasting, where, on the other hand, a tangent linear code [14] is available that provides a means of computing the matrix vector product [M.sub.k]x.

5. Conclusions. The standard implementations of KF and EKF become exceedingly time and memory intensive as the dimension of the underlying state space increases. Several variants of KF and EKF have been proposed to reduce the dimension of the system, thus making implementation in high dimensions possible. The Reduced Rank Kalman Filter or Reduced Order extended Kalman filter (see, e.g., [4,7, 30]) project the dynamical state vector of the model onto a lower dimensional subspace. The success of this approach depends on a judicious choice of the reduction operator. Moreover, since the reduction operator is typically fixed in time, they can suffer from "covariance leaks" [9]. A typical cause of this is that nonlinear systems do not generally leave any fixed linear subspace invariant.

[FIGURE 4.7 OMITTED]

In this paper, we propose the use of the limited memory BFGS (LBFGS) minimization method in order to circumvent the computational complexity and memory issues of standard KF and EKF. In particular, we replace the n x n, where n is the dimension of the state space, covariance matrices within KF and EKF with low storage approximations obtain using LBFGS. The large-scale matrix inversions required inKF and EKF implementations are also approximated using LBFGS. The resulting methods are denoted LBFGS-KF and LBFGSEKF, respectively.

In order to test these methods, we consider two test cases: large-scale linear and small scale nonlinear. LBFGS-KF is applied in the large-scale linear case and is shown to be effective. In fact, our method exceeds the speed of standard KF by an order of magnitude, and yields comparable results when both methods can be applied. Furthermore, it can be used on much larger scale problems. In the nonlinear, small scale case, LBFGS-EKF is implemented and is also shown to give results that are comparable to those obtained using standard EKF. We believe that these results suggest that our approach deserves further consideration.

The symmetric rank one (SR1) quasi-Newton method for minimizing (3.1) could be another attractive method for use within KF and EKF, since it also yields estimates of both the minimizer and inverse Hessian. The main drawback of using SR1, however, is that the inverse Hessian approximations are not guaranteed to be positive definite.

Appendix A. In this appendix, we give a general description of the LBFGS method for minimizing

q(u) = 1/2(Au,u) - (b, u);

where A is an n x n symmetric positive definite matrix and b is an n x 1 vector. It is given by.

The LBFGS method for quadratic minimization

v := 0;

[u.sub.0] := initial guess;

[B.sup.-1.sub.0] := initial inverse Hessian approximation; begin quasi-Newton iterations

[g.sub.v] := [nabla]q([u.sub.v]) = A[u.sub.v] - b;

[v.sub.v] = [B.sup.-1.sub.v][g.sub.v];

[[tau].sub.v] = ([g.sub.v], [v.sub.v]}/([v.sub.v], A[v.sub.v]};

[u.sub.v+1] := [u.sub.v] - [[tau].sub.v][v.sub.v];

[B.sup.-1.sub.v] := LBFGS approximation to [A.sup.-1]; end quasi-Newton iterations

A.1. The limited memory approximation for [A.sup.-1]. The BFGS matrix [B.sup.-1] is computed using recursion

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII],

However, for large-scale problems the storage of the full matrix [B.sup.-1.sub.v] is infeasible, which motivates the limited storage version of the algorithm. At iteration v, suppose that the j vector pairs [{[s.sub.i], [d.sub.i]}.sup.v-1.sub.i=-1-j] are stored. Then we the LBFGS approximation of the inverse Hessian is given by

(A.1) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Assuming exact arithmetic and that j, we have that [u.sub.v] converges to the unique minimizer of q in at most n iterations, and if n iterations are performed [B.sup.-1.sub.n+1] = [A.sup.-1] [22]. In the implementation in this paper, however, j << n and LBFGS iterations are truncated before convergence is obtained.

It is proven in [21] that for quadratic minimization problems and exact line searches LBFGS is equivalent to preconditioned conjugate gradient with fixed preconditioner B0. Thus its convergence rate is given by [22]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

A.2. A low storage approximation of A. The required formulas are given in [3] and take the following form. Let

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

We note that when [[xi].sub.v] = 1 for all v in (A.2), we have an exact equality between [B.sub.v] in (A.2) and (A.1). However, we have found that a more accurate Hessian approximation is obtained if, following [22], we use the scaling [[xi].sub.v] = [d.sup.T.sub.v-1] [d.sub.v-1][s.sup.T.sub.v-1][d.sub.v] instead.

We note that the middle matrix in (A.2) has size 2j x 2j, which is of reasonable size provided j is not too large, and its inversion can be carried out efficiently using a Cholesky factorization that exploits the structure of the matrix; see [3] for details.

Acknowledgments. The authors would like to thank the referees and ETNA editors; the paper is better because of their efforts. This work was supported by the Academy of Finland (application number 213476, Finnish programme for Centres of Excellence in research 2006-2011). We are thankful for M. Leutbecher for providing us with the Scilab version of the Lorenz95 code, that served as starting point for the Matlab model implementation. The second author would like to thank both the University of Montana and the University of Helsinki for their support during his stay in Finland during the 2006-07 academic year, where this work was undertaken.

REFERENCES

[1] H. Auvinen, H. Haario, and T. Kauranne, Optimal approximation of Kalman filtering with a temporally local 4D-Var in operational weather forecasting, in Proceedings of the 11th ECMWF Workshop on Use of High-Performance Computing in Meteorology, W. Zwiefelhofer and G. Mozdzynski, eds., World Scientific, London, 2005.

[2] H. Auvinen, J. M. Bardsley, H. Haario, and T. Kauranne, The variational Kalman filter and an efficient implementation using limited memory BFGS, Internal J. Numer. Methods Fluids, to appear.

[3] R. H. Byrd, J. Nocedal, and R. B. Schnabel, Representations of quasi-Newton matrices and their use in limited memory methods, Math. Program., 63 (1994), pp. 129-156.

[4] M. A. Cane, R. N. Miller, B. Tang, E. C. Hackert, and A. J. Busalacchi, Mapping tropical Pacific sea level: data assimilation via reduced state Kalman filter, J. Geophys. Res., 101 (1996), pp. 599-617.

[5] Y. Cao, J. Zhu, I. M. Navon, and Z. Luo, A reduced order approach to four-dimensional variational data assimilation using proper orthogonal decomposition, Internat. J. Numer. Methods Fluids, 53 (2007), pp. 1571-1583.

[6] R. Daley, Atmospheric Data Analysis, Cambridge University Press, Cambridge, UK, 1991.

[7] D. P. Dee, Simplification of the Kalman filter for meteorological data assimilation Q. J. R. Meteorol. Soc., 117(1990), pp. 365-384.

[8] M. Fisher, Minimization algorithms for variational data assimilation, in Recent Developments in Numerical Methods for Atmospheric Modelling, Reading, UK, September 7-11, 2008, European Centre for Medium-Range Weather Forecasts, Reading, UK, 1998, pp. 364-385.

[9] M. Fisher and E. Andersson, Developments in 4D-var and Kalman filtering, European Centre for Medium-Range Weather Forecasts Technical Memorandum 347, Reading, UK, 2001.

[10] M. Fisher and P. Courtier, Estimating the covariance matrices of analysis and forecast error in variational data assimilation, European Centre for Medium-Range Weather Forecasts Technical Memorandum 220, Reading, UK, 1995.

[11] M. Fisher, M. Leutbecher, and G. Kelley, On the equivalence between Kalman smoothing and weak-constraint four-dimensional variational data assimilation, Q. J. R. Meteorol. Soc., 131 (2005), pp. 3235-3246.

[12] M. Fisher, J. Nocedal, Y. TrEmolet, and S. J. Wright, Data assimilation in weather forecasting: a case studyin PDE-constrained optimization, Optim. Eng., 10 (2009), pp. 1389-4420.

[13] I. Y. Gejadze, F.-X. Le Dimet, and V. Shutyaev, On analysis error covariances in variational data assimilation, SIAM J. Sci. Comput., 30 (2008), pp. 1847-1874.

[14] F. -X. LeDimet and O. Talagrand, Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects, Tellus Ser. A, 38 (1986), pp. 97-110.

[15] M. Leutbecher, A data assimilation tutorial based on the Lorenz-95 system, European Centre for Medium-Range Weather Forecasts Web Tutorial, Reading, UK. Available at www.ecmwf.int/newsevents/training/meteorological_presentations/pdf/DA/ ToyModel_notes.pdf.

[16] Z. Li and M. Navon, Optimality ofvariational data assimilation and its relationship with the Kalman filter and smoother, Q. J. R. Meteorol. Soc., 127 (2008), pp. 661-683.

[17] A. C. Lorenc, Modelling of error covariances by 4D-Var data assimilation, Q. J. R. Meteorol. Soc., 129 (2003), pp. 3167-3182.

[18] E. N. Lorenz, Predictability: A problem partly solved, in Predictability of Weather and Climate, Tim Palmer and Renate Hagedorn, eds., Cambridge University Press, Cambridge, UK, 2006, pp. 40-58. Originally presented in a 1996 European Centre for Medium-Range Weather Forecasts workshop.

[19] E. N. Lorenz and K. A. Emanuel, Optimal sites for supplementary weather observations: simulation with a small model, J. Atmospheric Sci., (1998), pp. 399-414.

[20] L. Nazareth, A relationship between the BFGS and conjugate gradient algorithms and it implications for new algorithms, SIAM J. Numer. Anal., 16 (1979), pp. 794-800.

[21] J. Nocedal, Updating quasi-Newton matrices with limited storage, Math. Comp., 35 (1980), pp. 773-782.

[22] J. Nocedal and S. Wright, Numerical Optimization, Springer, New York, 1999.

[23] F. Rabier, H. Jarvinen, E. Klinker, J. F. Mahfouf, and A. Simmons, The ECMWF operational implementation offour-dimensional variational assimilation. Part I: experimental results with simplified-physics, Q. J. R. Meteorol. Soc., 126 (2000), pp. 1143-1170.

[24] C. D. Rodgers, Inverse Methods for Atmospheric Sounding: Theory and Practice, World Scientific, London, 2000.

[25] X. Tian, Z. Xie, and A. Dai, An ensemble-based explicit four-dimensional variational assimilation method, J. Geophys. Res., 113 (2008), pp. D21124 (1-13).

[26] J. Tshimanga, S. Gratton, A. T. Weaver, and A. Sartenaer, Limited-memory preconditioners, with application to incremental four-dimensional variational data assimilation, Q. J. R. Meteorol. Soc., 134 (2008), pp. 751-769.

[27] F. Veerse, Variable-storage quasi-Newton operators as inverse forecast/analysis error covariance matrices in variational data assimilation, INRIA Report 3685, April 26, 1999.

[28] F. Veerse, Variable-storage quasi-Newton operators for modeling error covariances, in Proceedings of the Third WMO International Symposium on Assimilation of Observations in Meteorology and Oceanography, 1999, Quebec City, Canada, World Meteorological Organization, Geneva.

[29] F. Veerse, D. Auroux, and M. Fisher, Limited-memory BFGS diagonal predonditioners for a data assimilation problem in meteorology, Optim. Eng., 1 (2000), pp. 323-339.

[30] A. Voutilainen, T. Pyhalahti, K. Kallio, J. Pulliainen, H. Haario, and J. Kaipio, A filtering approach for estimating lake water quality from remote sensing data, Int. J. Appl. Earth Observation & Geoinformation, 9 (2007), pp. 50-64.

[31] W. Yang, I. M. Navon, and P. Courtier, A new Hessian preconditioning method applied to variational data assimilation experiments using NASA general circulation models, Monthly Weather Rev., 124 (1996), pp. 1000-1017.

H. AUVINEN ([dagger]) J. M. BARDSLEY ([double dagger]), H. HAARIO ([dagger]) AND T. KAURANNE ([dagger])

([dagger]) Department of Mathematics and Physics, Lappeenranta University of Technology, Lappeenranta, Finland ({harri.auvinen,heikki.haario, tuomo.kauranne}@lut.fi).

([double dagger]) Department of Mathematical Sciences, University of Montana, Missoula, Montana 59812, USA (bardsleyj@mso.umt.edu).

* Received April 11, 2008. Accepted for publication August 29, 2009. Published online on November 30, 2009. Recommended by P. Van Dooren. This work was supported by a faculty exchange grant from the University of Montana.

Printer friendly Cite/link Email Feedback | |

Title Annotation: | Broyden-Fletcher-Goldfarb-Shanno |
---|---|

Author: | Auvinen, H.; Bardsley, J.M.; Haario, H.; Kauranne, T. |

Publication: | Electronic Transactions on Numerical Analysis |

Article Type: | Report |

Geographic Code: | 1USA |

Date: | Jan 1, 2009 |

Words: | 7216 |

Previous Article: | Gaussian direct quadrature methods for double delay Volterra integral equations. |

Next Article: | An efficient generalization of the rush-Larsen method for solving electro-physiology membrane equations. |

Topics: |