In recent years, several ensemble-based Kalman filter algorithms have been developed that have been classified as ensemble square root Kalman filters. Parallel to this development, the singular “evolutive” interpolated Kalman (SEIK) filter has been introduced and applied in several studies. Some publications note that the SEIK filter is an ensemble Kalman filter or even an ensemble square root Kalman filter. This study examines the relation of the SEIK filter to ensemble square root filters in detail. It shows that the SEIK filter is indeed an ensemble square root Kalman filter. Furthermore, a variant of the SEIK filter, the error subspace transform Kalman filter (ESTKF), is presented that results in identical ensemble transformations to those of the ensemble transform Kalman filter (ETKF), while having a slightly lower computational cost. Numerical experiments are conducted to compare the performance of three filters (SEIK, ETKF, and ESTKF) using deterministic and random ensemble transformations. The results show better performance for the ETKF and ESTKF methods over the SEIK filter as long as this filter is not applied with a symmetric square root. The findings unify the separate developments that have been performed for the SEIK filter and the other ensemble square root Kalman filters.
The original ensemble Kalman filter (EnKF; Evensen 1994) has been developed with the aim to enable the application of sequential data assimilation algorithms based on the Kalman filter with large-scale numerical models. Burgers et al. (1998) and Houtekamer and Mitchell (1998) clarified that the EnKF requires an ensemble of perturbed observations for statistical consistency. The EnKF represents the state estimate by the mean of an ensemble of model state realizations, while the ensemble covariance matrix represents the corresponding error covariance matrix. The prediction of the error covariance matrix is computed by propagating each model state of the ensemble with the full, usually nonlinear, numerical model.
Alternative filter algorithms have been developed that perform the analysis without perturbed observations. These filters use an explicit transformation of the state ensemble. Among these developments are the ensemble transform Kalman filter (ETKF; Bishop et al. 2001), the ensemble adjustment Kalman filter (EAKF; Anderson 2001), and the ensemble square root Kalman filter with sequential processing of observations (EnSRF; Whitaker and Hamill 2002). These filters also have been reviewed by Tippett et al. (2003) in a uniform way as ensemble square root Kalman filters. Another ensemble square root Kalman filter has been derived by Evensen (2004).
The ensemble-based singular “evolutive” interpolated Kalman (SEIK) filter has been introduced by Pham et al. (1998) a few years before the introduction of the ensemble square root Kalman filters. The behavior of SEIK filter for nonlinear models was examined by Pham (2001). Comparison studies between the SEIK filter and the EnKF (Brusdal et al. 2003; Nerger et al. 2005a) argue that the SEIK filter can be more efficient than the EnKF because a smaller ensemble could be used to achieve comparable estimation errors. In addition, the computations used in the SEIK filter are much less costly than those of the EnKF (Nerger et al. 2007).
Overall, the developments in the SEIK filter and the ensemble square root Kalman filters have been independent. In publications considering ensemble square root filters, the SEIK filter is only occasionally mentioned. For example, Sakov and Oke (2008) note that the SEIK and SEEK filters “essentially represent another flavor” of the ensemble square root filter. Similarly, publications using the SEIK filter, describe it as an efficient alternative to the EnKF (e.g., Triantafyllou et al. 2003; Nerger et al. 2005a). Thus, while there are indications that the SEIK filter is an ensemble square root filter, there is yet no clear classification of the SEIK filter or an identification of the square root used in this algorithm.
The aim of this work is to examine the relation of the SEIK filter to the ensemble square root Kalman filters in detail. For this task, the ETKF and the SEIK filter will be reviewed in section 2. In section 3 it is shown that the SEIK filter is an ensemble square root filter and its relation to the ETKF is discussed. A variant of the SEIK filter that results in identical ensemble transformations to those of the ETKF, which we term the error subspace transform Kalman filter (ESTKF), is derived in section 4. The computational cost of the filters as well as a possible reduction of the cost of the ETKF is discussed in section 5. Numerical experiments are performed in section 6 to compare the filter behavior for different variants of the ensemble transformation matrix.
2. Filter algorithms: ETKF and SEIK
In this section, the mathematical formulations of the ETKF and the SEIK filter are reviewed and the square root in the ETKF is identified in analogy to Tippett et al. (2003). Only the global analysis formulation is considered. A localization (see Nerger et al. 2006; Hunt et al. 2007) can be formulated in an identical way for both filters.
The ETKF and the SEIK filter are ensemble-based Kalman filters. The state of a physical system, like the ocean or atmosphere, is estimated at time tk by the state vector xk of size n and the corresponding error covariance matrix . An ensemble of m vectors x(α), α = 1, … , m, of model state realizations represents these quantities. The state estimate is given by the ensemble mean:
With the ensemble matrix
is given as the ensemble covariance matrix:
where with is the matrix of ensemble perturbations.
A forecast is computed by integrating the state ensemble using the numerical model until observations become available. The observations are available in form of the vector of size p. The model state is related to the observations by where is the observation operator, which is assumed to be linear. The vector of observation errors, εk, is assumed to be a white Gaussian distributed random process with covariance matrix .
The analysis equations of the ETKF and the SEIK filter are discussed separately below. As all operations are performed at the same time tk, the time index k is omitted.
a. Analysis step of the ETKF
The computations performed in the ETKF are based on a square root of the state covariance matrix given by the ensemble perturbations . The analysis state covariance matrix can be written as a transformation of the forecast ensemble perturbations as
Here, is an m × m matrix defined by
The is frequently denoted as the “transform matrix.” The factor γ is used to inflate the forecast covariance matrix to stabilize the filter performance.
The state estimate is updated according to
with the weight vector
The square root of the forecast state covariance matrix is given by the perturbation matrix up to the scaling by (m − 1)−1. To obtain the square root of the analysis state covariance matrix, is transformed as
The weight matrix is computed from the square root with as
Here, Λ is an arbitrary orthogonal matrix of size m × m or the identity. To preserve the ensemble mean, the vector (1, … , 1)T has to be an eigenvector of Λ.
When the ETKF was introduced by Bishop et al. (2001), the form of the square root was not further specified. Studies about the properties of the ensemble transformation in different square root filters (e.g., Wang et al. 2004; Sakov and Oke 2008) have shown that a symmetric matrix ensures that the ensemble mean is preserved during the ensemble transformation. The use of the symmetric square root:
has been proposed also for the localized version of the ETKF (LETKF; Hunt et al. 2007). Equation (10) can be obtained from the singular value decomposition (SVD) . The use of matrix from Eq. (10) provides a minimum transformation of the ensemble because the distance of the square root from the identity matrix is minimized in the Frobenius norm (see Yang et al. 2009).
with . This formulation leads directly to the analysis ensemble, without explicitly updating the state estimate by Eq. (6).
b. Analysis step of the SEIK filter
The SEIK filter has been introduced by Pham et al. (1998) and was described in more detail by Pham (2001). This review follows Nerger et al. (2006). The original separation of the analysis step into the state update (“analysis”) and ensemble transformation (“resampling”) is followed here. The SEIK filter is then explicitly reformulated as an ensemble square root filter analogously to the ETKF in section 3. Quantities that are similar but not identical to those of the ETKF are marked using a tilde. It is assumed that the forecast ensemble is identical to that used in the ETKF.
The computations of the analysis step update the state estimate and implicitly update the state covariance matrix from the forecast to the analysis matrix.
In the SEIK filter, the forecast covariance matrix is treated in terms of the forecast state ensemble by
Here, is an m × (m − 1) matrix with full rank and zero column sums. Previous studies have always defined matrix as
where 0 represents the matrix whose elements are equal to zero and is the identity. The elements of the matrix 1 are equal to one. Matrix implicitly subtracts the ensemble mean when the matrix is computed. In addition, removes the last column of , thus is an n × (m − 1) matrix that holds the first m − 1 ensemble perturbations.
The analysis update of the state estimate is given as a combination of the columns of the matrix by
Here, the vector of size m − 1 is given by
and the transform matrix of size (m − 1) × (m − 1) is defined by
In the SEIK filter, with is referred to as the “forgetting factor.” It is the inverse of the inflation factor γ used in Eq. (5) of the ETKF. The analysis covariance matrix is given in factorized form by
but does not need to be explicitly computed.
For efficiency, the term is typically computed as . Thus, operates on the p × m matrix , while operates on each ensemble state.
After the analysis step, the resampling of the ensemble is performed. Here, the forecast ensemble is transformed such that it represents and . The transformation is performed according to
In previous studies, the SEIK filter was always described to use a Cholesky decomposition of the matrix to obtain . However, other forms of the square root, like the symmetric square root used in the ETKF, could be chosen. Section 6 will test the influence of the chosen square root on the performance of the filter. The matrix Ω is an m × (m − 1) matrix whose columns are orthonormal and orthogonal to the vector (1, … , 1)T. Traditionally, Ω is described to be a random matrix with these properties. However, using a deterministic Ω is also valid. The procedure to generate a random Ω (Pham 2001; Hoteit 2001) and a procedure for generating a deterministic variant are provided in the appendix.
For efficiency, the matrix can be replaced by [Eq. (13)]. Then, the matrix can be applied from the left to smaller matrices like the weight vector or the matrix .
The original formulation of the SEIK filter used the normalization m−1 for the matrix instead of using the sample covariance matrix that is normalized by (m − 1)−1. For consistency with other ensemble-based Kalman filters, Nerger and Gregg (2007) introduced the use of the sample covariance matrix in SEIK, which is also used here. In the SEIK filter, the ensemble is generated to be consistent with the normalization of . Hence, the normalization acts only as a scaling factor that influences Eqs. (3) and (20) as well as the definition of in Eq. (14).
3. SEIK as an ensemble square root filter
To identify the SEIK filter as an ensemble square root filter, the analysis and resampling steps of SEIK are combined as a transformation of the square root of . Equation (20) can be written as
Equation (23) performs a transformation of the matrix analogous to the ensemble transformation of the ETKF [Eq. (11)]. Matrix is the square root of the covariance matrix used in the SEIK filter. With this, the SEIK filter is clearly an ensemble square root filter.
It is particular for the SEIK filter that the matrix has only m − 1 columns, while other filters use a square root with m columns. Using m − 1 columns is possible because the rank of is at most m − 1. The SEIK filter utilizes this property by accounting for the fact that the sum of each row of the perturbation matrix is zero. Thus, while the columns of are linearly dependent, the columns of are linearly independent if the rank of is m − 1. In this case, they build a basis of the error subspace estimated by the ensemble of model states (for a detailed discussion of the error subspace, see Nerger et al. 2005a). In contrast, can be regarded as a transformation from its m-dimensional column space to the error subspace of dimension m − 1 (see Hunt et al. 2007).
While the equations of the SEIK filter are very similar to those of the ETKF this does not automatically imply that their state and error estimates are identical, in particular because the analyses use matrices of different size. However, if the same forecast ensembles are used in the ETKF and the SEIK filter, the analysis state and the analysis state covariance matrix will be identical. This identity is due to the fact that the analysis formulations of both methods refer to the same error subspace to compute the optimal combination of ensemble perturbations. A basis of this space is given by . It is used directly by the SEIK filter. In contrast, the ETKF utilizes the ensemble representation of the error subspace given by . Nonetheless, the matrices [Eq. (5)] and [Eq. (18)] both describe the same quantity—an error covariance matrix—in the same space represented by either or . Therefore, the optimization computed in the analysis steps results in the same state and error estimates.
While the identity of xa and for both filters can be established by the argumentation above, the ensembles that represent these quantities are only unique up to a unitary matrix [i.e., ; see, e.g., Livings et al. (2008)]. For example, this is the case when random rotations are used to generate Ω or Λ. However, for deterministic transformations and in the use of the symmetric square root of , the experiments discussed in section 6 indicate that the differences between the transformation matrices of the SEIK filter and the ETKF are very small with differences in the matrix entries below 2%.
4. Identical transformations in SEIK and ETKF
The ensemble transformation in the square root formulation of SEIK, which was discussed in section 3, generally exhibits very small deviations from the transformation performed by the ETKF. As the transformation in the ETKF has been described to be the minimum transformation, it should be desirable to obtain the same transformation with the SEIK filter. This goal is achieved by a modification of the SEIK filter that is described in this section.
The modification of the SEIK filter is motivated by the properties of the matrix Ω. In general, Ω is an m × (m − 1) matrix that regenerates m ensemble perturbations in combination with an ensemble transformation matrix of size (m − 1) × (m − 1). For a deterministic ensemble transformation, a deterministic form can be used whose elements are defined by
Geometrically, is the Householder matrix associated with the vector m−½(1, … , 1)T (see the appendix). Thus, projects vectors in the ensemble space spanned by onto the error subspace spanned by . Like , has a full rank and zero column sums. In addition, the columns of are orthonormal, which is not the case for . Using , one can replace Eqs. (12)–(14) by
Now, matrix from Eq. (18) is computed as
Finally, the ensemble transformation [Eq. (20)] becomes
where is the square root of . Here is the projection of from the error space onto the ensemble space. If the symmetric square root is used to compute , the projected transformation matrix is identical to the matrix used in the ETKF. In case of random ensemble transformations, only the rightmost in Eq. (29) is replaced by the random matrix Ω, while is used at all other places.
This reformulation of the SEIK filter is consistent with its original motivation to compute the ensemble transformation matrix in the error space and to project the required matrices onto this space and finally back onto the ensemble space. The choice of is arbitrary as long as its column sums are zero and the matrix is of full rank. However, only the application of results in consistent projections, because it is symmetrically applied in the computation of as well as in the ensemble transformation [Eq. (29)]. Because the ensemble transformation is performed in the error subspace, the new filter variant is referred to as the ESTKF. The main difference between the SEIK filter and the ESTKF is that the application of in Eq. (13) subtracts the ensemble mean and drops the last ensemble member. The resulting matrix actually depends on the order of the ensemble members in the ensemble matrix , which is arbitrary. In contrast, matrix defined by Eq. (26) will be independent of the order of the ensemble members. This is evident from the action of when computing not only subtracts the ensemble mean, but also subtracts the value of the last column of divided by from each column. The columns of are then normalized by an additional division by . These operations ensure that the value of the last column of is implicitly contained in matrix .
The use of instead of does not change the computational cost of the filter. The matrix needs also to be initialized in the previous formulation of the SEIK filter. In addition, the multiplication of a matrix by has the same cost as the multiplication by .
5. Comparison of the computational costs and algorithmic enhancement of the ETKF
The computational cost of the SEIK filter is very similar to that of the ETKF. The leading costs of both filters are summarized in Table 1. The leading computational cost of both filter algorithms scales in the same way. However, the cost of the SEIK filter is slightly lower because of the use of matrix with m − 1 columns instead of with m columns.
One second-order term that does not appear explicitly in Table 1 is the computation of in the ETKF with a cost of O(nm). The SEIK filter applies the matrix to and to [Eq. (21)]. In the ESTKF, the matrix is applied analogously. These operations have a cost of O[p(m − 1) + m(m − 1)2]. In the typical situation, where the state dimension n is much larger than the observation dimension p and the ensemble size m is smaller than p, this alternative will be computationally less costly.
The ETKF can be modified to use an analog to matrix . The computation of the perturbation matrix can be formulated as
where the m × m matrix is defined by
Now, the equations of the ETKF that involve can be reformulated. Equation (5) becomes
and Eq. (7) is written as
Further, the transformation Eq. (11) becomes
As in the SEIK filter, this formulation avoids the explicit computation and storage of the ensemble perturbation matrix . Instead, the matrix is applied to of size p × m and to the sum of the weight matrices in Eq. (34) of size m × m. This changes the computational cost to O(pm + m3) instead of O(nm) for the direct computation of . This formulation can also be applied with domain localization, but here should be computed globally, before performing the local analyses.
6. Numerical experiments
a. Experimental setup
In this section, the behavior of the ETKF will be compared with the explicit square root formulation of the SEIK filter using the symmetric square root introduced in section 3 (referred to as SEIK-sqrt) and with the ESTKF. In addition, the original SEIK filter with a square root based on Cholesky decomposition from section 2b is applied (referred to as SEIK-orig). To compare the filters in the standard configuration of the ETKF, experiments with deterministic ensemble transformations are conducted. Experiments including a random rotation are then performed to compare the filters in the standard configuration of the SEIK filter.
The algorithms are applied in identical twin experiments using the model by Lorenz (1996), denoted below as the L96 model, which has been further discussed by Lorenz and Emanuel (1998). The L96 model is a simple nonlinear model that has been used in several studies to examine the behavior of different ensemble-based Kalman filters (e.g., Anderson 2001; Whitaker and Hamill 2002; Ott et al. 2004; Sakov and Oke 2008). Here, the same configuration as used by Janjić et al. (2011) is applied. The model state dimension is set to 40. It is small enough to allow for the successful application of the filters without localization for reasonably small ensemble sizes (see e.g., Sakov and Oke 2008). In our experiments, the localization mainly allowed for the use of smaller ensemble sizes compared to the global analysis, while the relative behavior of the filters was the same as without localization. Thus, for simplicity, only results for global filters are discussed below. The model as well as the filter algorithms are part of the release of the Parallel Data Assimilation Framework (PDAF; Nerger et al. 2005b, available online at http://pdaf.awi.de).
For the twin experiments, a trajectory over 60 000 time steps is computed from the initial state of constant value of 8.0, but with x20 = 8.008 (see Lorenz and Emanuel 1998). This trajectory represents the “truth” for the data assimilation experiments. Observations of the full state are assimilated, which are generated by adding uncorrelated random normal noise of unit variance to the true trajectory. The observations are assimilated at each time step with an offset of 1000 time steps to omit the spinup period of the model.
The initial ensemble for all experiments is generated by second-order exact sampling from the variability of the true trajectory (see Pham 2001). Identical initial ensembles are used for all filter variants.
All experiments are performed over 50 000 time steps. The ensemble size, as well as the forgetting factor, are varied in the experiments. For the ETKF, the covariance inflation is also expressed in terms of the forgetting factor [i.e., γ = ρ−1 is used in Eq. (5)]. Following the motivation of the SEIK filter as a low-rank filter, the ensembles used here are of a size that is at most equal to the state dimension.
Ten sets of experiments with different random numbers for the initial ensemble generation are performed for each combination of ensemble size and forgetting factor to assess the dependence of the results on the initial ensemble. The performance of the filters is assessed using the root-mean-square (RMS) error averaged over the 50 000 time steps of each experiment. The RMS errors are then averaged over each set of 10 experiments with different random numbers for the ensemble generation. We refer to this mean error as MRMSE. Note that the full length of the true trajectory is only used to generate the initial ensemble. For the computation of the RMS errors, only the time steps 1001 to 51 000 of the true trajectory are used.
b. Results with deterministic ensemble transformations
First, the performance of the filters is studied when deterministic ensemble transformations are used. This is the common configuration for the ETKF. In this case, the rotation matrix Λ in Eq. (9) of the ETKF is the identity. In the SEIK-orig, SEIK-sqrt, and ESTKF formulations, the deterministic matrix defined by Eq. (24) is used. For the SEIK-orig filter, a Cholesky decomposition of in Eq. (20) is applied, while the other filters use the symmetric square root.
The left column of Fig. 1 shows the MRMSE for the four filter variants as a function of the forgetting factor and the ensemble size. Filter divergence is defined for an MRMSE larger than one. A white field indicates a parameter set for which the filter diverges in at least one of the 10 experiments.
The ETKF and SEIK-sqrt methods provide almost identical results, with some differences mostly close to the edge to filter divergence. The differences between the results from the ETKF and the ESTKF are even smaller. While mathematically, both variants are identical, the numerical results differ slightly close to the edge to filter divergence. Here, the results of each set of 10 experiments with different random numbers show a larger variability. Thus, the behavior of the filters is less stable in this region and small differences can lead to significant differences. For example in the case with m = 40 and a forgetting factor of 0.99, the ESTKF still converges, while the ETKF diverges. However, the divergence occurs only in 3 of the 10 experiments, which is counted as divergence in the computation of the mean MRMSE. The differences in the MRMSE for the ETKF and ESTKF result from the distinct analysis formulations of both filters. These become visible with the finite numerical precision of the computations over the long assimilation experiments of 50 000 analysis steps. When one considers only the first analysis step, the difference between the transformation matrices is of . The differences in the ensemble transformation matrices of ETKF and SEIK-sqrt are of . While these differences are small with a difference up to 2% of the actual values of the transformation matrix, they can lead to a slightly larger deviation of the MRMSE for the SEIK-sqrt from the MRMSE of ESTKF and ETKF.
The behavior of the SEIK-orig is distinct from the other filters. The filter diverges in most cases with a forgetting factor of 0.97 and above. In contrast, the other filters diverge only for a forgetting factor of at least 0.99. In addition, the minimum MRMSE obtained with SEIK-orig using the deterministic is 0.192 in contrast to the MRMSE of about 0.180 obtained with ETKF and SEIK-sqrt. This difference is statistically significant.
c. Results with random ensemble transformations
The original SEIK filter was always described using a random transformation matrix Ω that preserves the ensemble mean and covariance matrix. Here, the performance of the four filter methods is examined using random rotations. Thus, Λ in Eq. (9) is now used as a mean-preserving random matrix. In SEIK-orig and SEIK-sqrt, a random matrix Ω is used (see the appendix for its construction). In the ESTKF a random matrix Ω is only used for the computation of the weight matrix in Eq. (22). Because Λ and Ω have distinct sizes and are generated by different schemes, the random rotations applied in the ETKF will be distinct from those used in the SEIK filters and the ESTKF.
The MRMSE for the four filter variants with random transformations is shown in the right column of Fig. 1. The randomization results in almost identical MRMSE for all four methods. This indicates that the ensembles of the four methods are statistically of equal quality. Significant differences between the four filters only occur close to the edge to filter divergence, where the filters’ behavior is less stable. The fact that the results of SEIK-orig are comparable to those of the other filters shows that the traditional use of the Cholesky decompostion of in Eq. (20) in SEIK-orig does not deteriorate the state estimate.
The smallest obtained MRMSE is 0.1754. Thus, the MRMSE is slightly smaller with random than with deterministic transformations. This behavior is consistent with the findings by Sakov and Oke (2008). The difference to the MRMSE obtained with deterministic transformations is statistically significant.
d. Ensemble quality
The inferior behavior of SEIK-orig in case of deterministic ensemble transformations can be related to a suboptimal representation of the ensemble. The analysis equations of the filter algorithms based on the Kalman filter assume that the errors are Gaussian distributed. Lawson and Hansen (2004) discussed the effects of nonlinearity on the example of the classic EnKF with perturbed observations and the deterministic ensemble square root filter (Whitaker and Hamill 2002). They found that the ensemble distributions remain closer to Gaussian in the case of the stochastic EnKF.
The ensemble quality can be assessed on the basis of the skewness and kurtosis of the ensembles. These statistical moments will be nonzero if the ensembles are non-Gaussian. Table 2 shows the median and the semi-interquartile range (SIQR) of the skewness and kurtosis for experiments with m = 40 and a forgetting factor of ρ = 0.97. The median of the skewness is about equal for all four filters. However, the SIQR is larger for SEIK-orig than for the other filters. Thus, it is more likely that the ensemble is skewed when applying SEIK-orig. Furthermore, the median and SIQR of the kurtosis are much larger for SEIK-orig than for the filters using the symmetric square root. Thus, the ensemble distributions of SEIK-sqrt, ESTKF, and ETKF are closer to Gaussian distributions than the distribution of SEIK-orig. The stronger deviation from Gaussianity of the ensemble for SEIK-orig is frequently caused by outliers.
When random ensemble rotations are applied, the statistics of skewness and kurtosis are almost identical for all four methods. The median of the skewness is about zero with an SIQR of 0.24. The kurtosis has a median of −0.26 with an SIQR of 0.37. Thus, the values of SIQR and median are closer to zero than in the case of deterministic transformations. This behavior can be attributed to the removal of ensemble outliers by the random rotation (see Sakov and Oke 2008; Anderson 2010).
This study examined the singular “evolutive” interpolated Kalman (SEIK) filter. It was shown that the SEIK filter belongs to the class of ensemble square root Kalman filters. In addition, a variant of the SEIK filter was developed that results in ensemble transformations that are identical to those of the ETKF, but has at a slightly lower computational cost. The variant is referred to as error subspace transform Kalman filter (ESTKF) because it explicitly projects the ensemble onto the error subspace and computes the ensemble transformation in this space.
Numerical twin experiments with the Lorenz-96 model and deterministic ensemble transformations showed very similar results for the SEIK filter with symmetric square root and the ETKF. The differences in the results of the ESTKF and the ETKF are significantly smaller except in the parameter region where both filters exhibit unstable behavior. The variations in the results are related to the ensemble transformations performed in the filters. The differences in the ensemble transformations of SEIK and ETKF are very small. The transformations of the ESTKF and ETKF are analytically identical and at the initial time of the experiments also identical up to numerical precision. However, in the full twin experiments the tiny differences grow because of the finite precision of the computations in combination with the nonlinearity of the model.
Using a Cholesky decomposition in the original SEIK filter with deterministic ensemble transformation resulted in higher errors than the application of the symmetric square root. This effect was caused by an inferior ensemble quality. Accordingly, the experiments indicate that for deterministic ensemble transformations, the symmetric square root should be used in the SEIK filter.
The assimilations with random ensemble transformations provided results that were superior to those using deterministic transformations. This effect was due to the fact that with randomization the ensemble statistics were closer to Gaussian distributions, which are assumed in the analysis step of the Kalman filter. In the case of random transformations, the original SEIK filter with Cholesky decomposition provided state estimates of the same quality as the other filter methods. The numerical results are particular for the specific implementation of the filter algorithms as well as the Lorenz-96 model. However, following the analytical considerations, other implementations of the SEIK filter, the ESTKF, and the ETKF should provide similar results.
The findings of this study unify the developments of the SEIK filter with the class of ensemble square root Kalman filters. Furthermore, the newly introduced ESTKF variant of the SEIK filter provides consistent projections between the ensemble space and the error subspace. Together with the ETKF, the ESTKF has the advantage to provide minimum transformations of the ensemble members. If the minimum transformation is not required, the original SEIK filter is also well suited for practical data assimilation applications.
We are grateful to the editor, Dr. Herschel Mitchell, as well as three anonymous reviewers whose comments helped to improve the text. Also, we like to thank Dr. Marc Taylor for carefully proofreading the manuscript.
Generation of Matrix Ω
The generation of the matrix Ω based on random numbers has been discussed by Hoteit (2001) and Pham (2001) as “second-order exact sampling.” With respect to generating a particular deterministic form of Ω, we review its proposed generation. Note that the algorithm to generate Ω results in spherical sigma points discussed by Wang et al. (2004).
Matrix Ω is required to have orthonormal columns. In addition, the columns need to be orthogonal to the vector whose elements are all one. A Householder matrix associated with the vector ai = (ai,1, … , ai,i)T of size i can be used to generate Ω. It is given by
Here, is identical to ai except for the last element, which is .
Using h(ai), the following recursion (see Hoteit 2001) generates a random matrix Ω:
Set Ω1 = a1, where a1 is 1 or −1 with equal probability.
- Recursion: for i = 2, … , m − 1 initialize a random vector ai of unit norm. Then use the first i − 1 columns of the Householder matrix h(ai) in Eq. (A1), denoted by h−, to compute the i × i matrix:
A simple deterministic variant of Ω can be obtained by taking
with am = m−1/2(1, … , 1)T. This is equivalent to choosing in Eq. (A3).