Decomposition of likelihoods and techniques for multi-scale data assimilation

Figure 1.  Illustration of Theorem 3.3. Each image displays a matrix, and the shading of individual cells suggests the value of the corresponding matrix element. Our main result is as follows. Given an arbitrary observation covariance $ {\bf R} $ (top left) and projection $ {\bf \Pi}_1 $ (here, $ {\bf \Pi}_1 $ projects onto a 2d subspace), we compute the projected covariance $ {\bf R}_1 $ (top middle, shown in data space). Using only knowledge of $ {\bf R} $, $ {\bf \Pi}_1 $, and $ {\bf U}_1 $, we discover an $ {\bf R}_2 $ (top right) that can be expressed in model or data space. The unique property of $ {\bf R}_2 $ is that it allows one to reconstruct the original observation precision $ {\bf R}^{-1} $ as, essentially, the sum of $ {\bf R}_1^\dagger $ and $ {\bf R}_2^\dagger $. The visualization above is one way to understand our main result. Theorem 3.3 states the equivalent result $ (y - {\bf H} u)^T{\bf R}^{\dagger}(y - {\bf H} u) = (y_1 - {\bf H}_1 u)^T {\bf R}_1^{\dagger} (y_1 - {\bf H}_1 u) + (y_2 - {\bf H}_2 u)^T {\bf R}_2^{(g)} (y_2 - {\bf H}_2 u) $

Figure 2.  Consider an ensemble of ten evenly weighted particles $ w^i = 1/10 $ for $ i = 1, ..., 10. $ The image plots the output of resampling methods after application to this ensemble. The number of resampled particles coming from the $ i $th original particle are shown by the number of markers in the $ i $th cell. The multinomial resampler (red crosses) returns two copies of the second, sixth, and eighth particles, and none of the first, seventh, or tenth. The stratified resampler (blue circles) returns one of each particle, which is clearly the optimal result

Figure 3.  Error and condition number ratios for repeated experiments decomposing $ {\bf R} $ into $ {\bf R}_1, \, {\bf R}_2 $. Left: the reconstruction error, which is just the distance between $ (y - {\bf H} u)^T{\bf R}^{-1}(y - {\bf H} u) $ and $ \left[(y_1 - {\bf H}_1 u)^T {\bf R}_1^{\dagger} (y_1 - {\bf H}_1 u) + (y_2 - {\bf H}_2 u)^T {\bf R}_2^{(g)} (y_2 - {\bf H}_2 u)\right] $, is bounded by $ 10^{-13} $. Right: the condition number for $ {\bf R}_2 $ is typically equal to the condition number for $ {\bf R} $, but can be up to $ 100 $ times larger. The process of generating these results is fully described at the beginning of Section 5

Figure 4.  Prior and data for the first two variables of the proof of concept in Section 5.1. The 20 particles comprising the prior ensemble are shown by the center of each circle. The radius of the circles indicates the prior uncertainty in each particle. The red line shows the data: we observe the horizontal coordinate, but do not observe the vertical coordinate. As the vertical coordinate $ u_2 $ is not correlated with any other variable, the two intersections between the red line and the prior ensemble suggest that the posterior would be bimodal

Figure 5.  Posteriors for the first two variables of the proof of concept in Section 5.1. In the top row, we show the posteriors for a small ensemble. For both ensemble sizes $ N = 20 $ and $ N = 1000 $, the Pf-EnKF algorithm captures the bimodal posterior. It has an ESS around $ 0.22 N $, retains about $ 30\% $ of particles after resampling, and estimates the remaining variables accurately (shown in Table 1). In contrast, the Op-PF requires the larger ensemble size to capture the bimodal posterior. The EnKF (at either $ N $) approximately captures the posterior mean and variance, but concentrates the majority of its' mass in between the two modes of the posterior

Figure 6.  Performance of DA algorithms for the standard L96 problem described in Section 5.2. To better distinguish between the three stable schemes, the RMSE for time 50-550 plots the moving average over 50 assimilation cycles. Mean RMSE and spread are: 0.19 and 0.21 for ETKF; 0.26 and 0.25 for the stochastic EnKF; 0.28 and 0.23 for Pf-EnKF. The mean ESS in the Pf-EnKF is $ 19.85 $ (out of twenty). The horizontal dotted line shows the expected observation error. As all three stable schemes have spread approximately equal to RMSE and significantly less than the observation standard deviation of 1, they appear well tuned. The Op-PF RMSE is similar to not performing assimilation at all

Figure 7.  Space-time plot of typical solution to "storm track" Lorenz '96 model

Figure 8.  RMSE plots for the 'storm track' model. The Pf-EnKF and Op-PF are stable and more accurate than ETKF or the stochastic EnKF. Mean RMSE and spread are: 0.37 and 0.37 for PF-EnKF; 0.40 and 0.31 for Op-PF; 1.36 and 0.13 for stochastic EnKF. The mean ESS in the Pf-EnKF is $ 18 $ (out of twenty). RMSE is plotted unaltered for the burn-in (time 0 to 50), and for clarity is plotted as a moving average over 50 points for time 50-550

Figure 9.  Rank histograms for the plots in Figure 8. In optimal filtering scenarios, the rank histogram would be flat. The peaks at 0 and 1 indicate observations that are outside the span of the ensemble. With the model error (forcing $ F = 5 $ in the forecast model, while $ F = 8 $ generated data), all methods consistently under-estimate the true state, but Pf-EnKF does so the least