wiki:OMI_nondiagonal_observation_error_covariance_matrices

Version 9 (modified by lnerger, 2 months ago) (diff)

--

Using non-diagonal R matrices with OMI

Contents of this page

  1. Routines to perform the analysis step
  2. Call-back routines handling observations
    1. prodRinvA_pdafomi
    2. likelihood_pdafomi
    3. add_obs_err_pdafomi
    4. init_obscovar_pdafomi
    5. prodRinvA_l_pdafomi
    6. likelihood_l_pdafomi
    7. prodRinvA_hyb_l_pdafomi
    8. likelihood_hyb_l_pdafomi
  3. Hints

This feature was introduced with PDAF V2.3.

The default mode of PDAF-OMI is to use a diagonal observation error covariance matrix R and specifying the observation error variances, i.e. the diagonalof R only. This is in line with the common choice in data assimilation to assume that observation errors are uncorrelated.

However, there are also observation types with significant observation error correlations, which should be represented by a non-diagonal observation error covariance matrix. With PDAF V2.3 support for such nondiagonal R matrices was added to OMI.

PDAF-OMI's support for nondiagonal R-matrices consists in given the user access to the routines that perform operations involve R. This differs with the filter type, e.g. in LESTKF and LETKF a produce of some matrix with the inverse of R ha sto be computed, while in the traditional, perturbed observations, EnKF the matrix R has to be added to some other matrix. For the particle filter and the NETF the computation of the likelihood involves R.

Routines to perform the analysis step

For using a nondiagonal R there is a modified variant of the PDAF-OMI routines PDAFomi_assimilate_* or PDAFomi_put_state_*. Compared to the case of a diagonal R, the routines are slightly less generic. Below we provide an overview of the routines including the links to the interface description of each. The last column gives the name of the additional routine(s) provided in the interface to implement the non-diagonal R-matrix.

Global Filter diagonal R non-diagonal R additional routine(s)
EnKF PDAFomi_assimilate_global
PDAFomi_put_state_global 		PDAFomi_assimilate_enkf_nondiagR
PDAFomi_put_state_enkf_nondiagR 		add_obs_err_pdafomi
init_obscovar_pdafomi
ESTKF PDAFomi_assimilate_global
PDAFomi_put_state_global 		PDAFomi_assimilate_global_nondiagR
PDAFomi_put_state_global_nondiagR 		prodRinvA_pdafomi
ETKF PDAFomi_assimilate_global
PDAFomi_put_state_global 		PDAFomi_assimilate_global_nondiagR
PDAFomi_put_state_global_nondiagR 		prodRinvA_pdafomi
PF PDAFomi_assimilate_global
PDAFomi_put_state_global 		PDAFomi_assimilate_nonlin_nondiagR
PDAFomi_put_state_nonlin_nondiagR 		likelihood_pdafomi
NETF PDAFomi_assimilate_global
PDAFomi_put_state_global 		PDAFomi_assimilate_nonlin_nondiagR
PDAFomi_put_state_nonlin_nondiagR 		likelihood_pdafomi
SEIK PDAFomi_assimilate_global
PDAFomi_put_state_global 		PDAFomi_assimilate_global_nondiagR
PDAFomi_put_state_global_nondiagR 		prodRinvA_pdafomi


Local Filter diagonal R non-diagonal R additional routine(s)
LEnKF PDAFomi_assimilate_lenkf
PDAFomi_put_state_lenkf 		PDAFomi_assimilate_lenkf_nondiagR
PDAFomi_put_state_lenkf_nondiagR 		add_obs_err_pdafomi
init_obscovar_pdafomi
LESTKF PDAFomi_assimilate_local
PDAFomi_put_state_local 		PDAFomi_assimilate_local_nondiagR
PDAFomi_put_state_local_nondiagR 		prodRinvA_l_pdafomi
LETKF PDAFomi_assimilate_local
PDAFomi_put_state_local 		PDAFomi_assimilate_local_nondiagR
PDAFomi_put_state_local_nondiagR 		prodRinvA_l_pdafomi
LSEIK PDAFomi_assimilate_local
PDAFomi_put_state_local 		PDAFomi_assimilate_local_nondiagR
PDAFomi_put_state_local_nondiagR 		prodRinvA_l_pdafomi
LNETF PDAFomi_assimilate_local
PDAFomi_put_state_local 		PDAFomi_assimilate_lnetf_nondiagR
PDAFomi_put_state_lnetf_nondiagR 		likelihood_l_pdafomi
LKNETF PDAFomi_assimilate_local
PDAFomi_put_state_local 		PDAFomi_assimilate_lknetf_nondiagR
PDAFomi_put_state_lknetf_nondiagR 		likelihood_l_pdafomi
likelihood_hyb_l_pdafomi
prodRinvA_l_pdafomi
prodRinvA_hyb_l_pdafomi


3D-Var diagonal R non-diagonal R additional routine(s)
3DVar PDAFomi_assimilate_3dvar
PDAFomi_assimilate_3dvar 		PDAFomi_put_state_3dvar_nondiagR
PDAFomi_put_state_3dvar_nondiagR 		prodRinvA_pdafomi
En3DVar ESTKF PDAFomi_assimilate_en3dvar_estkf
PDAFomi_assimilate_en3dvar_estkf 		PDAFomi_put_state_en3dvar_estkf_nondiagR
PDAFomi_put_state_en3dvar_estkf_nondiagR 		prodRinvA_pdafomi
En3DVar LESTKF PDAFomi_assimilate_en3dvar_lestkf
PDAFomi_assimilate_en3dvar_lestkf 		PDAFomi_put_state_en3dvar_lestkf_nondiagR
PDAFomi_put_state_en3dvar_lestkf_nondiagR 		prodRinvA_l_pdafomi
hyb3DVar ESTKF PDAFomi_assimilate_hyb3dvar_estkf
PDAFomi_assimilate_hyb3dvar_estkf 		PDAFomi_put_state_hyb3dvar_estkf_nondiagR
PDAFomi_put_state_hyb3dvar_estkf_nondiagR 		prodRinvA_pdafomi
hyb3DVar LESTKF PDAFomi_assimilate_hyb3dvar_lestkf
PDAFomi_assimilate_hyb3dvar_lestkf 		PDAFomi_put_state_hyb3dvar_lestkf_nondiagR
PDAFomi_put_state_hyb3dvar_lestkf_nondiagR 		prodRinvA_l_pdafomi

Call-back routines handling observations

The tables above show the additional routine(s) that need to be implemented the nondiagonal R matrix for the different filters and 3D-Var variants. Here we describe the interface of these routines. Generally the required routine(s) for the chosen DA method should be added to callback_obs_pdafomi.F90 and to each observation module, where thisobs or thisobs_l can be included.

Note that a mixed implementation is not possible. If one observation type has a non-diagonal R matrix, the additional routine(s) have to be implemented for all observations, even if an observation type uses a diagonal R. However, for observation wit a diagonal R, the corresponding PDAF-OMI routine to be called in the observation module can be called. Below we provide the information where this routine can be found.

prodRinvA_pdafomi

The interface for this routine is: 

SUBROUTINE prodRinvA_pdafomi(step, dim_obs_p, rank, obs_p, A_p, C_p)

  INTEGER, INTENT(in) :: step                ! Current time step
  INTEGER, INTENT(in) :: dim_obs_p           ! PE-local dimension of obs. vector
  INTEGER, INTENT(in) :: rank                ! Rank of initial covariance matrix
  REAL, INTENT(in)    :: obs_p(dim_obs_p)    ! PE-local vector of observations
  REAL, INTENT(in)    :: A_p(dim_obs_p,rank) ! Input matrix from analysis routine
  REAL, INTENT(out)   :: C_p(dim_obs_p,rank) ! Output matrix

The routine computes the product of the inverse of the observation error covariance matrix with matrix A_p. For the ESTKF this matrix holds the observed part of the ensemble perturbations. The product has to be given as C_p.

For a model with domain decomposition, A_p contains the part of the matrix that resides on the model sub-domain of the calling process. The product has to be computed for this sub-domain, too.

Hints:

  • The matrix A_p relates to the observations of all observation types. Thus, one needs to take the offset of an observation in the observation vector of all observation types into account. This offset is given by thisobs%off_obs_f.
  • The routine does not require that the product is implemented as a real matrix-matrix product. Rather, the product can be implemented in its most efficient form. For example, if the observation error covariance matrix is diagonal, only the multiplication of the diagonal with matrix A_p has to be implemented.
  • The observation vector obs_p is provided through the interface for cases where the observation error variance is relative to the actual value of the observations.
  • The interface has a difference for ESTKF and ETKF: For ETKF the third argument is the ensemble size (dim_ens), while for the ESTKF it is the rank (rank) of the covariance matrix (usually ensemble size minus one). In addition, the second dimension of A_p and C_p has size dim_ens for ETKF, while it is rank for the ESTKF. (Practically, one can usually ignore this difference as the fourth argument of the interface can be named arbitrarily in the routine.)
  • For diagonal R PDAF-OMI uses the routine PDAFomi_prodRinvA in /src/PDAFomi_obs_f.F90 as the routine that is called within the observation module. This routine can serve as a template for an implementation by the user.

likelihood_pdafomi

The interface for this routine is: 

SUBROUTINE likelihood_pdafomi(step, dim_obs_p, obs_p, residual, likely)

  INTEGER, INTENT(in) :: step                ! Current time step
  INTEGER, INTENT(in) :: dim_obs_p           ! PE-local dimension of obs. vector
  REAL, INTENT(in)    :: obs_p(dim_obs_p)    ! PE-local vector of observations
  REAL, INTENT(in)    :: residual(dim_obs_p) ! Input vector holding the residual y-Hx
  REAL, INTENT(out)   :: likely              ! Output value of the likelihood

The routine computes the likelihood of the observations for each ensemble member. The likelihood is computed from the observation-state residual according to the assumed observation error distribution. Commonly, the observation errors are assumed to be Gaussian distributed. In this case, the likelihood is exp(-0.5*(y-Hx)T*R-1*(y-Hx)).

For a model with domain decomposition, resid contains the part of the matrix that resides on the model sub-domain of the calling process. The likelihood has to be computed for the global state vector. Thus some parallel communication might be required to complete the computation.

Hints:

  • The matrix residual relates to the observations of all observation types. Thus, one needs to take the offset of an observation in the observation vector of all observation types into account. This offset is given by thisobs%off_obs_f.
  • The routine is very similar to the routine U_prodRinvA. The main addition is the computation of the likelihood after computing R-1*(y-Hx), which corresponds to R-1*A_p in U_prodRinvA.
  • The information about the inverse observation error covariance matrix has to be provided by the user. Possibilities are to read this information from a file, or to use a Fortran module that holds this information, which one could already prepare in init_pdaf.
  • The routine does not require that the product is implemented as a real matrix-vector product. Rather, the product can be implemented in its most efficient form. For example, if the observation error covariance matrix is diagonal, only the multiplication of the inverse diagonal with the vector resid has to be implemented.
  • The observation vector obs_p is provided through the interface for cases where the observation error variance is relative to the actual value of the observations.
  • For diagonal R PDAF-OMI uses the routine PDAFomi_likelihood in /src/PDAFomi_obs_f.F90 as the routine that is called within the observation module. This routine can serve as a template for an implementation by the user.

add_obs_err_pdafomi

The interface for this routine is: 

SUBROUTINE add_obs_err_pdafomi(step, dim_obs, C)

  INTEGER, INTENT(in) :: step                ! Current time step
  INTEGER, INTENT(in) :: dim_obs             ! Dimension of obs. vector
  REAL, INTENT(inout) :: C(dim_obs, dim_obs) ! Matrix to that the observation 
                                             !    error covariance matrix is added

During the analysis step of the EnKF, the projection of the ensemble covariance onto the observation space is computed. This matrix is provided to the routine as C_p. The routine has to add the observation error covariance matrix to C_p.

The operation is for the global observation space. Thus, it is independent of whether the filter is executed with or without parallelization.

Hints:

  • Matrix C relates to the observations of all types. Thus, for a single obsevation type one has to take the offset of this observtion type in the observation vector into account. For this, PDAF-OMI initializes an array obsdims, which can also be used by a user-implemented routine. See the routine PDAFomi_add_obs_error in /src/PDAFomi_obs_f.F90 for how this can be implemented.
  • The routine does not require that the observation error covariance matrix is added as a full matrix. If the matrix is diagonal, only the diagonal elements have to be added.
  • For diagonal R, PDAF-OMI uses the routine PDAFomi_add_obs_error in /src/PDAFomi_obs_f.F90 as the routine that is called within the observation module. This routine can serve as a template for an implementation by the user.

init_obscovar_pdafomi

The interface for this routine is: 

SUBROUTINE init_obscovar_pdafomi(step, dim_obs, dim_obs_p, covar, m_state_p, &
     isdiag)

  INTEGER, INTENT(in) :: step                ! Current time step
  INTEGER, INTENT(in) :: dim_obs             ! Dimension of observation vector
  INTEGER, INTENT(in) :: dim_obs_p           ! PE-local dimension of observation vector
  REAL, INTENT(out) :: covar(dim_obs, dim_obs) ! Observation error covariance matrix 
  REAL, INTENT(in)  :: m_state_p(dim_obs_p)  ! PE-local observation vector 
  LOGICAL, INTENT(out) :: isdiag             ! Whether the observation error covar. matrix is diagonal

The routine has to initialize the global observation error covariance matrix covar. In addition, the flag isdiag has to be initialized to provide the information, whether the observation error covariance matrix is diagonal.

The operation is for the global observation space. Thus, it is independent of whether the filter is executed with or without parallelization.

Hints:

  • Matrix covar relates to the observations of all types. Thus, for a single obsevation type one has to take the offset of this observtion type in the observation vector into account. For this, PDAF-OMI initializes an array obsdims, which can also be used by a user-implemented routine. See the routine PDAFomi_init_obscovar in /src/PDAFomi_obs_f.F90 for how this can be implemented.
  • The local observation vector m_state_p is provided to the routine for the case that the observation errors are relative to the value of the observation.
  • For diagonal R, PDAF-OMI uses the routine PDAFomi_init_obscovar in /src/PDAFomi_obs_f.F90 as the routine that is called within the observation module. This routine can serve as a template for an implementation by the user.

prodRinvA_l_pdafomi

The interface for this routine is: 

SUBROUTINE prodRinvA_l_pdafomi(domain_p, step, dim_obs_l, rank, obs_l, A_l, C_l)

  INTEGER, INTENT(in) :: domain_p             ! Current local analysis domain
  INTEGER, INTENT(in) :: step                 ! Current time step
  INTEGER, INTENT(in) :: dim_obs_l            ! Dimension of local observation vector
  INTEGER, INTENT(in) :: rank                 ! Rank of initial covariance matrix
  REAL, INTENT(in)    :: obs_l(dim_obs_l)     ! Local vector of observations
  REAL, INTENT(inout) :: A_l(dim_obs_l, rank) ! Input matrix from analysis routine
  REAL, INTENT(out)   :: C_l(dim_obs_l, rank) ! Output matrix

The routine is called during the loop over the local analysis domains. In the algorithm, the product of the inverse of the observation error covariance matrix with some matrix has to be computed. For the SEIK filter this matrix holds the observed part of the ensemble perturbations for the local analysis domain of index domain_p. The matrix is provided as A_l. The product has to be given as C_l.

This routine is also the place to perform observation localization. To initialize a vector of weights, the routine PDAF_local_weight can be called. The procedure is used in the example implementation and also demonstrated in the template routine.

Hints:

  • The routine is a local variant of the routine U_prodRinvA. Thus if that routine has been implemented before, it can be adapted here for the local filter.
  • The matrix A_p relates to the observations of all observation types. Thus, one needs to take the offset of an observation in the observation vector of all observation types into account. This offset is given by thisobs_l%off_obs_l.
  • The routine does not require that the product is implemented as a real matrix-matrix product. Rather, the product can be implemented in its most efficient form. For example, if the observation error covariance matrix is diagonal, only the multiplication of the diagonal with matrix A_l has to be implemented.
  • The observation vector obs_l is provided through the interface for cases where the observation error variance is relative to the actual value of the observations.
  • The interface has a difference for LESTKF and LETKF: For LETKF the third argument is the ensemble size (dim_ens), while for LESTKF it is the rank (rank) of the covariance matrix (usually ensemble size minus one). In addition, the second dimension of A_l and C_l has size dim_ens for LETKF, while it is rank for LESTKF. (Practically, one can usually ignore this difference as the fourth argument of the interface can be named arbitrarily in the routine.)
  • For diagonal R PDAF-OMI uses the routine PDAFomi_prodRinvA_l in /src/PDAFomi_obs_k.F90 as the routine that is called within the observation module. This routine can serve as a template for an implementation by the user.

likelihood_l_pdafomi

The interface for this routine is: 

SUBROUTINE likelihood_l_pdafomi(domain_p, step, dim_obs_l, obs_l, resid_l, likely_l)

  INTEGER, INTENT(in) :: domain_p             ! Current local analysis domain
  INTEGER, INTENT(in) :: step                 ! Current time step
  INTEGER, INTENT(in) :: dim_obs_l            ! Dimension of local observation vector
  REAL, INTENT(in)    :: obs_l(dim_obs_l)     ! Local vector of observations
  REAL, INTENT(inout) :: resid_l(dim_obs_l)   ! Input vector holding the local residual y-Hx
  REAL, INTENT(out)   :: likely_l(dim_obs_l)  ! Output value of the likelihood

The routine is called during the loop over the local analysis domains. In the NETF, as in other particle filters, the likelihood of the local observations has to be computed for each ensemble member. The likelihood is computed from the observation-state residual according to the assumed observation error distribution. Commonly, the observation errors are assumed to be Gaussian distributed. In this case, the likelihood is exp(-0.5*(y-Hx)T*R-1*(y-Hx)).

This routine is also the place to perform observation localization. To initialize a vector of weights, the routine PDAF_local_weight can be called. The procedure is used in the example implementation and also demonstrated in the template routine.

Hints:

  • The routine is a local variant of the routine U_likelihood. Thus if that routine has been implemented before, it can be adapted here for the local filter.
    • The matrix resid_l relates to the observations of all observation types. Thus, one needs to take the offset of an observation in the observation vector of all observation types into account. This offset is given by thisobs_l%off_obs_l.
    • The routine is very similar to the routine U_prodRinvA_l. The main addition is the computation of the likelihood after computing R-1*(y-Hx), which corresponds to R-1*A_p in U_prodRinvA_l.
    • The information about the inverse observation error covariance matrix has to be provided by the user. Possibilities are to read this information from a file, or to use a Fortran module that holds this information, which one could already prepare in init_pdaf.
    • The routine does not require that the product is implemented as a real matrix-vector product. Rather, the product can be implemented in its most efficient form. For example, if the observation error covariance matrix is diagonal, only the multiplication of the inverse diagonal with the vector resid has to be implemented.
    • The observation vector obs_l is provided through the interface for cases where the observation error variance is relative to the actual value of the observations.
    • For diagonal R PDAF-OMI uses the routine PDAFomi_likelihood_l in /src/PDAFomi_obs_l.F90 as the routine that is called within the observation module. This routine can serve as a template for an implementation by the user.

prodRinvA_hyb_l_pdafomi

The interface for this routine is: 

SUBROUTINE prodRinvA_hyb_l_pdafomi(domain_p, step, dim_obs_l, dim_ens, obs_l, gamma, A_l, C_l)

  INTEGER, INTENT(in) :: domain_p             ! Current local analysis domain
  INTEGER, INTENT(in) :: step                 ! Current time step
  INTEGER, INTENT(in) :: dim_obs_l            ! Dimension of local observation vector
  INTEGER, INTENT(in) :: dim_ens              ! Ensemble size
  REAL, INTENT(in)    :: obs_l(dim_obs_l)     ! Local vector of observations
  REAL, INTENT(in)    :: gamma                ! Hybrid weight
  REAL, INTENT(inout) :: A_l(dim_obs_l, dim_ens) ! Input matrix from analysis routine
  REAL, INTENT(out)   :: C_l(dim_obs_l, dim_ens) ! Output matrix

The routine is called during the loop over the local analysis domains. In the algorithm, the product of the inverse of the observation error covariance matrix with some matrix has to be computed. This matrix holds the observed part of the ensemble perturbations for the local analysis domain of index domain_p. The matrix is provided as A_l. The product has to be given as C_l.

This routine is also the place to perform observation localization. To initialize a vector of weights, the routine PDAF_local_weight can be called. The procedure is used in the example implementation and also demonstrated in the template routine.

The routine is analogous to prodRinvA_l_pdafomi, but also has to apply the hybrid weight gamma of the hybrid filter LKNETF. This is a simple multiplication with the input value in the loop where C_l` is initialized.

Hints:

  • This routine is a simple extension of prodRinvA_l_pdafomi. One can implement the hybrid variant by copying this routine and adapting it. gamma is computed inside PDAF and provided to the routine.
  • Please see the further implementation hints for prodRinvA_l_pdafomi.
  • For diagonal R PDAF-OMI uses the routine PDAFomi_prodRinvA_hyb_l in /src/PDAFomi_obs_l.F90 as the routine that is called within the observation module. This routine can serve as a template for an implementation by the user.

likelihood_hyb_l_pdafomi

The interface for this routine is: 

SUBROUTINE likelihood_hyb_l_pdafomi(domain_p, step, dim_obs_l, obs_l, resid_l, gamma, likely_l)

  INTEGER, INTENT(in) :: domain_p             ! Current local analysis domain
  INTEGER, INTENT(in) :: step                 ! Current time step
  INTEGER, INTENT(in) :: dim_obs_l            ! Dimension of local observation vector
  REAL, INTENT(in)    :: obs_l(dim_obs_l)     ! Local vector of observations
  REAL, INTENT(inout) :: resid_l(dim_obs_l)   ! Input vector holding the local residual y-Hx
  REAL, INTENT(in)    :: gamma                ! Hybrid weight
  REAL, INTENT(out)   :: likely_l(dim_obs_l)  ! Output value of the likelihood

This routine is a variant for U_likelihood_l_pdafomi. See the description of this routine for its functionality. The routine is called during the loop over the local analysis domains. The likelihood of the local observations has to be computed for each ensemble member. The likelihood is computed from the observation-state residual according to the assumed observation error distribution. Commonly, the observation errors are assumed to be Gaussian distributed. In this case, the likelihood is exp(-0.5*(y-Hx)T*R-1*(y-Hx)).

This routine is also the place to perform observation localization. To initialize a vector of weights, the routine PDAF_local_weight can be called. The procedure is used in the example implementation and also demonstrated in the template routine.

The routine also has to apply the hybrid weight gamma. This is a simple multiplication with 1-gamma in the loop where Rinvresid_l is initialized.

Hints:

  • This routine is a simple extension of U_likelihood_l. One can implement the hybrid variant by copying this routine and adapting it. gamma` is computed inside PDAF and provided to the routine.
  • Please see the further implementation hints for prodRinvA_l_pdafomi.
  • For diagonal R PDAF-OMI uses the routine PDAFomi_likelihood_hyb_l in /src/PDAFomi_obs_l.F90 as the routine that is called within the observation module. This routine can serve as a template for an implementation by the user.

Hints

For now, we do not have an example code that we could provide. However, for the implementation the follonig hints might be useful:

  • The matrix R could be constructed e.g. from decorrelation length scales or by direct implementation.
  • The routine PDAF_correlation_function can be used to obtain values of a correlation function.
  • The matrix R cannot be stored fully for larger numbers of observations. E.g. for 100 000 observations, the full matric R would have 1010 entries, which corresponds to 80 GB memory at double precision. Thus, one would need compressed ways of storing the matrix or one can implement the products with the inverse of R or the addition of R in form of oeprations.
  • For the localized operations in prodRinvA_l_pdafomi and likelihood_l_pdafomi it should be possible to handle the part of R that relates only to the local observations used in a local analysis domain. However, the inverse of this local matrix part is not identical to the local part of the global inverse. The different might be mitigated by applying the localization weight to R.
  • Examples of implementations are also available without PDAF-OMI, e.g. in the tutorial and in /models/lorenz96

.