Adaptive covariance inflation in the ensemble Kalman filter by Gaussian scale mixtures
Patrick N. Raanes Marc Bocquet Alberto Carrassi
physics.data-an nlin.CD stat.AP



		
		This paper studies multiplicative inflation:
		the complementary scaling of the state covariance in the ensemble Kalman filter (EnKF).
		
		
		Firstly, error sources in the EnKF are catalogued and discussed in relation to inflation;
		nonlinearity is given particular attention as a source of sampling error.
		In response,
		the “finite-size” refinement known as the EnKF-N is re-derived
		via a Gaussian scale mixture,
		again demonstrating how it yields adaptive inflation.
		Existing methods for adaptive inflation estimation are reviewed,
		and several insights are gained from a comparative analysis.
		One such adaptive inflation method is selected to complement the EnKF-N to make a hybrid that is
		suitable for contexts where model error is present and imperfectly parameterized.
		Benchmarks are obtained from experiments with the two-scale Lorenz model
		and its slow-scale truncation.
		The proposed hybrid EnKF-N method of adaptive inflation is found to yield systematic accuracy improvements
		in comparison with the existing methods,
		albeit to a moderate degree.














[                                                                              
  \begin@twocolumnfalse                                                                
    Unmixing urban hyperspectral imagery with a Gaussian mixture model
on endmember variability
    Yuan Zhou, Student Member, IEEE, Erin B. Wetherley,
and Paul D. Gader, Fellow, IEEEY. Zhou and P. D. Gader are with the Department of Computer and Information
Science and Engineering, University of Florida, Gainesville, FL, USA.
E-mail: mailto:yuan,pgader@cise.ufl.eduyuan,pgader@cise.ufl.edu.E. B. Wetherley is with the Department of Geography, University of
California Santa Barbara, Santa Barbara, CA, USA. E-mail: mailto:wetherley@umail.ucsb.eduwetherley@umail.ucsb.edu.
    Received: date / Revised version: date
==============================================================================================================================================================================================================================================================================================================================================================================================================================================================================================


                                                                      
     
                                                                       
  \end@twocolumnfalse                                                                  
]                                                                                        












§ INTRODUCTION



Consider the problem of


estimating the state _k ∈^M
given the observation _k ∈^P,

as generated by:

	
    equation]eqn:state_sde_discrete_k   = (_k-1) + _k   ,        _k ∼(0,_k)   , 
    equation]eqn:obs_eqn_k        = _k + _k   ,        _k ∼(0,_k)   ,

	subeqns]eqn:HMM

for sequentially increasing time index k,

where the Gaussian noise processes, _k and _k,
are independent in time and from each other.

More specifically, 
the Bayesian filtering problem consists of computing and representing (_k|_1:k),
namely the probability density function ()
of the current state, _k, given the current and past observations,
_1:k = {_l}_l=1^k.


In data assimilation (DA) for the geosciences,

the state size, M, and possibly the observation size, P, may be large,
and the dynamical operator, , may be nonlinear
(observation operators that are nonlinear are implicitly included
by state augmentation <cit.>).

These difficulties necessitate approximate solution methods such as the ensemble Kalman filter (EnKF),
which is simple and efficient <cit.>.

 
The EnKF computes an ensemble of N realizations, or “members”,
to represent (_k|_1:k) as a (supposed) sample thereof.
It consists of a forecast-analysis “cycle”
for each sequential time window of the DA problem.
The forecast step simulates the dynamical forecast
eqn:state_sde_discrete for each individual member.
This paper is focused on the analysis step.

Since the analysis only concerns a fixed time, k,
this subscript is henceforth dropped,
as is the explicit conditioning on 
_1:k-1.
Thus, the prior at time k is written (),
and the analysis (posterior) at time k becomes (|) ∝(|) (),
per Bayes' rule.




Denote {_n}_n=1^N the forecasted ensemble
representing (),
and define the prior sample mean and covariance:

	
    equation]eqn:bx   = 1/N∑_n=1^N _n  ,   
    equation]eqn:barB   = 1/N-1∑_n=1^N ( _n - )
    			( _n - )  .

	subeqns]eqn:bx_bP_def_0


The EnKF analysis update can be derived by assuming that

 and  exactly equal the true moments of (),
labelled  and ,
and carefully dealing with rank issues <cit.>.
The posterior then arises as in the Kalman filter,
described by the analysis moments ^a and ,
or a (deterministic, “square-root”) ensemble transformation to match these.



Multiplicative inflation is an auxiliary technique to adjust (typically increase)
the ensemble spread and thereby covariance, initially studied by
<cit.>.


Here, the specific variant studied is that of
multiplying the prior state covariance matrix, ,
by the inflation factor, α>0, ahead of the analysis:

    ↦α  .






The need for inflation may arise from intrinsic deficiencies of the EnKF:
errors due to non-Gaussianity or the finite size of the ensemble.

The technique of localization should be applied as the primary remedy,
but inflation is still generally necessary and beneficial <cit.>.
Inflation may also be necessary as a heuristic but pragmatic treatment for
extrinsic deficiencies, i.e. model and observational errors,
meaning any misspecification of <ref>.


Again, however, it is advisable to exploit any prior knowledge of errors
(bias, covariance, subspace, etc.)
with more advanced treatments before employing multiplicative inflation.
Examples include additive noise <cit.>,
relaxation <cit.>,
and square-root transformations <cit.>.


It is difficult to formulate directives for the tuning configurations of the EnKF with any generality.

Concerning α,
it may be that the accuracy of the EnKF is improved either by
well-tuned inflation (α>1) or deflation (α<1).
For example,
as detailed in <ref>, sampling error promotes the use of inflation. 
By contrast, the consequences of non-Gaussianity are less transparent.



     




Nevertheless, it generally seems reasonable to inflate
because non-Gaussianity yields an error (intrinsic to the EnKF)
adding to other errors.

Similarly, inflating is typically required in conditions of
extrinsic error such as model error
<cit.>.




Further specificity and quantitative guidelines are difficult to deduce.
Therefore, the inflation parameter typically requires application-specific,
off-line tuning for a fixed value,
sometimes at significant expense.

As an alternative strategy,
adaptive inflation aims to estimate the inflation factor on-line.
This also naturally promotes the use of time-varying values.

The EnKF-N
<cit.>
is a refinement of the analysis step of the EnKF
that explicitly accounts for sampling error in  and ,
meaning their discrepancy from the true moments  and ,
which are seen as uncertain, hierarchical “hyperparameters”.
The derivation proceeds from the rejection of the assumption that  and  are exact
<cit.>.
Moreover, when using a non-informative hyperprior for  and ,
the EnKF-N has been shown to yield a “dual” form which
can be straightforwardly identified as a scheme for adaptive inflation
<cit.>.
Its implementation only requires minor add-ons to the (square-root) EnKF,
with negligible computational cost.
In the idealistic context
where model error is absent or accurately parameterized by the noise process,
as detailed by <ref>,
the EnKF-N nullifies the need for inflation tuning,
making it opportune for synthetic experiments.



However, (i) wider adoption of the EnKF-N has been limited by some technically challenging aspects of its derivation.
Moreover, (ii) the idealism of the above context means that the EnKF-N would still be reliant on
ad-hoc inflation tuning in real-world, operational use.





This paper addresses both of the above issues of the EnKF-N.
Firstly, by re-deriving it with a focus on inflation,
<ref> further elucidates its workings.
Then, <ref> reviews and analyses the literature on adaptive inflation estimation.
In contrast to the EnKF-N, these adaptive inflation methods have hyperpriors that are
time-dependent (as opposed to being “reset” at each analysis time) making them suitable
for realistic contexts where model error is present and imperfectly parameterized.
Then, <ref> uses one such method to complement the EnKF-N
and create a new, hybrid method.
Lastly, <ref> presents benchmark experiments
of the various adaptive inflation methods. 
Expressions and properties of the standard parametric s
in use in this paper,
, , , , ,,
can be found in <ref>.









§ IDEALISTIC CONTEXTS AND SAMPLING ERROR IN THE ENKF



Model-error adaptive inflation is considered from <ref> and onward.
By contrast, this section is focused on the effects of sampling error,
as well as its causes, especially nonlinearity.
<Ref> will show how sampling error is partially remedied by the EnKF-N.











































 §.§ Two univariate experiments


Consider the univariate (scalar) filtering problem where
the likelihood (y|x) = (0|x,2)
and dynamical model (x) = √(2) x
repeat identically for each time index,
and the initial prior is (x) = (x|0,2).
This is a computational (rather than estimation) problem for the posterior;
it is highly artificial, with its numeric values set so as to yield a simple solution.
Indeed, as is perfectly computed by the Kalman filter,
the initial posterior is then (x|y) = (x|0,1),
yielding a forecast prior that is identical to the initial prior.
The cycle thus repeats identically through time.

Now consider the same problem except with nonlinear dynamics,
(x), detailed in <ref>.
This model has been designed to preserve Gaussianity
despite being nonlinear: if (x|y) = (x|0,1)
then (x) has the distribution (0,2).


Hence the nonlinear DA problem has
exactly the same solution as the linear one.

However,
as illustrated in <ref>,
applying a deterministic square-root EnKF (without any inflation or other fixes) to the two problems
yields significantly contrasting results.

The initial ensemble is identical for both cases,
consisting of N=40 members drawn randomly from (x).
But, in the linear case, the resulting sampling errors are quickly attenuated,
and the ensemble statistics converge to the exact ones.


By contrast, in the nonlinear case, the jitteriness (sampling error) is chronic.
This demonstrates that sampling error may arise purely due to nonlinearity,
i.e. without actual stochasticity.

Furthermore, note that the true distributions are perfectly Gaussian,
and therefore the EnKF would compute the exact solution if N were infinite.


Thus, even though nonlinearity typically yields non-Gaussianity, this is not always the case.
Hence, the issue of sampling error, even if caused by nonlinear models,
can be analysed and addressed separately from the issue of non-Gaussianity.

An instructive scenario (not shown) of the nonlinear experiment is that in which

the initial ensemble has a mean of 0 and a variance of 2, exactly.
Despite the “perfect” initialization,
sampling errors will still be generated,
as predicted by <ref>.
However, this error is not immediately as big as if
the ensemble were actually randomly sampled from (0,2),
in which case ∼(2,N-1),
and its expected squared error is [ - 2]^2 = 8/,
per <ref>.
Indeed, repeated experiments indicate that it takes about

5 consecutive applications of 
for the ensemble to saturate at a noise level of 8/.
This gradual build-up also reflects the rule-of-thumb
that stronger nonlinearity breeds larger sampling error.








 §.§ Cataloguing the circumstances for inflation


This subsection is summarized in <ref>,
whose rows correspond to paragraphs, as numbered ().































1.
An important property of the EnKF is that
it is a consistent estimator in the linear-Gaussian case
<cit.>:
at each time k,

the EnKF statistics  and 
converge (in probability, as N →∞)



to the true moments,  and .
Clearly, in this context, inflating or deflating will degrade the ensemble estimates.






2.
Stochastic forms of the EnKF employ pseudo-random “observation perturbations”
for the analysis update step.
Similarly, the forecast step may simulate additive or more advanced stochastic parameterizations
of the forecast noise.
With N<∞, this introduces sampling error.

One cause of the typical need for α>1 is the negative bias
of the posterior ensemble covariance matrix <cit.>:

    [()] < (^a)
    	  ,


	
	
	
	
	

where the expectation is taken over the prior ensemble
(or equivalently the covariance, ), and

    = (^-1 + ^-1)^-1 ,
    	
    ^a   = (^-1 + ^-1)^-1  .


In other words, even though [] =,
the nonlinearity (concavity) of ,
as a function of ,
causes a bias.
A related but distinct bias applies for the Kalman gain matrix,
=  ( + )^-1.
Note, though, that the sampling error originates in the prior;
therefore, the prior covariance is the root cause,
and targeting (inflating) ,
rather than  and ,
is more principled.

There is a misconception



that this bias leads to ensemble “collapse”,
meaning that →0 and →0 as k →∞.
But no matter how acute the single-cycle bias is,
its accumulation will saturate,
because it is counteracted by reductions in .

The term “inbreeding” is sometimes used
to refer to the bias eqn:EnKF_bias_3.
However, inbreeding also encompasses two other issues,
namely the introduction of non-Gaussianity and of dependency between ensemble members.
These are caused by the cross-member interaction
that takes place through the EnKF update <cit.>.
It is not quite clear how these effects will
impact the need for inflation in later cycles.

















Analytical, quantitative results on the bias eqn:EnKF_bias_3
have been obtained for the general, multivariate case by
<cit.>.
However, the degree of the approximation is not entirely clear,
the assumption of the ensemble being truly stochastic is unreliable,
and the related correctional methods were only moderately successful.

An alternative approach is that of 15.3 of <cit.>,
where the bias is empirically estimated by using a companion ensemble of white noise.

However, as discussed below <ref>,
a significant drawback of the inflation methods targeting this bias
is that they do not establish a feedback mechanism through the cycling of DA.
Moreover, as shown by the theory of the EnKF-N in <ref>,
even in a single cycle, the observations, , contain information that can improve
estimates of prior hyperparameters “before” utilising  to update the state vector, ,
thereby reducing sampling error and biases.




3.
Deterministic, square-root update forms of the EnKF
(which may also be formulated for the forecast noise <cit.>)
do not introduce sampling error in the mean and covariance.
Yet, with N<∞, sampling errors will arise due to model nonlinearities.
This was illustrated in the experiments of <ref>,
and predicted by <ref>.
As in 2, sampling error will instigate the need for inflation.


Indeed, the bias eqn:EnKF_bias_3
is slightly visible in the nonlinear experiment
of <ref>,
where the covariances,  and , are on average lower
(long-run averages: 1.95 and 0.98) than the true values.



Filter “divergence” is the situation where
the actual error is far larger than expected from .
It cannot occur in the linear context,
except by extrinsic errors <cit.>.
It may, however, arise in nonlinear, chaotic contexts because,
heuristically,
(i) smaller covariances are prone to deficient (relative) growth by the forecast,
creating an instability that
(ii) might not be adequately controlled by the analyses.
Further, the deficiency in growth typically depends on the starting deficiency of the covariance,
a form of positive feedback that makes the cycle even more “vicious”.
The alarming prospect of divergence,
especially in light of the bias eqn:EnKF_bias_3,
favours “erring on the side of caution”, i.e. using α>1.



4.
With a deterministic, square-root EnKF in the linear context,
sampling error can only come from the initial ensemble and,
as was observed in the experiments of <ref>,
it will be attenuated through the filtering cycles.
Thus, except perhaps from an initial transitory period,
it is not advisable to use inflation.
This is not always true in experiments, however,
because numerical instabilities (or countermeasures such as regularization)
may allow for improved accuracy with some inflation.

The attenuation of sampling errors can be explained as follows.
Apart from the erroneous initial covariance,
the square-root EnKF is here analytically equivalent
to the Kalman filter <cit.>.
Thus, the covariance obeys the Riccati recurrence,
which forgets its initial (erroneous) condition,
also in the case of =0 <cit.>.
Hence, convergence (in time k) holds for any N ≥ M,
with a rate independent of N.

Interestingly, a similar analysis reveals that
the choice of normalization factor for the covariance estimator
(usually 1/N-1, or 1/N)
does not impact the asymptotic EnKF-estimated moments (in the linear context):
they always converge to the true moments as k→∞.
This means that the success of the EnKF
does not so much rely on some statistical, single-cycle optimality
or unbiasedness (in , , or ),
but rather on the above insensitivity to the choice of normalization factor.




5.

Decreasing the ensemble size, N,
increases the sampling error,

the bias eqn:EnKF_bias_3,
and the need for inflation.
Worse, if N ≤ M, then the ensemble is said to be rank-deficient;
this is a separate issue from sampling error,
with the grave consequence that
the truth, , will not lie entirely within the ensemble subspace (cf. <ref>).
By operating marginally,
“localization” <cit.>,
can mitigate the rank deficiency.



Localization also diminishes off-diagonal sampling errors (“spurious correlations”),
thus decreasing the need for inflation.

On the other hand, by eliminating prior correlations,
localization affects an overly uncertain prior,
yielding too strong a reduction of the ensemble spread[
	Formally, quantify the reduction via |_P - | = || / |  + |,
	the determinant of the reduction in the variance.
	Localization decreases the magnitude of the off-diagonals of +,
	provided the eigen-structures of the two terms are not too dissimilar.
	Thus, localization increases the denominator,
	hence reducing _P - and the posterior variance.
].

Another consequence of rank deficiency
is the possibility of the Bayesian uncertainty (i.e. potential error)
outside of the ensemble subspace “mixing in”, and adding to, the ensemble subspace uncertainty.


If = 0 and the context is linear,
this interaction is small and transitory.
It then does not seem beneficial to (inflate in order to)
have the ensemble spread match the total (as opposed to the subspace) uncertainty.

By contrast, if >0 <cit.>,
or in the nonlinear context <cit.>,
the interaction will occur, favouring the use of α>1.
In their section 4,  showed that
(scalar/homogeneous) inflation is well-suited to combat this type of error;
this applies for both multiplicative and additive treatments.





Assuming = 0, the long-run (k →∞) rank of the true state covariance, ,
is the number of non-decaying modes (non-negative Lyapunov exponents) of the dynamics,
i.e. the rank of the “unstable subspace”, 0 ≤ n_0 < M.
This correspondence also holds approximately in the nonlinear context,
and means that the rank deficiency of the ensemble may be much less severe than M-N+1
<cit.>.
If this is the case, a duplicate of <ref> applies,
with M replaced by n_0.




Filter divergence will (almost surely) occur if N ≤ n_0,
if localization is not used.
In contrast with 3,
inflation is then futile,
because the divergence is caused by rank deficiency,
regardless of the degree of nonlinearity of the growth.
It could be speculated that nonlinearity will
sequentially “rotate” the ensemble around in the unstable subspace,
and hence effectively encompass it.
However, twin experiments with the 40-dimensional Lorenz model,
such as the data point N=n_0=14 of Figure 6.6 of <cit.>,
do not give credence to this hypothesis.






§ RE-DERIVING THE DUAL ENKF-N VIA A GAUSSIAN SCALE MIXTURE


This sec:Sketch gives a new derivation of the dual EnKF-N.
Subsection sec:Sketch outlines the main ideas.
The details are filled in by the subsequent subsections.





 §.§ Overview of the derivation







Suppose
the Bayesian forecast prior for the “truth” is Gaussian,
with mean  and covariance ;
formally,
( | _1:k-1) = ( | , ),
where the conditioning on past observations has been made explicit again.
Furthermore, assume that the sample {_n}_n=1^N is an “ensemble”,
meaning that its members are independent and
statistically indistinguishable from the truth <cit.>,
having been drawn from the very same distribution.
In short,

    and _n ∼(, ) .















The assumption eqn:const_4 is convenient,
but may be too idealistic in case of severe
inbreeding, non-Gaussianity, and model error.
Conversely, it may be too agnostic in case the ensemble is not fully random,
as discussed in <ref>.
For convenience, assemble the ensemble into the matrix
= [ _1,   … _n,   …  _N ].


Even in the linear-Gaussian context,
computational constraints induce the use of an ensemble to carry the information on the state,
and thus the approximation

    ( | _1:k-1) ≈(|)   ,

meaning the reduction of the information of _1:k-1 to that represented by the forecast ensemble, .
Thus, while in principle (with infinite computational resources)
the “true moments”,  and , are known,
this is not so when employing the EnKF.
Here, all that is known about  and  comes from .


The appropriate response is to consider all of the possibilities;
indeed, since by the above assumptions
(, ,  | ) =  ( | ,  )  (,  | ),
marginalization yields:

	
	
	
	
	





    (|)
    	   =
    	
    	∫_∫_^M ( | , )
    	 (,  | )   ,

where  is the set of M × M (symmetric) positive-definite matrices[
	 is the Euclidean space ^M(M+1)/2
	corresponding to the M(M+1)/2 upper-triangular elements in ,
	restricted to positive-definite matrices (the conic subset wherein > 0).
].
<Ref> says that the “effective prior”,
(|),
is a (continuous) mixture:

	
	
	

the average of the “candidate priors”,
( | , ),
as weighted by the “mixing distribution”,
(,  | ).
Since the distribution of the state, , depends on
the abstract parameters  and  that are themselves unknown,
these are called hyperparameters and this layered structure is called hierarchical.



The standard EnKF may be recovered from the mixture eqn:hier_5
by assuming that the ensemble size is infinite (N = ∞),
in which case the sample mean and covariance,

		
		
		
	

 and  of <ref>,
are exact, implying a mixing distribution of Dirac delta functions:
δ( - ) δ( - ),
and hence the effective prior: ( | , ).

The EnKF-N does not make this approximation,
but instead acknowledges that N is finite (whence the “finite-size” moniker).
The mixing distribution is obtained with Gaussian sampling theory
and a non-informative hyperprior, (,).
For now, N>M is assumed,
in which case ^-1 exists <cit.>.

The connection to inflation comes from noting, as will be proven later,
that <ref> reduces to:

    (|)
    	   =
    	∫_α > 0 ( | , α)
    	 (α | )
    	α  ,

which is a mixture of candidate Gaussians over a scalar, scale parameter, only.


The mixture eqn:hier_14 is illustrated by the orange objects in <ref>.
The candidate (prior) Gaussians are distinguished solely by the scaling,
α, of the covariance, .
Only a finite selection of the continuous family of candidate priors is plotted,
the selection being representative of the mixing distribution, (α | ).

Interestingly, as detailed later,
this yields an effective prior, (|), which is not Gaussian,
but rather a (Student's) t distribution.


The effective posterior,
(|,) ∝(|) (|),
is given by Bayes' rule, i.e. pointwise multiplication.
But the likelihood,

    (|) = (|,)   ,

per <ref>, is Gaussian.
The posterior is then neither Gaussian nor t,
and does not simplify parametrically.
This poses a computational challenge in high-dimensional problems,
and the question of how the posterior (or an ensemble thereof)
is to be computed in practice.

Progress can be made by noting that
the averaging over the prior moments can be “delayed” until after application of Bayes' rule, i.e.

    (|,) ∝∫( | , )
    		 ( | , α)
    	_(,|α,) (α | ) α  .



Thus, the effective posterior can also be seen as the average of the (Gaussian) candidate posteriors,
(|α,,),
each of which is given by the Kalman filter formulae for a given α,

and computable essentially simultaneously for all α.

The by-product of Bayes' rule is the “evidence”, (|α,).
In this context, it is not a constant,
but instead constitutes the likelihood of the mixing parameter, α.

To reflect this,
the candidate posterior curves in <ref>
have not been normalized to integrate to 1,
but instead (|α,) · c.

The constant c has been inserted and set such that
the particular candidate posterior whose mode coincides with that of the effective posterior also shares its height.
This makes it visible that no candidate posterior is fully coincident with the effective posterior.
Nevertheless, it seems a reasonable approximation.
But this candidate posterior corresponds to a candidate prior,
which merely amounts to choosing a particular prior inflation, α_⋆.
The approximation can thus be written:

    (|,)
    	   ≈(|α_⋆,,)
    	
    	  ,

meaning that the integral over the hyperparameter, α,
for the effective posterior eqn:prior_marg_y,
is replaced by using a particular value, α_⋆,
which is chosen after taking into account .
This approximation is a form of “empirical Bayes”,
known as such because the effective prior
is approximated in a way that depends on the observations, .
This may appear to over-use the observations, ,
but it is merely an artefact of the approximation.

Indeed, decomposing the integrand in <ref>
as (|α,,) (|α,)
makes it apparent that α depends on .





A posterior ensemble corresponding to the approximate posterior eqn:post_EB
may be computed using standard EnKF formulae, except with  replaced by the selected value,
α_⋆.
Provided that the 
choice among the approximating Gaussian posteriors is judicious,
it stands to reason that the resulting ensemble yields
an improved analysis compared to that of the standard EnKF.
After all, the standard EnKF chooses its covariance estimate (=) before taking into account .
By contrast, the EnKF-N lets  inform this choice (= α_⋆).

For the same reason,
even though the EnKF-N does not target any particular unbiasedness,
improvement could be achieved compared to the methods targeting “single-cycle unbiasedness”,
described below <ref>.

However, the main asset of the EnKF-N is that its secondary dependence in 
implicitly establishes a negative feedback loop
via the sequential cycling of DA:
if the covariance estimate was too small at time k,
this will likely be detected and adjusted for at k+1.
Moreover, this feedback is “theoretically tuned”:
parameters that may be tuned exist (cf. <ref>),
but none strictly require it.

As will be shown, the inflation prior is centred on 1,
conferring important advantages to the EnKF-N.
However, this anchoring to 1 also reflects the main drawback of the EnKF-N:
the hyperprior is static,
so that no explicit accumulation of past information takes place for the inflation factor,
which otherwise could have been used to account for model error.
Redressing this is the subject of <ref> and onwards.




 §.§ The mixing distribution

This subsection and the next further describe <ref> for the effective prior, (|).
They are largely sourced from textbooks on Gaussian sampling theory and inference,
under the heading of “predictive posterior”:
the probability of another draw, ,
from the same distribution as the sample, 
<cit.>.
The presentation is didactic, giving meaning to intermediate stages.
A concise version is provided by .



The mixing distribution in <ref> is given by:

    (,  | ) ∝(|,) (,)
    	  ,

where (,) is a hyperprior to be specified.
Here, as in ,
the Jeffreys priors are independently assigned to the hyperparameters:

    (, )
    	=
    	 ()  ()
    	∝
    	1 ·^-(M+1)/2  .

This is a prior designed to be as non-informative (agnostic) as possible. 

It may be derived by positing invariance in location and scale
<cit.>.
 also showed the utility of
using a highly informative hyperprior, suitable in contexts with little nonlinearity.
Examples were also given for encoding information such as
climatology or conditional statistics,
resulting in a form of localization.




By the Gaussian ensemble assumption eqn:const_4,

    (|,)    = ∏_n (_n|,)
    	
       ∝^-N/2  e^-∑_n _n - ^2_ /2  ,

where _^2 = ^-1 = (^-1).




	
	
	

Now, writing _n -  = (  -  ) + ( _n -  ),
it can be shown that

    ∑_n _n - ^2_   =
    	
    	
    	
    	N  - ^2_ + (^-1)
    	  .



Combining <ref>
for the mixing distribution eqn:mix_1,
the resulting factors may be identified as:

    (,  | )
    	   =
    	 ( | , /N)
    	_ ( | , ) ( | ,  )
    	_ ( | )  ,

where  is the inverse-Wishart distribution (cf. <ref> of <ref>).



 §.§ Integrating over the mean

Writing the integrand of <ref> as
(,,|) =  ( | , , )  ( | , )  ( | ),
the integral over  becomes trivial, leaving just the latter two factors:

    (|)
    	   =
    	∫ ( | , )
    	 ( | )
    	  ,

of which ( | ) was obtained in <ref>.

Meanwhile, recalling ( | , ) 
and ( | , ) from <ref> respectively,
it may be shown by completing the square in  that

    ( | , )
    	 ( | , )
    	=
    (| N + /N+1, /(N+1))
    	_ ( | , , ) ( | , )
    	_ ( | , )  ,

where = 1 + 1/N.
The underbraces follow by identification
and provide the first factor in <ref>.

Thus,

    (|)
    	   =
    	∫ ( | , )
    	 ( | ,  )
    	  .

It should be appreciated that 
<ref> would be unchanged
if = had been assumed from the start,
except for the slight adjustment of 
and the reduction from N to  in the “certainty” parameter of .
By contrast, as shown in the following,
the uncertainty in  has significantly more interesting consequences.





 §.§ Reduction to a scale mixture




This section derives the scale mixture <ref>.

While conventional,

the assumption “∝” is ill-suited for inflation targeting sampling error,

as it yields an inflation prior with an overpowering confidence,
to the detriment of the likelihood <cit.>.
This assumption is therefore not made.
But then merely defining the inflation parameter becomes challenging.
Clearly, it must be some scalar summary statistic on the “ratio” of  versus ;
possibilities include using the determinant, trace, or matrix norms.
However, the subsequent assignment “=” would represent an artificial approximation.

By contrast, the following definition and developments make no approximations.

Consider a fixed ,
and define the (squared) inflation:


    =  - ^2_/ - ^2_  .





Now, given the ensemble, , the sample moments
 and  are known (fixed),
while ( | ) =  ( | ,  )
per <ref>.
Thus, by the reciprocity of the Wishart distribution
(<ref> of <ref>),
^-1∼(^-1,).
<Ref> can then be applied to yield
1/∼(1, ).
Thus, again by reciprocity (<ref>),

    (|)
    	   =
    	( | 1, )
    	  ,

meaning that  is inverse-chi-square (cf. <ref>),
with location parameter 1 and certainty N-1.






But the  (|)
could also have been derived
by marginalizing (|) over ∈_,


where  denotes (any parameterization of)
the degrees of freedom in  not fixed by ,
i.e.
_ = {∈^M(M+1)/2-1  ;  (,) ∈,  -^2_ = -^2_}.
Formally,

    ∫__(|)
    	
    	=
    	(|)
    	  ,

with  denoting the Jacobian determinant of (,) ↦.
Inserting the s from <ref>:

    ∫__(|, N-1)
    	   =
    	(|1,N-1)
    	  .

Now, the covariance mixture eqn:h50 can be rearranged as:


    (|)
    	∝
    ∫_exp(-^2_)
    	
    	(| N-1/N,  N )
    	  .


The same change of variables



then yields:

    (|)
    	∝∫_>0exp(-^2_)
    	
    (
    	∫__(| N-1/N,  N )
    	)
    	  .

The inner integral can be substituted by comparing it to 
<ref>, yielding:

    (|) 
       ∝∫exp(-^2_)
    	(| N-1/N, N )
    	
       ∝∫ ( - _ |   0, )
    	( | 1, )
    	  .

In conclusion, the covariance mixture of <ref> reduces to a scale mixture.






An alternative, direct proof,
using tricks from complex analysis instead of <ref>,

was given in a preprint version of this paper.



Note that the scale mixture eqn:h61 has been written
using the notational trick
where  acts as a univariate function.
Also, since  is defined via ,
the integrand of <ref> cannot be read as
“(|,)(|)”.














By contrast,
the mixture eqn:hier_14 is obtained by undoing the trick, and defining α =.















 §.§ Ensemble subspace parameterization


Let  be the vector of ones of length N, and _N the N × N identity matrix.
Then the sample moments,
given in <ref>, may be conveniently expressed as:

    = /N  ,      =   ,

where
	= [ _1 -,     … _n -,     …  _N - ] = _^⊥
is the ensemble “anomalies”, 
with _^⊥ =
the orthogonal projector

	
	

onto ()^⊥, the orthogonal complement space to ().




So far it has been assumed that N>M
so that  is invertible (almost surely) and that the ensemble spans the entire state space.
This is unrealistic for geoscientific DA,
where N rarely exceeds 100, while M may exceed 10^9.
More reasonably, it is henceforth assumed
that the support of the forecast  is confined to the ensemble subspace,
i.e. the affine space {∈^M  [ - ] ∈()}.
This assumption is actually conventional,
as it is implied by the standard EnKF's assumption that = and = along with Gaussianity.
The assumption means that the ensemble has sufficient rank.
Thus, one may expect tolerable accuracy of the filter,
even without localization <cit.>.


It is preferable to work with variables that embody the restriction of the assumption
<cit.>;

	
	
	
	
	
	
	

therefore,
with ∈^N,
the following change of variables is done:


    () =  +   .


In terms of the new variable,
the likelihood eqn:lklhd_x may be succinctly written as:

    (|) =  ( | , )   .

with =  - /N the average innovation,
and =  = the corresponding observation anomalies.


For the effective prior eqn:h61,
note that the ensemble members expressed in the coordinate system of  are merely
the coordinate vectors (_n =  + _n, with _n being the n-th column of _N).
Hence, in this coordinate system,
the sample mean is /N, replaceable by zero since =0,
and the sample covariance matrix is _N.
Substituting these for  and  in <ref>
is a shortcut to obtain the effective prior for ;
with α =,

    (|)
    	∝
       ∫α^-g/2(_1/N-1_N |   0, α)
    	(α | , )
    	α  ,

where the presence of α^-g/2 is explained in the following.

Denote g the dimensionality of the nullspace of .
Due to  it holds that g = max(1, N-M), almost surely.
Thus, typically g=1, and the parameterization in  has one direction of redundancy,
warranting careful attention.
The issue is analogous to expressing 1 random variable as the sum of 2,
or indeed expressing N-g random variables as a linear combination of N.
The principle is that regardless of how the probability space is augmented with the redundant degrees of freedom,
once these are marginalized out, one should be left with the original distribution.
 showed that the adjustment of g in <ref> is then required.




 §.§ The saddlepoint form











Denote f the integrand
of the scale mixture eqn:h71.
Expanding the parametric s yields:


    (α,z)
    	   =
    	α^-(N+g)/2-1  e^ - z / 2 α  ,

where all of the dependency in  is contained in

    z()    = (  + ^2 )
    	  .

Defining F(z()) ∝(|) for the effective prior,
<ref> may be restated as:

    (z) =
    	∫(α,z)
    	
    	
    	
    	
    	
    	
    	
    	α  .

It can be seen that the change of variables u = α/z
factors z out of the integral,
yielding


    (z)
    	   =
    	z^-(N+g)/2(1)
    	  ,

or, reverting to the full  notation,

    (|)
    	   ∝(
    		 + ^2
    	)^-(N+g)/2  ,

which is a t distribution

(cf. <ref>),
also called a Cauchy distribution when g=1.
The t distribution is elliptical, like the Gaussian
<cit.>,
but has thick tails,

making it suited for robust inference
<cit.>.




	
	
	
	
	






Unlike ,
here the t distribution form eqn:effective_prior_w of the effective prior
will not be used directly.
Instead, the effective prior (F) will again be expressed
as a Gaussian (G) with inflation.

To that end, note that,
for general functions F and G with ⊷[F] ⊆⊷[G],
there will always exist a function (z) such that
F(z) = G((z)).



To find a suitable G,
consider applying the mean-value theorem to the integral eqn:sm_l1 for a fixed z,
denoting ζ^-1 the particular point for α .
This will not work because the integrational interval,
[0, ∞), is of infinite length.
In place of the length, therefore,
substitute c_1 ^-1 to form:
G (ζ, z)
	=
	c_1 ζ^-1(ζ^-1,z),
where the constant c_1

ensures that the height (and hence image) of G(  .  ,z) is sufficient.


Inserting  from <ref> yields

    G(ζ,z)
    	   =
    	c_1
    	ζ^(N+g)/2  e^
    			- z ζ / 2
    		  ,

This G works well;
indeed, equating <ref> to eqn:pw_zF
immediately yields the associated function (z) = c_2/z.

In summary, the effective prior
may be expressed by G eqn:Gzz which,
similarly to the integrand, f,
is Gaussian in .
For later optimization purposes,
c_2=N+g is set,
and the logarithm is taken:

	
    (|)
    		   = 
    		c_1
    		
    		
    		exp( L )
    		  ,
    		equation]eqn:Lap_aprx_0
    
    		L(ζ,)
    		   =
    		(  + ^2 ) ζ
    		- (N+g) logζ  ,
    		equation]eqn:g_fawef
    ζ()
    		   =
    		N+g/ + ^2equation]eqn:zeta_w_0  ,

	subeqns]eqn:lap_aprox_all



The value of c_2 was chosen so as to yield the property that
Lζ(ζ(),) = 0 for any ,
as can be directly verified.
Conversely, this means that ζ may be treated as a free variable to be optimized for,
because <ref> is satisfied wherever Lζ = 0.
This tactic becomes useful in the following section.



The above form of the effective prior
may be derived (as in a preprint version of this paper)
as a “saddlepoint approximation”.
Here, however, there is no approximation.
Its exactitude is a remarkable feature
known to arise in a few cases <cit.>.



































































































































 §.§ The posterior and its mode


































Define J = -2 log(|,) plus a constant,
where the posterior is given by Bayes' rule
with the likelihood eqn:lklhd_w
and the effective prior eqn:lap_aprox_all.

The log posterior reads:

	
    J(ζ,)
    		=
    		L(ζ,)
    		-2 log ( | , )
    		+ c
    		equation]eqn:post_4
       =
    		
    		ζ - (N+g) logζ
    		+ ζ^2
    		
    		+	 - _^2
    		equation]eqn:post_7  .

	subeqns]eqn:post_47


Completing the square in  yields:

    J(ζ,)
    	   =  - ^a(ζ)^2__(ζ) + D(ζ)
    	  ,

where the quadratic form is specified
by the usual EnKF subspace analysis formulae:

	
    _(ζ)    = ( ζ_N + ^-1)^-1equation]eqn:EnKF-N-w-KF-Cov  , 
    ^a(ζ)    = _(ζ)  ^-1  ,
    		equation]eqn:EnKF-N-w-KF-Mu

	subeqns]eqn:EnKF_subspace

and D should be recognized as
the “dual”, as in <cit.>:

    D(ζ)
    	   = ζ - (N+g) logζ
    	+ ^2_ + /ζ  .


Now, ζ depends on ,
and so the maximization of (|,) is not as obvious as <ref> suggests.
Fortunately, to find a critical point, it suffices to satisfy

	
    J   = 0equation]eqn:Jw0  , 
    Jζ   = 0 equation]eqn:Jz0  .

	subeqns]eqn:J0

This is because
the criteria eqn:Jw0,eqn:Jz0
imply J = J + Jζζ = 0,
where ζ() is given by <ref>,
which is enforced since Lζ=0,
as follows from <ref>.




Now, the first criterion eqn:Jw0

is trivially satisfied by setting
= ^a(ζ) for a given ζ, as seen from <ref>.
But it can also be seen that J = D along the constraint = ^a(ζ), and so 

    Dζ = Jζ = Jζ + Jζ = Jζ  .

Hence, finding ζ such that Dζ = 0 will satisfy the second criterion eqn:Jz0.


In conclusion, = ^a(ζ_⋆)
is a critical point of the effective posterior
if and only if ζ_⋆ is a local minimizer of D(ζ).
Since both of the terms D(ζ) and - ^a^2__ of
<ref> are here individually minimized,

this critical point
must be a minimum,
as was originally shown using Lagrangian duality theory by <cit.>.
Hence the N-dimensional, non-Gaussian mode-finding problem for (|,)
may be exchanged for the scalar optimization problem in ζ.

The optimization of D(ζ) requires iterating,
but each evaluation of eqn:Dual_def and its derivative
is computationally negligible,
given the singular value decomposition (SVD),

    (σ̅_1,…,σ̅_P)  = []  ,

has been obtained beforehand,
as is typical to compute <ref>.
Multiple minima are a rarity;
in such cases ζ_⋆ will depend on the optimizer and initial guess,
here Newton's method and ζ = N-1.

To obtain an analysis posterior ensemble,
a Gaussian approximation to the effective posterior is chosen.
In addition to its simplicity,
twin experiments 
have provided solid support to using that of ζ_⋆:

    (|,) ≈(| ^a(ζ_⋆),_(ζ_⋆))
    	  .

Notably, this approximation matches the mode and Hessian
of the exact (|,).

The corresponding ensemble is constructed as:

    ^a
    	=
    	[+^a(ζ_⋆)] 
    	+ √()   ,

where = (_(ζ_⋆))
is readily computed using the same SVD eqn:svd_nRY as before,
and may be appended by a mean-preserving orthogonal matrix <cit.>.
Note that this is “just” the symmetric square-root EnKF,
i.e. the ensemble transform Kalman filter (ETKF) of
<cit.>,
except for a prior (squared) inflation of

    α_⋆ = /ζ_⋆  .


As discussed by ,
the choice of candidate posterior is not unassailable,
and certain modifications of the dual function can be argued
on the grounds of modifying this choice.
Indeed, if the influence of the likelihood,
quantified by

    σ̅^2 = (^-1)/P   ,

and computed using the SVD eqn:svd_nRY,
is small:  σ̅^2 → 0,
then it becomes crucial to adjust the choice of inflation factor towards 1.

Moreover, weakly nonlinear contexts create relatively little sampling error.
In such cases, the hyperprior may be too agnostic.
This can be corrected for by increasing (tuning)
the certainty of  to a higher value than ,
yielding a Dirac-delta in the limit.

Finally,
since ζ is not actually constant in ,
the Hessian, (_(ζ_⋆))^-1,
should be corrected by subtracting e ^a^a,
where e = 2ζ_⋆^2/N+g,
.
Since this is but a rank-one update,
a corrected transformation can be
computed without significant expense as
(_N + γ),
where = ^a, γ = e/(τ^1/2 + τ),
and τ = 1 - e,
as may be deduced from equation (B2) of <cit.>.
However, this is typically a very minor correction,
while the Gaussian posterior remains but an approximation.

None of these adjustments were deemed necessary to use in the numerical experiments of this paper.






§ OVERVIEW: ADAPTIVE INFLATION



The filtering problem as formulated with <ref>
uses additive noise processes _k and _k
to represent model and observation errors.
“Primary” filters <cit.> assume the noise/error covariances,  and ,
are perfectly known.
In reality, these system statistics are rarely well known;
this deteriorates the state estimates and can cause filter divergence.
Self-diagnosing and correcting filters
which estimate and modify the given noise statistics are called “adaptive”.
Adaptive Kalman filter literature include <cit.>.



Adaptive particle filter literature include <cit.>.


This section is concerned with adaptive filtering for DA,
focusing on the EnKF and the on-line[
	“On-line” means that the estimation is included in the
	DA cycling loop, so as to be updated each time new data  is received.
	Some off-line approaches not further reviewed include
	<cit.>.
]
estimation of the model error, i.e. .[
The methods may also possess some skill in dealing with errors
due to non-Gaussianity, biases, and even misspecification of .
]
Among the literature using full
<cit.>
and/or diagonal
<cit.> covariance parameterizations,
some success has been noted.
However, the scope of this paper is restricted to the estimation of a single
multiplicative inflation factor, β, for .

This assigns a structure to the covariance and reduces the dimensionality and complexity
of the problem, thus regularizing it.
The tradeoff is a bias,
but this drawback is largely offset by spatialization.

Since this section presents adaptive inflation for model error,
the inflation factor is here labelled β.
This contrasts it to α of the EnKF-N, which targets sampling error.
Formally, the unknown, β, is defined by the modelling assumption
that the ensemble is now drawn with a covariance
that is β times too small,

    _n ∼(, /β)   ,

as compared to the distribution eqn:const_4 of the truth.
Moreover, in order to neglect sampling errors,
the ensemble is assumed infinite (N=∞).
Again, this is commonly tacitly assumed in the adaptive EnKF literature.
The reason, as discussed below <ref>,
is that it yields a prior that is
(the limit of the t distribution which is)
Gaussian.
In summary,[
For familiarity,
the state distributions are written in terms of the original variable ,
rather than the subspace variable  of <ref>.
Note that most of the following inflation estimators
also implicitly address the rank deficiency issue,
as they do not ignore components of the innovation outside of the (observed) ensemble subspace.
This may not be generally beneficial (cf. 5 of <ref>).


]

    (|β, )
    	   =
    	(| , β)
    	  .



The N=∞ assumption is rolled back in the bias study of <ref>.
More pragmatically, <ref> makes a hybrid of the EnKF-N inflation
with a method of adaptive inflation for model error.
The following review and analyses serve to choose a method with which to make this hybrid.


The review is split between methods that may be termed
“marginal”, working with β separately from ,
and “joint”, working with (β,) simultaneously.
The joint approach, including “variational” and “hierarchical” methods,
is theoretically appealing.
However, on closer inspection, including numerical testing,
it was found to be less advantageous.
Therefore, the marginal methods take centre stage,
while the review of the joint methods has been placed
in <ref>.

	


 §.§ Survey: marginal estimation


Recall the average innovation,
=  -.
For brevity, the explicit conditioning on the ensemble, , is henceforth dropped.




  §.§.§ Taking the trace of covariances

Writing = ( - ) + ( - ) + ( - )
it can be seen that

    [ |  ] =  +   .

A departure from non-ensemble works <cit.>
is the adjustment = 1 + 1/N,
resulting from (necessarily) using
the prior ensemble mean, ,
rather than the true prior mean, ,
to define the innovation, .
However, since N=∞ is here assumed, = 1.

Substitute the expectation [ |  ]
by the observed value  in <ref>,
and  by β, per <ref>.
Then,

    ≈β + (β)
    	  ,

which suggests matching (some univariate summary of) 
by adjusting the inflation, β.
For example, <cit.>
working in the framework of the ETKF,
use the estimator:

    β̂_
    	=
    	^2 - ()/()   ,

which follows directly from the trace of <ref>.
Alternatively, the estimator of <cit.>
can be obtained by transforming <ref> by ^-1 before taking the trace,
thus allowing for heterogeneous observations and different units. This yields:

    β̂_ =
    	
    	
    	
    	
    	^2_/P - 1/σ̅^2  ,

where (^-1) = _^2
has been used and σ̅^2 is defined in <ref>.

	
	
	


<cit.> 
proposes a variant using the Schur product ∘^-1.
However, the regular algebra behind β̂_ is preferable,
as it decorrelates the diagonals of 
and hence diminishes the variance of their trace.
More variants are studied in <ref>.
However, <ref> show them
to have more bias than β̂_
because of further exposure to the uncertainty in ,
which is present when the N=∞ assumption is not made.




  §.§.§ Maximum likelihood estimation


As an alternative to the trace-based estimators,
consider some results of likelihood maximization.
Recalling <ref>,
it can be shown that the likelihood for β is:

    (|β)
    	   =
    	( | , )
    	=
    	( | 0, )
    	  ,

where = (β) of <ref>.
By this Gaussianity, 
 
can be seen as the maximum of the likelihood for ,
but only in the case of univariate observations (P=1).
In the multivariate case, the maximization is only defined
restricted to the direction of .







Less artificially,
the maximization can be rendered well-posed
by restricting  to the scaling: (θ) = θ_0 for some _0,
in which case the maximum likelihood (ML) estimate of θ is ^2__0/P,
as noted by <cit.>.



Unfortunately, with the parameterization
(β) = (β) of <ref>,
the ML estimate, β̂_ML, is not analytically available,
and will require iterations
<cit.>.

Also note that <ref> may be derived in the variational framework,
where it is seen as a consistency criterion on the cost function
<cit.>.
These references are also known for deriving further diagnostics,
employing distances labelled “observation-analysis” and “analysis-background”,
which enable the simultaneous estimation of the observation error covariance.
<cit.> make use of this in an ensemble framework,



as do <cit.>,
who makes use of the scheme in a relaxation variety.





  §.§.§ Secondary filters

The accumulation and memorization of past information has gradually become more sophisticated.
<cit.> take the geometric average of β̂_ over time,
while <cit.> use the median of instantaneous maximum likelihood estimates.
Lacking the fuller Bayesian setting, neither yields consistency (convergence to the true value) in time,
because a temporal average of some point estimate β̂_k based on (_k|β)
is generally not the value to which β|_1:K converges as K →∞.
Nevertheless, this mismatch is not likely to be severe,
and the simplicity of the approach is an advantage.


The approach of temporally averaging likelihood estimates has since
been replaced by the more rigorous approach of filtering.
It should be noted that this “secondary” filter is valid
because the innovations are supposedly independent in time <cit.>.
<cit.> and <cit.>
assign a Gaussian prior (β) = (β|β^f, V^f),
where the mean β^f is a persistence or relaxation of the previous analysis mean,
and where the variance V^f is a tuning parameter.
The likelihood is also assumed Gaussian:
(|β) = (β̂|β,V̂),
where β̂ = β̂_ in <cit.>
and β̂ = β̂_ in <cit.>.
In the former, V̂ is a tuning parameter,
while the latter suggests using the variance of β̂_ of <ref>:

    V̂ =
    	[ (^-1)/P + 1/σ̅^2 √(P/2)]^2
    	≈[ β^fσ̅^2 + 1/σ̅^2 √(P/2)]^2
    	  ,

where the approximation comes from using ≈β^f,
which is consistent with operative assumption that N=∞.
Importantly, <ref> has the logical consequence that
the observations are given no weight when
they carry no information, i.e. when
σ̅^2 → 0 or P → 0.

The method is spatialized by
associating each local analysis domain with its own inflation parameter.
Localization tapering promotes smoothness of the inflation field.


Also working in the framework of a square-root EnKF,
<cit.> use diffusive forecasts
for the inflation distribution.
They do not sequentially approximate the posteriors,
instead expressing them explicitly in terms of all of the past innovations,
progressively diffused by a “forgetting exponent”, ϕ<1:

    (β | _1:K)
    	   ∝((β))^ϕ^K∏_k=1^K
    	(_k  | 0, ϕ^k-K_k(β)) 
    	  .

The chosen estimate, which maximizes (β | _1:K), is found iteratively.
The requisite evaluations of the cost function is not prohibitive
because the square-root formulation means that diagonalizing matrix decompositions
have already been computed.



<cit.>, working in the framework of the EAKF,
explicitly considers forecasting the (parameters of the) inflation distribution,
but only tests persistence forecasts, i.e. (β_k | β_k-1) = δ(β_k - β_k-1).
The posteriors (and hence the forecast priors) are approximated by Gaussians:
for a given time,

    (β|y_i)    ≈(β|β̂_MAP,V^a)
    	  ,

which are computed serially for each component y_i of the observation .
The univariate-observation likelihood
allows the maximum  (MAP) estimate β̂_MAP to be found exactly
and without iterations, via a cubic polynomial;
V^a is fitted by requiring that <ref>
be exact for another particular value of β,
as normalized by its value at β̂_MAP.

<cit.>
spatializes the method,
associating each state variable with its own inflation parameter.
Here, correlation coefficients complicate the likelihood,
so that β̂_MAP must be found via a Taylor expansion of the likelihood
and a resulting quadratic polynomial.






  §.§.§ Prior family

<cit.> also notes that a Gaussian prior
does not constrain the inflation estimates to positive values.
Thus, nonsensical negative estimates will occasionally occur for small innovations.
Imposing some lower cap (bound), e.g., β̂_MAP>0 or even β̂_MAP>1,
is done as a quick-fix.
As <cit.> note, such capping may cause a bias.
However, this bias should be avoided by employing the un-capped values in the averaging.


Instead of ad-hoc mechanisms, using a better tailored prior will pre-empt the problem.
<cit.> use an exponential 

for the initial prior, but acknowledge that it is not appropriate for small values.

Indeed, a better choice is the inverse-chi-square distribution, ,
also called inverse-Gamma;
it is the conjugate prior to the variance parameter of a Gaussian sample,
and was shown in <ref> to be intimately linked to inflation.
It has been employed in adaptive filtering by e.g.,
<cit.>,
and in an EnKF contexts by <cit.>
and <cit.>,
who used it to enhance the scheme of <cit.>.
Similarly, the following section formulates an inflation filter
based on β̂_ using  distributions.





 §.§ Renouncing the Gaussian framework



This subsection improves the formulation of the estimator β̂_;
in particular, it abstains from assuming Gaussianity for
the distributions for β.

Recall the definition eqn:bsig_def of σ̅^2

and make the approximation ^-1≈σ̅^2 _P.
This proportionality renders the likelihood symmetric and hence reducible;
indeed, the likelihood eqn:mult_norm_lklhd can now be written:

    (|β)
    	   (^-1/2 | 0, (1 + σ̅^2 β)_P)
    	
       ∝
    	(1 + σ̅^2 β)^-P/2 e^-_^2/2(1 + σ̅^2 β)
       ∝(_^2/P  |  (1 + σ̅^2 β), P)
    	  .

Remarkably, the value of β that maximizes this approximate likelihood is
	
	
	
	
β̂_ of <ref>.

The likelihood eqn:uni_norm_lklhd_3 may be further approximated
by fitting the following shape to it:

    (|β)
    	   ≈(β̂_|β,ν̂)
    	  .

Irrespective of the certainty parameter ν̂,
the mode (in β) is then the same as for <ref>, namely β̂_.
Also fitting the curvature at the mode to that of <ref>
yields ν̂ = P [σ̅^2 β̂_/(1 + σ̅^2 β̂_)]^2.
Remarkably,
this yields a variance (cf. <ref>) equal to V̂ of <ref>,
except with β̂_ in place of β^f.


The benefit of making second approximation eqn:b_lklhd_approx
in addition to the first eqn:uni_norm_lklhd_3
is the resulting conjugacy with an  prior.
Indeed, suppose the forecast prior for the inflation parameter is:

    (β)
    	   =
    	(β | β^f, ν^f)
    	  ,

for some (β^f, ν^f).
The posterior is then:

    (β| )    = (β|β^a, ν^a)
    	  .

where

	
    equation]eqn:nu_postν^a   = ν^f + ν̂  ,
    		
    equation]eqn:beta_postβ^a   = (ν^fβ^f + ν̂β̂_)/ν^a  .

	

This weighted-average update for the parameters is the same as that of the Kalman filter,
except with slightly different meaning to the parameters.
As such,
it constitutes a natural and original derivation of the inflation filter of <cit.>.







 §.§ Potential improvements


Rather than approximating the likelihood as ,
an improved approximation could be obtained by
reverting to the likelihood (|β) of <ref>
and directly fitting (β|β^a, ν^a)
to the posterior (|β) (β|β^f, ν^f)
by matching their modes and local curvatures.
As described for <ref>,
fitting posterior criteria is the approach taken by <cit.>.
With the  prior, though,
an immediate benefit is that of precluding negative and nonsensical
values for β^a.
Moreover, it can be shown that the posterior mode satisfies a cubic equation.
However, the curvature is more complicated; 
alternatives include setting ν^a = ν^f+ν,
where ν = ν̂
or where ν is such such that a likelihood
(β̂_|β,ν) would yield the aforementioned cubic-root mode.

The numerical approach offers another set of options:
locating the mode by optimization and computing the curvature by finite differences.
The latter could also be exchanged for the ratio of two points,
as in <cit.>.
If such avenues are pursued,
it is then not necessary to make the approximation of
reducing the likelihood eqn:mult_norm_lklhd to eqn:uni_norm_lklhd_3;
indeed, it is feasible to find the required statistics with the full likelihood,
provided the pre-computed SVD eqn:svd_nRY.
However, this is not necessarily advisable, because it yields substantial bias
due to the uncertainty in , as discussed in <ref>.

An alternative approximation is to fit (β^a, ν^a)
via the mean and variance of the posterior, as computed by quadrature.
This requires careful implementation to avoid a drift due to truncation errors,
which otherwise accumulate exponentially through the DA cycles.
It is also important to judiciously define the extent of the grid.

Another option is to abandon the parametric approach altogether,
and instead represent the posterior on a grid <cit.>.
In the case of a single, global, inflation parameter, as in this paper,
this approach is affordable for any purposeful precision.
Lastly, a Monte-Carlo representation can also be employed <cit.>.


For all variants, a choice must be made as to which point value of β to use to inflate the ensemble.
Rather than using the parameter β^a directly, one can use the mean or the mode (cf. <ref>).
Typically, though, ν^a is so large that this does not matter much.
Similarly, although β̂_ has a slight bias (<ref>),
it errs on the side of caution.
For this reason or other,
de-biasing did not generally yield gains in testing by twin experiments.
Therefore, and for simplicity, de-biasing is not further employed.




 §.§ Forecasting


The following is a simple and pragmatic
modelling of the forgetting of past information,
achieved by “relaxing” towards
some background (and initial) prior, (β|β_0, 0):

    (β)
    	   =
    	(β | β^f, ν^f)
    	  ,

where 

	
    equation]eqn:nu_fν^f   = e^-Δ/Lν^a  ,
    		
    equation]eqn:s_fβ^f   = e^-Δ/L (β^a  - β_0) + β_0 
    		  ,

	

with Δ as the time between analyses.
The time scale, L ≥ 0, controls the rate of relaxation,
and could be set as a multiple of the time scale of the model.
Alternatively, it can be set by solving the stationarity condition
ν_∞ = e^-Δ/Lν_∞ + ν̂,
derived from <ref>.


The forecast eqn:s_f,eqn:nu_f was engineered
– not derived from dynamics.
While the issue is largely academic,
this lack of formality has been perceived as a difficulty
<cit.>.
The following construes a few possible resolutions.

One possibility is diffusion.
This would yield similar results to eqn:s_f,eqn:nu_f,
though how the  shape may be maintained is not clear.
However, multiplicative, positive noise seems preferable to avoid negative values.
Computationally, if the parametric distribution is not preserved under the forecast,
then gridded, Monte-Carlo, or kernelized inverse-transform approaches may be employed.
Instead of diffusion, <cit.> suggested
exponentiating (β) in the forecast step (cf. <ref>),
also called “annealing” <cit.>.
This maintains the  shape, but is difficult to motivate physically.
Alternatively, the desired effect of “forgetting”
is perhaps most naturally modelled with an autoregressive model with limited correlation length.
This can be rendered Markovian (to fit with filtering theory) using state augmentation
<cit.>.

It may also be argued that the search for physical dynamics is misguided.
After all, the inflation, β, is already a hyperparameter.
So why should it be more natural to forecast β
rather than the hyper-hyperparameters,
β^a,f and ν_k^a,f as in <ref>?






 §.§ Specification of the variants shown in the comparative benchmarks



Several of the improvements of <ref> were tested with twin experiments.
The results (not shown) indicate that most schemes, including the original one of <ref>,
perform surprisingly well (in terms of filter accuracy)
provided they have reasonable settings of their tunable parameters.
The most plausible explanation is that the inflation filters
are consistent estimators:
as the DA cycles build up, they all eventually converge towards a near-optimal value of inflation.

	
	
	
	
	
	



In view of this parity
it seems logical to opt for the simplest scheme,
namely that of <ref>.


Similarly, instead of <ref>,
the numerical experiments use the fixed value ν^f=10^3,
corresponding to a variance of approximately 2(β^f)^2/10^3, according to <ref>.
Moreover, instead of fitting the likelihood's certainty, ν̂,
it was simply set to 1.
Keeping ν^f fixed is suboptimal in the spin-up phase of the twin experiments,
but this part of the experiment is not included in the time-averaged statistics.
Keeping ν^f fixed also foregoes the interesting possibility of
actually having the certainty decrease through an update,
something that will occur with the posterior fitting or non-parametric approaches.
Nevertheless, this simplification was done
(i) to facilitate reproduction;
(ii) because in the twin experiments of this paper,
the models and observational networks are homogeneous,
and therefore using localized and variable ν̂ and ν^f is not crucial;
(iii) it was found that using the fitted ν̂
sometimes yielded worse filter accuracy than keeping it fixed;
(iv) to provide fair comparisons between all of the methods
by equalizing the sophistication of their forecast step.

In addition to the adaptive inflation,
<cit.> also uses a fixed inflation of 1.015.
This was tested and found to make little difference in the experiments herein.

Without directly impacting its distribution,
the value of the inflation actually applied to the ensemble is capped below by 0.9
for all methods.
However, this clipping never occurred after the spin-up time of the experiment.

The above scheme is practically identical to that of <cit.>,
the differences having proven largely irrelevant and “cosmetic”, as discussed above.
It is therefore labelled “ETKF adaptive”.
Another scheme that was tested, labelled “EAKF adaptive”,
is that of <ref>.
Instead of fitting V^a, however, the value of V^f was fixed, and set (tuned) to 0.01.
As described by <cit.>,
the actual inflation applied to the ensemble is damped, using the value 0.9 β̂_MAP.
The above two schemes are the established standard in the literature,
and have featured in many operational studies,
although mainly in their spatialized formulations.






The proliferation of tunable “hyper-hyperparameters” (the hyperparameters of the inflation filter),
is due to the hierarchical nature of adaptive filters,
and may make the adaptive approach appear counterproductive.
But, considering the breadth of error sources targeted,
it should be recognized that
such methods will necessarily be ad-hoc,
and that the existence of tuning parameters is to be expected.
One should not be dismayed, however,
because the performances of the adaptive filters are largely insensitive to the hyper-hyperparameter settings.
Indeed, intuitively, more abstract parameters (further up the hierarchy) should have less impact,
as is illustrated by the results of <cit.>.
Indeed, the given values for the hyper-hyperparameters
(i) were only tuned for a single experimental context,
(ii) seem reasonable,
and (iii) yield satisfactory filter accuracy almost universally across the experiments.
Point (iii) corroborates previous findings <cit.>,
and suggests that these values may be used in vivo.











§ A HYBRID-INFLATION ENKF-N




<Ref> showed that
the EnKF-N may yield a form of adaptive inflation, α_⋆.
However, the EnKF-N is built on <ref>,
with the assumption that the truth is statistically indistinguishable from the ensemble.
Therefore, it only targets sampling error,
and the prior eqn:lam3 always has the location parameter 1.
The EnKF-N is therefore not robust in the context of model error,
analysed in <ref>.

If tuned inflation is used in concert with the EnKF-N,
then the tuned value could be seen as a measure of model error
disentangled from sampling error <cit.>.
Here, however, the aim is to hybridize the EnKF-N with an
adaptive inflation scheme that estimates an inflation factor, β, targeted at model error.
Notably, such a scheme has a prior that is time-dependent
and, generally, not with location parameter 1.
The scheme used is the “adaptive ETKF” specified in <ref>.

Again, the explicit conditioning on the ensemble, , is here dropped for brevity.
Consider

    (,β | )
    	∝( | ) ( | β) (β)
    	  .

In contrast with <ref> and <ref>,
N will here not be assumed infinite.
Re-deriving, therefore, the EnKF-N prior eqn:hier_14,
but with <ref> in place of eqn:const_4,
reveals that ( | β) is a scale mixture over α,
but now with  also scaled by β,
i.e. ( | α, β) = (|,αβ).
Further, sampling error is assumed independent from model error:
(α | β) = (α).
Hence,

    (,β | )
    	∝( | ) ( ∫( | α, β) (α) α) (β)
    	  .

Moving the likelihood inside as in the EnKF-N eqn:prior_marg_y
yields a mixture over (,  | α, β),
which can be re-factorized to obtain:

    (,β | )
    	   ∝( ∫( | , α, β) ( | α, β) (α) α) (β)

As in the EnKF-N eqn:post_EB, the mixture is approximated by empirical Bayes:

    (,β | )
    	   ( | , α_⋆, β)
    	_≈( | , β)( ∫( | α, β) (α) α_( | β)) (β)
    	  ,

meaning that ( | α_⋆, β, ) is given by <ref>,
with the change of variables eqn:CVar_w,

except that the prior covariance is now also scaled by β,
which also impacts the selection of α_⋆ through the dual cost function.

The remaining integral of <ref> is again approximated using a particular value of α:

    ∫( | α, β) (α) α   ≈( | α^f, β)
    	  .

In this instance, however, the optimizing value of α is not conditioned on ,
and is therefore denoted α^f.
In practice, α^f = 1 is used for simplicity.
Thus, ( | α^f, β) becomes the same as in <ref>,
yielding the posterior of <ref>.

A final approximation is made to decouple the joint posterior:
replacing β in the conditional distribution, ( | , α_⋆, β),
by β_⋆, some point estimate from (β | , α^f).
This may again be seen as empirical Bayes,
except that β is not a latent variable.
Other reasons for only using a single inflation value for the ensemble
are discussed in <ref>.
The particular point used is the mean:

    β_⋆ = ν^a/ν^a-2β^a  .

Thus, <ref> becomes:

    (,β | )
    	   ≈( | , α_⋆, β_⋆)
    	( | α^f, β) (β)
    	
       ∝( | , α_⋆, β_⋆)
    	(β | , α^f)
    	  .


The algorithm of the analysis update of the EnKF-N hybrid
can be stated as follows.

	
  * Update the general-purpose inflation, β, according to <ref>.
	
  * Conditional on β_⋆ eqn:beta_star, find the EnKF-N inflation using the dual
		eqn:Dual_def:
		α_⋆ = _α D(/β_⋆α).
	
  * Update the ensemble by the ETKF with a prior inflation of α_⋆β_⋆.
		In other words, implement <ref> with
		/ α_⋆β_⋆ in place of ζ_⋆.

Note that the second step is the only essential difference to the ETKF adaptive inflation scheme
of <cit.>.

The forecast is greatly facilitated by
the decoupling of the posterior eqn:post_x_beta_24,
which means that the ensemble and general-purpose inflation β are independent.
As discussed by <ref>,
it is pragmatic to separate the forecast of the ensemble
from the forecast of β,
which should be carried out as in <ref>.
Instead, however, for the reasons described in <ref>,
here ν^f = 10^4,
and β^f is set to the previous β^a.
The EnKF-N inflation, α, is not forecasted,
as its prior is static.
























































































§ BENCHMARK EXPERIMENTS


The standard methods of <ref>
and the hybrid of <ref> are tested with
twin experiments:
a synthetic truth and observation thereof are simulated
and subsequently estimated by the DA methods.
Contrary to the meaning of “twin”, however,
the setup is deliberately one of model error:
the model provided to the DA system is different from the one
actually generating the truth.


The main system used is the
two-scale/layer Lorenz model <cit.>.
It constitutes a surrogate system for synoptic weather,
exhibiting similar characteristics,
and enables studying the impact of unresolved scales on filter accuracy.
The autonomous dynamics are given by:
ψ_i^±()
	=
	u_i ∓ 1 (u_i ± 1 - u_i ∓ 2) - u_i,
where the indices apply periodically. 

Then,

    x_it   =
    	ψ_i^+()
    	+ F - h c/b∑_j=1^10 z_j + 10(i-1)  ,
    	
    z_jt   =
    	c/bψ_j^-(b)
    	+ 0
    	+ h c/b x_1+(j-1)//10  ,

for i = 1,…,36, j = 1,…,360, and
where // means integer division.
Unless otherwise stated, 
the constants are set as in <cit.>:
time-scale ratio: c=10,
space-scale ratio: b=10,
coupling: h=1,
forcing: F=10.
The resulting dynamics,
illustrated in <ref>,
are chaotic and have a leading Lyapunov exponent of 1.3775 <cit.>.




The truth is composed of both fields:
[ , ].
By contrast, 
the model provided to the DA methods is a truncated version of the full one:

    x_it   =
    	ψ_i^+()
    	+ F - [A + B x_i]
    	  ,        i = 1,…,36.

The term in the brackets parameterize (compensate for)
the missing coupling to the small scale.
Ahead of each experiment, the constants A and B are determined
by linear regression to unresolved tendencies,
as described by <cit.>.
The error parameterization removes the linear bias of the truncated model;
this has been done because inflation is not well suited to deal with systematic errors.
Higher-order parameterizations were tested and found to yield little improvement to the filter accuracies,
while a 0-order parameterization yielded too much model error,
overly dominating the dynamical growth in the filter error.

Another source of model error is that
the full model is integrated with a time step of 0.005, as necessitated by stiffness <cit.>,
while the truncated model uses 0.05 to lower computational costs.
However, this source of error was found to be negligible compared to the truncation itself.

The time between observations is Δ = 0.15.
Direct observations are taken of the full  field with error covariance = _M.
There is no model noise: = 0.

The filters are assessed by their accuracy
as measured by root-mean squared error:

    = √(1/M - ^2)  ,

which is recorded immediately following each analysis.
This instantaneous RMSE
is then averaged in time (3300 cycles following a spin-up of 40 cycles),
and over 32 repetitions of each experiment.
A table of RMSE averages is compiled for a range of experiment settings

and plotted as curves for each method.
The figures also present (thin, dashed lines) the root-mean variances (RMS spread) of each method.

All of the experiment benchmark results can be reproduced using Python-code scripts hosted online at
<https://github.com/nansencenter/DAPPER/tree/paper_AdInf>.





















































































<Ref> shows benchmarks obtained with the two-scale system as a function of the forcing, F.
The increasing RMSE averages of all filters
reflect the fact (not shown) that the system variability and chaoticity both increase with F.
The same applies for decreasing c in <ref>, where F is fixed at 10.
Note that all of the adaptive filters are largely coincident at F=10 and c=10,
with RMSE scores almost as low as fixed, tuned inflation.
This is because the hyper-hyperparameters for each method,
described in <ref>,
were tuned at this point (and this point only).
By contrast, the fixed inflation of the “ETKF tuned” filter
is determined for each experiment setting by selecting for the lowest RMSE
among 40 inflation values between 0.98 and 3, most of which are close to 1.

It is not surprising that tuning an adaptive filter
will make it about as accurate any other.
The objective, however, is to avoid tuning.
In that regard, it is surprising is how well all of the adaptive filters perform overall.
Indeed, except for the fairly extreme contexts of F>15 or c<4,
the difference in RMSE is small
in the sense that the adaptive filters are all superior to “ETKF excessive”:
a fixed-inflation filter with a suboptimal inflation factor that adds 0.1 to the optimal value.



Benchmarks were also obtained with the single-scale system,
where both the truth and the DA systems are given by:

    x_it   =
    	ψ_i^+()
    	+ F
    	  ,        i = 1,…,40,

and there is no  field.
As in <cit.>,
the model error consists in using
a different value of F for the truth than for the DA system.
The setup is otherwise repeated from above.
As shown in <ref>,
the benchmarks plotted as a function of F are V-shaped,
with the lowest scores obtained in the absence of error (F=8).
The adaptive filters score very similar RMSE averages,
which are generally significantly in excess of the RMS spread scores.
The mismatch can be explained by the well known bias-variance decomposition
of RMSE, and the fact that the model error contains significant bias.
The presence of bias is also a likely cause for the closeness of the
RMSE averages of the adaptive methods,
because inflation is not well suited to treat bias.

Tests were also run with the 3-variable Lorenz-63 system
<cit.>,
where the model error consists in adding independent white noise to the truth.
The setup is the same as above,
except that = 2 _3.

<Ref> shows the corresponding benchmarks.
These are obtained with a small ensemble (N=3);
using a larger ensemble, the relative advantage of the hybrid disappears.


The hybrid EnKF-N obtains slightly superior accuracy relative to
the adaptive ETKF and EAKF for nearly all experiments.
This is as expected from theory: separate, dedicated treatment of
sampling and model errors yield improved accuracy.
The practical advantage of the hybrid is illustrated in the
time series of <ref>.
Notably, the inflation of the hybrid EnKF-N has much more volatility
(shorter time scale).
This is made possible by the static prior which “anchors” the inflation to 1.
By contrast, similar volatility in the adaptive ETKF and EAKF would
require much more lenient settings of ν^f (or V^f),
which would yield excessive longer-term volatility (i.e. variance).



On the other hand, the volatility also means that larger spikes in the inflation will occur.
Since inflation is not a physical or especially gentle way of increasing spread,
this could potentially cause trouble.
For example, the pure EnKF-N with the full model and large F values,
blows up due to stiffness, as illustrated by the grey curves in <ref>.
This could have been prevented by increasing (doubling) its certainty parameter,
something that would also slightly improve its accuracy across all of the experiment settings.

Different hyper-hyperparameter values for the adaptive EAKF and ETKF
will penalize their RMSE scores for some settings,
and reward it for other settings.
This sensitivity was observed to be much reduced for the hybrid,
which is in line with the principal objective:
to avoid tuning across a multitude of contexts.

A secondary objective is to obtain improved accuracy compared to fixed, tuned inflation,
as has been previously observed for the pure EnKF-N
in the perfect-model context .
<Ref> shows that this is sometimes achieved,
but by a very small margin.

Further experiments (not shown) were carried out,
using different ensemble sizes and other types of model errors.
The trends were similar to the benchmarks already shown,
but typically with less relative difference between the filters.






§ CONCLUDING REMARKS



This paper has developed an adaptive inflation scheme
as a hybrid of (α) the finite-size ensemble Kalman filter inflation

and (β) the inflation estimation conventionally associated with the ensemble transform Kalman filter (ETKF).
In so doing, it has provided several novel theoretical insights on the EnKF and adaptive inflation estimation.

The first part of the paper is focused on idealistic contexts,
with sampling error being the main concern.
Using two univariate toy experiments,
<ref> illustrated the generation of sampling error
by nonlinearity, as predicted by <ref>.
<Ref> then discussed the circumstances for inflation,
cataloguing them according to linearity, stochasticity, and ensemble size.
The discussion revealed
why sampling errors are attenuated in the linear context,
why the choice of normalization factor (e.g., 1/N-1) is not crucial,
and also touched on topics such as ensemble collapse and filter divergence.

Next, <ref> gave a birds-eye view of the EnKF-N,
showing how (e.g., empirical Bayes) and why (e.g., feedback) it works.
The following sections filled in the details;
in particular,
<ref> showed how the effective prior

reduces to a Gaussian scale mixture,
again demonstrating the relationship between sampling error and inflation.
The mixture parameter, α, is shown to be .
<Ref> derived a saddlepoint form
to retain the inflation-explicit expressions all the way up to the posterior.
Without recourse to Lagrangian duality theory,
<ref> then finalized the re-derivation of the (dual) EnKF-N
by showing how the mode of the posterior may be found by optimizing for the inflation factor, α.

In contrast to the above, <ref> is focused on model error,
neglecting sampling error.
A formal and unifying survey of the existing adaptive inflation literature is presented;
particular attention is given to the schemes
conventionally associated with the EAKF and the ETKF (β̂_).
The ETKF scheme is given a new and natural derivation in <ref>,
again yielding the  distribution.
Several potential improvements, some novel, were discussed in <ref>,
but generally found to be less rewarding than hoped for.
<Ref> gives some new results on biases, including the maximum likelihood estimator.
The survey is supplemented by
<ref> on joint schemes,
including a suggestion for inflation estimation by variational Bayes.
<Ref> commented on the forecasting of hyperparameters such as the inflation parameter, β.


Combining the above, <ref> developed
a hybrid between the EnKF-N and the adaptive ETKF inflation scheme.
The hybrid employs two inflation factors, α and β,
separately targeting sampling and model error, respectively.
The EnKF-N component (α) adds
negligible computational cost and no further tuning parameters
to the ETKF and its adaptive inflation (β),
yet increases the inflation volatility and thus ability.

The experiments of <ref> showed
that the hybrid generally yields similar filter accuracy as fixed inflation,
even in bias-dominated contexts,
but without the costly need for tuning.
It also yields improved filter accuracy in comparison with the standard, pre-existing
adaptive inflation schemes of the ETKF and the EAKF.

Unless the ensemble size was small and the context strongly nonlinear,
however, the gains were found to be relatively modest,
as was the difference in between the existing methods.
This is somewhat surprising in view of
the essential importance of inflation in many configurations of the EnKF.
Part of the explanation may be that, as a hyperparameter,
the accuracy of the inflation estimates is not as important as that of the (primary) state variables and that,
instead, the main importance of the inflation scheme consists in its capacity to avoid divergence occurrences,
which is a matter of a more boolean character.
Another cause is that the inflation estimates converge and become nearly constant within a relatively short span of time,
and that these asymptotic estimates are sufficiently accurate for all of the methods.


While the experimental results clearly demonstrated
the improvements of the hybrid adaptive inflation scheme,
extrapolating these findings to other, larger applications is non-trivial.

Spatialization of the inflation parameter will likely be necessary;
it may be implemented without considerable complexity as in

<cit.>.
Still, the relative modesty of the above experimental results
does not promise great, general gains.
On the other hand it suggests the conclusion,
aided by the rigour and scope of this study,
that further sophistication of single-factor adaptive inflation estimation schemes
is unlikely to yield significant, further improvements.



  
  
    
      
    
      
      
    
      
		
      
    
    
    
      
      
      
    
  

































§ STANDARD DISTRIBUTIONS


<Ref> specifies the distributions in use in this paper.
The following properties are useful.







	
	The (“scaled”) chi-square distributions are equivalent to the Gamma distributions:
	
    ^± 2(β | s, ν)    = Gamma^± 1(β|ν/2,ν s^∓ 1 /2)   ,

	where the switch sign ± has been used to represent both the regular and inverse distributions.
The  parameterization has been preferred
for the parameter interpretations offered by <ref>,
and the notational simplicity of <ref>.


	
	Asymptotic normality.
	If β∼^± 2( s, ν),
	then the distribution of √(ν) (β - s) converges to (0, 2 s^2)
	as ν→∞.
	
	
	
	
	
	
	
	


	Since it describes the sum of squared Gaussians,
	the asymptotic result for  is a consequence of the central limit theorem.
	The result for 
	can be shown by through the pointwise convergence of the  of √(ν) (β - s),
	normalized by its value at 0.


	Note that the same limit would have applied if β∼(s, 2 s^2/ν).
	This shows that s plays the role of a location parameter in ^± 2( s, ν),
	while 2 s^2/ν plays the role of variance,
	and explains why “certainty” is preferred to “degree of freedom” for ν in this paper.


	
	In the univariate case (M=1),
	
    ^± 1(β | s,  ν)    = ^± 2(β | s, ν)   .



	
	Reciprocity. With t = 1/β,
	
    (β)     = (β|s, ν) 
    iff. (t)    = (t|1/s, ν)   .



	
	Reciprocity. With = ^-1,
	
    ()    = (|, ν) 
    iff. ()    = (|^-1, ν)   ,

	as follows by the change of variables and the Jacobian
	^-(M+1) <cit.>.


	
	Let ≠0 be any M-dimensional vector,
	or an (almost never zero) random vector.
	If ∼(,ν) is independent of ,
	
	then
	
    /   ∼(1,ν)
    		  .

	
	Moreover, this statistic is also independent of .
	Proof: theorem 3.2.8 of <cit.>.







§ NONLINEARITY AND SAMPLING ERROR


This discussion complements that of <ref>.



 §.§ Why does nonlinearity generate sampling error?


First, consider what is meant by “sampling error”.
A sample does not per se have a sampling error;
it is by definition random, i.e. subject to variation.
By contrast, estimators, or rather their realized estimates, have sampling error:
the difference between the estimate and its expected value.
By extension, any statistic (any function of the sample) may be said to have sampling error;
for simplicity, however, the discussion below is limited to the non-central sample moments,
i.e. μ̂_m({x_n}_n=1^N) = N^-1∑_n=1^N x_n^m, in the univariate case.
If the sample is drawn from the same distribution as x,
then μ̂_m is an unbiased estimate of
the m-th moment of x, i.e. μ_m = [x^m],
and the sampling error is the difference:

    Error_m    = μ̂_m - μ_m
    	
    	
    	  .




It is a well known property of the Kalman filter that
the covariance does not depend on the mean.
In the forecast step, this is due to the fact that
their evolutions are entirely decoupled.
Indeed, in the case of linear dynamics (d=1),
the m-th forecast moment is given by:


μ_m^f = ^m μ_m,
where  is the (scalar) linear model:

x^f =  x.
For nonlinear forecast dynamics , however,
the moments will be coupled through .



For example, if (locally to the support of (x))
the model  can be represented by a polynomial of degree d,

then the m-th moment of the random variable (x) is a linear combination
of moments of x of order 1 through m d:


    μ_m^f   =
    	∑_i=1^m d
    	C_m,iμ_i
    	  .

Thus, for d>1 the moments get mixed and, in particular,
impacted by moments of higher order.
This is known as the “closure problem” <cit.>.

A similar analysis reveals that the same coupling takes place for the sample moments, μ̂_m.
Therefore the sampling errors are also coupled:

    Error_m^f   =
    	∑_i=1^m d
    	C_m,iError_i
    	  .

But an N-sized ensemble can only match N moments,
e.g., Error_i = 0 for i=1,…,N,
and so there will always be some sampling error present for i>N.
Then, by the mixing of <ref>,
this error will cascade into the lower-order forecast errors.




In summary, nonlinearity causes sampling error
in (e.g.) the mean and covariance
by pulling in the inevitable (with finite N)
sampling error from higher-order moments.

The above analysis is concerned with the generation of sampling error.
Another reason for sampling error in the context of nonlinearity
is that chaos prevents its elimination
by limiting the effect of far-past observations
(as opposed to the linear case illustrated in <ref>,
where the initial sampling error is quickly attenuated).
It is not immediately clear whether this is a separate cause
or, rather, a different perspective on the same phenomenon.



 §.§ A nonlinear model preserving 


The nonlinear model of <ref> was
designed using “inverse transform sampling”.
It is specified by:
(x) = √(2)^-1(  (x^2) ),
where  is the cumulative distribution function (CDF) for (1,1),
and ^-1 is the inverse CDF for (0,1).

	For context, note that :
	(i) is V-shaped, with a singularity at 0,
	(ii) is closely approximated as
	(x) ≈√(2){0.88x + 0.23 log(x^2) -0.4 },
	(iii) applied to a density symmetric about 0, it may be visualized
	as folding it up in the middle before smearing it back out again,
	(iv) may be generalized to higher dimensionality
	by expressing  in polar coordinates,
	since (^2) = (^2|M,M)
	if () = (|0,_M).
	
	


Presumably, any nonlinear model that preserves Gaussianity must include a singularity.
This renders the example using  to generate sampling error without non-Gaussianity
somewhat artificial,
but does not jeopardize the utility of considering the two issues separately.
















































































































































































































§ JOINT STATE-COVARIANCE ESTIMATION



As mentioned immediately above <ref>,
the joint approach approximates (,β|) simultaneously in  and β.
Thus, their analyses impact each other.
This is particularly relevant for state-inflation (or state-covariance) estimation problem,
because of the non-Gaussianity of (,β|).



 §.§ Variational methods


A common approximate solution to this non-Gaussianity
is to use variational methods for parametric fitting.
As detailed below, this leads to iterative schemes where
the update to the hyperparameter, β,
is computed in terms of the updated ensemble (for ),
and .
Some of the literature below is concerned with estimating  (jointly with ),
but is still pertinent by the proximity of the problem to that of estimating .


<cit.>
introduce the Variational Bayes (VB) method in the framework of the Kalman filter.
<cit.> extend it to the EnKF.
The VB method imposes an approximate posterior with
two factor distributions, (|) (|),
fitted by minimizing the
Kullback-Leibler divergence from the correct posterior, (,|).
This yields the condition that each factor be
the (geometric) marginal of (,|) with respect to the other,
i.e. two coupled  equations.
Independence is assumed for the prior:
(,) = (|,)(|^f, ν^f).


It is then shown that the distributions of the approximate posterior are again 
 and ,
with parameters computable by fixed point iteration.
Iteration i+1 consists of the EnKF equations
using the i-th estimate of ^a to compute
the analysis mean, ^a, and covariance, ,
whose updated values are then used to in the update
^a = {ν^f^f + }/ν^a
where ν^a = ν^f + 1 and

    =
    	( - ^a)( - ^a) +   .




<cit.>, also estimating  in the framework of the EnKF,
use the expectation maximization (EM) method to
maximize the marginal posterior (|), with an  prior.
The expectation is over , given  and the current estimate of .
The empirical distribution is assumed for the prior:
() ≈ N^-1∑_n δ(-_n),
yielding iterations involving weighted statistics.
Our investigation indicates that if, instead,
the standard EnKF assumption had been used:
() ≈(|,),
then the method would have yielded iterations as in <ref>.
This is unsurprising in view of the close connection between EM and VB.


An original result is obtained by applying the VB method to estimate the inflation parameter.
Consider the prior
(,β) = (|,β) (β|β^f,ν^f);
note that  and β are not assumed independent.
It can then be shown that the resulting VB scheme consists of using the i-th iterate of β^a
to compute the i+1 iterates of ^a and 
which, in turn, are used to compute the i+1 iterate
β^a = {ν^fβ^f + β̂}/ν^a
where ν^a = ν^f + M and

    β̂   =
    	^a-^f^2_ + (^-1)
    	  ,

which should be computed in the ensemble subspace if N ≤ M.
<Ref> may be interpreted using
the trigonometric relations of <cit.>.


	
	
	
	
	

Numerical experiments indicate that the bias of this VB method (with a flat prior)
is significantly higher than for β̂_, typically also with a larger variance.
It seems likely that the reason is similar to that of β̂_ML,
analysed in <ref>.







 §.§ Hierarchical methods


Instead of the variational approach,
a more principled approach is to represent
the marginal (β|) with a Monte-Carlo sample.
This yields a hyper-ensemble
of distributions (|β,) and their ensembles.
The label “hierarchical” is sometimes reserved for adaptive filters
treating the hyperparameter in this more Bayesian (i.e. full-) manner.
However, the number of realizations quickly becomes exorbitant.
On the other hand, if only a single sample is drawn from (|β,) for each β
(i.e. not using an ensemble of ensembles),
then it seems especially prone to spurious correlations.
Another criticism is that it might cause discrete jumps between localization domains.

<cit.> provided one example of a hierarchical EnKF,
but without conditioning the hyperparameter on ,
without accounting for model error,
and with a forecast of (the hyperprior of)  that decays towards 0.
<cit.> present an ambitious hierarchical EnKF,
explicitly considering both model and sampling error.
However, the filter is only tested in experiments with a novel, univariate model.

An intriguing approach is that of <cit.>.
They peg the scaling of  and  together,
so as to estimate only a single parameter, β.
This is difficult to justify, but simplifies the problem significantly
because then  also scales with β (provided linear dynamics)
so that (,β) = (|,β) (β|…)
is conjugate to (|,β),
and the full, joint posterior is available in closed form,
with trivial parametric updates,
hence avoiding Monte-Carlo.
However, without pegging  to , the conjugacy is lost.

The EnKF-N can also be said to be hierarchical
because of its careful treatment of the full marginal, (α|),
before the variational approximation eqn:post_EB.








§ MORE ON THE MARGINAL INFLATION ESTIMATORS



This section complements <ref>.




 §.§ Other trace-based estimators




























It is common to form chi-square diagnostics
by measuring  by its Mahalanobis norm
<cit.>.
This provides the motivation to use (1)^-1
to transform <ref>, yielding:

    β̂_ =
    	
    	
    	^2_(1) - ((1)^-1)/((1)^-1)  .


Alternatively, <ref> can be transformed
by ()^-1,
yielding
(  - ) ()^-1
=
β_P
For the purpose of inflation,
this seems like the best option because then the trace
consists of terms with the same expected magnitude,
yielding the lowest aggregate variance.
The estimator becomes:

    β̂_   =
    	1/P{^2_
    		-( ()^-1)
    		}  ,

If  is rank-deficient,
then β̂_ must be defined using the pseudo-inverse.
The estimate would then not be impacted by components of  outside of the ensemble subspace.


Lastly, note that none of the trace-based estimators are computationally costly,
because they can all be computed via the SVD eqn:svd_nRY.







 §.§ Single-cycle bias



This subsection considers
the properties of the estimators within a single analysis,
based only on the inflation estimators' sampling distributions
Here, the assumption of <ref> that N=∞ is undone,
so that  is also random (before conditioning on the ensemble), in addition to , and ,
and subject to sampling errors.
As will be shown, this
yields biases in the inflation estimators.




Now, using variables defined via the SVD eqn:svd_nRY,
it can be shown that each of the following estimators of β
satisfy the condition:

    0 =
    	∑_i=1^P γ_i
    	( 1 + β̂σ̅_i^2 - d_i^2 )
    	  ,

where

d_i is the i-th component of the transformed innovation,
,
and

    γ_i =
    	
    		1                                   for β̂_
    
    		(1+σ̅_i^2)^-1   for β̂_
    
    		(σ̅_i^2)^-1   for β̂_
    
    		
    			
    			
    		
    																		 
    		(1+σ̅_i^2)/(1+β̂σ̅_i^2)^2    for β̂_ML  .





<Ref> may be solved explicitly for β̂,
except in the case of β̂_ML.










Nevertheless, <ref> may be used to provide an insight on the bias of β̂_ML,
a subject of study since <cit.>.

Indeed, note that while the full matrix, , is an unbiased estimator of ,
the spectrum of  is a biased estimate of the spectrum of 
<cit.>.
Thus, the spectrum of ^-1,
namely {σ̅_i^2}_i=1^P, is also biased.
Hence, generally, functions of the spectrum will be biased.
An important exception is that
(∑_i σ̅_i^2) = (^-1),
meaning that the expectation of <ref> holds for β̂_.
Note that β̂_ is still biased, as it requires inverting ∑_i σ̅_i^2.
Nevertheless, considering the expressions eqn:s_conds for γ_i,
it seems logical that the bias of {σ̅_i^2}_i=1^P
will significantly carry over into the more complicated ones, such as β̂_ML.
Numerical experiments confirm this,
and show that the bias of β̂_ML is worse than it is for β̂_,
but less than for β̂_ and β̂_.


Why is the bias of β̂_ of <ref> the least?
Loosely speaking, because the trace of  is taken before dividing.
Moreover, in case  and  have the same structure,
i.e. ^-1/2 = σ^2 _P,
the bias of β̂_ can be obtained analytically, as follows.
Due to the modified assumption eqn:fawefhawe,
<ref> now yields ∼(/β,).
Thus, the diagonal elements of ^-1/2 are ,
with distribution (σ^2/β,N-1).
Taking the trace increases the certainty to ν = P:

    σ̅^2 = (^-1/2) ∼(σ^2/β,ν)
    	  .

Then, according to <ref> and <ref>,
[1/σ̅^2] = ν/ν-2β/σ^2.
Meanwhile,
[ ] =  + (1 + 1/β N)
so that
[ ^2_/P - 1 ] = (1 + 1/β N) σ^2.
Now, the sample mean and variance of Gaussian samples are independent.
Hence, the nominator and denominator of β̂_
of <ref> are independent, so that


    [β̂_]
    	   = [ ^2_/P - 1 ]  [1/σ̅^2]
    	
       = (1 + 1/β N) ν/ν-2β  ,

which is close to β for large N and P.






















































































































































































































































































 §.§ Acknowledgements

The authors thank
the two anonymous reviewers for their constructive comments;
A. Farchi for technical help and discussions;
L. Bertino for fostering the collaboration and for his patience;
C. Grudzien for insight on the behaviour of the Lorenz system;
F. Counillon, for reading and opinion;
A. Karspeck and M. Gharamti for stimulating questions.

Author P. N. Raanes has been funded by the EmblAUS project of the Nordic countries funding agency NordForsk,
and by DIGIRES, a project sponsored by PETROMAKS2 of the Research Council of Norway and industry partners.
CEREA is thanked for providing office space and community for visiting researcher P. N. Raanes. 
CEREA is a member of the Institut Pierre-Simon Laplace (IPSL).
Author A. Carrassi has been partly funded by the project REDDA of the Norwegian Research Council. 






wileyqj







