Asymptotic forecast uncertainty and the unstable subspace in the presence of additive model error
Colin Grudzien Alberto Carrassi Marc Bocquet
math.OC



Correlations between synapses in pairs of neurons slow down dynamics in randomly connected neural networks
    Srdjan Ostojic
    December 30, 2023
==========================================================================================================





It is well understood that dynamic instability is among the primary drivers of forecast uncertainty in chaotic, physical systems.  Data assimilation techniques have been designed to exploit this phenomena, reducing the effective dimension of the data assimilation problem to the directions of rapidly growing errors.  Recent mathematical work has, moreover, provided formal proofs of the central hypothesis of the Assimilation in the Unstable Subspace methodology of Anna Trevisan and her collaborators: for filters and smoothers in perfect, linear, Gaussian models, the distribution of forecast errors asymptotically conforms to the unstable-neutral subspace.  Specifically, the column span of the forecast and posterior error covariances asymptotically align with the span of backward Lyapunov vectors with non-negative exponents.  

Earlier mathematical studies have focused on perfect models, and this current work now explores the relationship between dynamical instability, the precision of observations and the evolution of forecast error in linear models with additive model error.  We prove bounds for the asymptotic uncertainty, explicitly relating the rate of dynamical expansion, model precision and observational accuracy.  Formalizing this relationship, we provide a novel, necessary criterion for the boundedness of forecast errors.  Furthermore, we numerically explore the relationship between observational design, dynamical instability and filter boundedness.  Additionally, we include a detailed introduction to the Multiplicative Ergodic Theorem and to the theory and construction of Lyapunov vectors. 

While forecast error in the stable subspace may not generically vanish, we show that even without filtering, uncertainty remains uniformly bounded due its dynamical dissipation.  
However, the continuous re-injection of uncertainty from model errors may be excited by transient instabilities in the stable modes of high variance, rendering forecast uncertainty impractically large.  In the context of ensemble data assimilation, this requires rectifying the rank of the ensemble-based gain to account for the growth of uncertainty beyond the unstable and neutral subspace, additionally correcting stable modes with frequent occurrences of positive local Lyapunov exponents that excite model errors.





 Kalman filter, data assimilation, model error, Lyapunov vectors, control theory




  93E11, 93C05, 93B07, 60G35, 15A03




§ INTRODUCTION


The seminal work of Lorenz <cit.> demonstrated that, even in deterministic systems, infinitesimal perturbations in initial conditions can rapidly lead to a long-term loss of predictability in chaotic, physical models.  In weather prediction, this understanding led to the transition from single-trajectory forecasts to ensemble-based, probabilistic forecasting <cit.>.  Historically, ensembles have been initialized in order to capture the spread of rapidly growing perturbations <cit.>.  Data assimilation methods have likewise been designed to capture this variability in the context of Bayesian and variational data assimilation schemes; see e.g., Carrassi et al. for a recent survey of data assimilation techniques in geosciences <cit.>.  The ensemble Kalman filter, particularly, has been shown to strongly reflect these dynamical instabilities <cit.>, and its performance depends significantly upon whether these rapidly growing errors are sufficiently observed and corrected.   

The Assimilation in the Unstable Subspace (AUS) methodology of Trevisan et. al. <cit.> has provided a robust, dynamical interpretation of these observed properties of the ensemble Kalman filter.  For deterministic, linear, Gaussian models, Trevisan et al. hypothesized that the asymptotic filter error concentrates in the span of the unstable-neutral backward Lyapunov vectors (BLVs), and this has recently been mathematically proven.  Gurumoorthy et. al. <cit.> demonstrated that the null space of the forecast error covariance matrices asymptotically contain the time varying subspace spanned by the stable BLVs.  This result was generalized by Bocquet et. al. <cit.>, proving the asymptotic equivalence of reduced rank initializations of the Kalman filter with the full rank Kalman filter: as the number of assimilations increases towards infinity, the covariance of the full rank Kalman filter converges to a sequence of low rank covariance matrices initialized only in the unstable-neutral BLVs.

The convergence of the Kalman smoother error covariances onto the span of the unstable-neutral BLVs, and stability of low rank initializations, was established by Bocquet & Carrassi <cit.>; this latter work also numerically extended this relationship to weakly nonlinear dynamics and ensemble-variational methods.  The works of Bocquet et al. <cit.> and Bocquet & Carrassi <cit.> relied upon the sufficient hypothesis that the span of the unstable and neutral BLVs remained uniformly-completely observed.  This hypothesis has recently been refined to a necessary and sufficient criterion for the exponential stability of continuous time filters, in perfect models, in terms of the detectability of the unstable-neutral subspace <cit.>.  

The present study is concerned with extending the limits of the results developed in deterministic dynamics (perfect models), now to the presence of stochastic model errors.  This manuscript and its sequel <cit.> seek to: (i) determine the extent to which stable dynamics confine the uncertainty in the sequential state estimation problem in models with additive noise, and (ii) to use these results to interpret the properties, and suggest design, of ensemble-based Kalman filters with model error.  This manuscript studies the asymptotic properties of the full rank, theoretical Kalman filter <cit.>, and the unfiltered errors in the stable BLVs.  The sequel <cit.> utilizes these results to interpret filter divergence for reduced rank, ensemble-based Kalman filters.

In <ref> we present a detailed introduction to the BLVs.  In section <ref> develop novel bounds on the forecast error covariance, describing the evolution of uncertainty as the growth of error, due to dynamic instability and model imprecision, with respect to the constraint of observations.  Together, the rate of dynamic instability and the observational precision form an inverse relationship which we use to characterize the boundedness of forecast errors.  In <ref> and <ref>, we prove a necessary criterion for filter boundedness in autonomous and time varying systems: the observational precision, relative to the background uncertainty, must be greater than the leading instability which forces the model error.  Our results derive from the bounds provided in <ref> for autonomous dynamics and <ref> for time varying systems.  An important consequence is that under generic assumptions, forecast errors in the span of the stable BLVs remain uniformly bounded independently of filtering.  Described in <ref>, this extends the intuition of AUS now to the presence of model errors: filters need only target corrections to the span of the unstable and neutral BLVs to maintain bounded errors.  

However, the intuition of AUS needs additional qualifications when interpreting the role of model errors in reduced rank filters. Unlike perfect models, uncertainty in the stable  BLVs does not generically converge to zero as a consequence of re-introducing model errors.  Moreover, while stability guarantees that unfiltered errors remain uniformly bounded in the stable BLVs, the uncertainty may still be impractically large due: even when a Lyapunov exponent is strictly negative, positive realizations of the local Lyapunov exponents can force transient instabilities which strongly amplify the forecast uncertainty.  The impact of stable modes on forecast uncertainty differs from similar results for nonlinear, perfect models by Ng. et. al. <cit.>, and Bocquet et. al. <cit.>, where the authors demonstrate the need to correct stable modes in the ensemble Kalman filter due to sampling errors induced by nonlinearity.  Likewise, this differs from the EKF-AUS-NL of Palatella & Trevisan <cit.>, that accounts for truncation errors in the estimate of the forecast uncertainty in nonlinear models.  In <ref>, we derive the mechanism for the transient instabilities amplifying perturbations as a linear effect in the presence of model errors.  We furthermore provide a computational framework to study the variance of these perturbations.    


In <ref>, we study the filter boundedness and stability criteria of Bocquet et al. <cit.> and Frank & Zhuk <cit.> in their relation to bounding forecast errors in imperfect models.  Likewise, we explore their differences in the context of dynamically selecting observations, similar to the work of Law et al. <cit.>.  With respect to several observational designs as benchmarks, we numerically demonstrate that the unconstrained growth of errors in the stable BLVs of high variance can be impractically large compared to the uncertainty of the full rank Kalman filter.  
These results have strong implications for ensemble-based filtering in geosciences and weather prediction, where ensemble sizes are typically extremely small relative to the model dimension.  In perfect models, an ensemble size large to correct the small number unstable and neutral modes might suffice.  However, our results suggest the need to further increase the rank of ensemble-based gains.  The significance of this result for ensemble-based Kalman filters and their divergence is further elaborated on in the sequel <cit.>.




§ LINEAR STATE ESTIMATION


The purpose of recursive data assimilation is estimating an unknown state with a sequential flow of partial and noisy observations; we make the simplifying assumption that the dynamical
and observational models are both linear and the error distributions are Gaussian.  In this setting, given a Gaussian distribution for the initial state, the distribution of the estimated state is Gaussian at all times.  Formulated as a Bayesian inference problem, we seek to estimate the distribution of the random vector _k ∈ℝ^n evolved via a linear Markov model,

    _k   = _k_k-1 + _k ,
 
with observations _k ∈ℝ^d given as

    _k   = _k _k + _k.

The model variables and observation vectors are related via
the linear observation operator _k: ℝ^n ↦ℝ^d.  Let _n denote the n × n identity matrix.  We denote the model propagator from time t_l-1 to time t_k as _k:l≜_k ⋯_l, where _k:k≜_n.

For all k, l ∈ℕ, the random vectors of model and observation noise, _k, _l∈ℝ^n and _k,_l∈ℝ^d, are assumed mutually independent, unbiased, Gaussian white
sequences.  Particularly, we define

    𝔼[_k_l^] = δ_k,l_k    and  𝔼[_k_l^] = δ_k,l_k,

where 𝔼 is the expectation, _k∈ℝ^d× d is the observation error covariance matrix at time t_k, and _k ∈ℝ^n× n stands for the model error covariance matrix. The error covariance matrix _k can be
assumed invertible without losing generality.  For simplicity we assume the dimension of the observations d≤ n will be fixed.  

For two positive semi-definite matrices,  and , the partial ordering is defined
≤ if and only if - is positive semi-definite. To avoid pathologies, we assume that the model error and the observational error covariance matrices are uniformly bounded, i.e., there are constants q_inf, q_sup, r_inf, r_sup∈ℝ such that for all k,

    ≤ q_inf_n ≤_k ≤ q_sup_n,
       < r_inf_d ≤_k ≤ r_sup_d.



Rather than explicitly computing the evolution of the distribution for _k, the Kalman filter computes the forecast and posterior distributions parametrically via recursive equations for the mean and covariance of each distribution. 

The forecast error covariance matrix _k of the Kalman filter satisfies the discrete-time
dynamic Riccati equation <cit.>

    _k+1 = _k+1_n+_k _k^-1_k_k+1^ +_k+1 ,

where _k ≜_k^_k^-1_k is the precision matrix of the observations.  


Equation <ref> expresses the error covariance matrix, _k+1, as the result of a two-step
process: (i) the assimilation at time t_k yielding the analysis error covariance,

    ^a_k = _n+_k _k^-1_k;

and (ii) the forecast, where the analysis error covariance is forward propagated by

    _k+1 = _k+1^ a_k_k+1^ + _k+1 .

Assuming that the filter is unbiased, such that the initial error is mean zero, it is easy to demonstrate that the forecast and analysis error distributions are mean zero at all times.  In this context, the covariances _k,^a_k represent the uncertainty of the state estimate defined by the filter mean.  As we will focus on the evolution of the covariances, we neglect the update equations for the mean state and refer the reader to Jazwinski <cit.> for a more complete discussion.


The classical conditions for the 
boundedness of filter errors, and the independence of the asymptotic filter behavior from its initialization, are given in terms of observability and controllability.  Observability is the condition that given finitely many observations, the initial state of the system can be reconstructed.
Controllability describes the ability to move the system from any
initial state to a desired state given a finite sequence of control actions — in our case the moves are the realizations of model error.  These conditions are described in the following definitions, beginning with the information and controllability matrices.

We define _k:j to be the time varying information matrix and _k:j to be the time varying controllability matrix, where

    _k:j   ≜∑_l=j^k_k:l^-_l _k:l^-1,       _k:j≜∑_l=j^k_k:l_l _k:l^.

For γ≥ 0 let us define the weighted controllability matrix as

    ^γ_k:j≜∑_l=j^k 1/1+γ^k-l_k:l_l _k:l^.


Note that, ^0_k:j≡_k:j.  We recall from section 7.5 of Jazwinski <cit.> the definitions of uniform complete observability (respectively controllability).

Suppose there exists N_Φ, a,b > 0 independent of k such that k> N_Φ implies

    0< a _n ≤_k:k - N_Φ≤ b _n,

then the system is uniformly completely observable.  Likewise suppose there exists N_Υ, a,b >0 independent of k for which k> N_Υ implies

    0< a _n ≤_k:k- N_Υ≤ b _n,

then the system is uniformly completely controllable.



Assume that the system of equations <ref> and <ref> is 


  * uniformly completely observable; 

  * uniformly completely controllable.
  


We will explicitly refer to <ref> whenever it is used.  When we refer <ref> alone, we refer to both parts (a) and (b).  At times, we will explicitly only use either part (a) or (b) of <ref>.


Suppose the system of equations <ref> and <ref> satisfies <ref> and _0 >0.  Then there exists constants p^a_inf and p^a_sup independent of k such that the analysis error covariance is uniformly bounded above and below,

    0 < p^a_inf_n ≤^a_k ≤ p^a_sup_n < ∞.
 
Given any two initializations of the prior error covariance _0,_0>0, with associated sequences of analysis error covariances ^a_k,^a_k, the covariance sequences converge, lim_k→∞ł|^a_k - ^a_k |̊ = 0,  exponentially in k.

These are classical results of filter stability, see for example Theorem 7.4 of Jazwinski <cit.>, or Bougerol's work with random matrices <cit.> for a generalization.    

The square root Kalman filter is a reformulation of the recurrence in equation <ref> which is used to reduce computational cost and obtain superior numerical precision and stability over the standard implementations see, e.g., <cit.>.  The advantage of this formulation to be used in our analysis is to explicitly represent the recurrence in equation <ref> in terms of positive semi-definite, symmetric matrices.


Let _k be a solution to the time varying Riccati equation <ref> and define _k ∈ℝ^n× n to be a Cholesky factor of _k, such that

    _k = _k ^_k.
  


The root _k in equation <ref> can be interpreted as an ensemble of anomalies about the mean as in the ensemble Kalman filter <cit.>.  In operational conditions, it is standard that the forecast error distribution is approximated with a sub-optimal, reduced rank surrogate <cit.>.  Using a reduced rank approximation, the estimated covariance and exact error covariance are not equal, and this can lead to the systematic underestimation of the uncertainty <cit.>.  However, in the following we will assume that _k is computed as an exact root.  The sequel to this work explicitly treats the case of reduced rank, sub-optimal filters <cit.>.


We order singular values σ_1 > ⋯ > σ_n such that,

    0 ≤σ_n_k^_k _k_n ≤_k^_k _k ≤σ_1_k^_k _k_n <∞.

We define

    α≜inf_k {σ_n_k^_k _k}≥ 0,       β≜sup_k {σ_1_k^_k _k}≤∞,

and we write 0≤α_n ≤_k^_k _k ≤β_n ≤∞ for all k.


Equation <ref> is closely related to the singular value analysis of the precision matrix by Johnson et al. <cit.> and the analysis of the conditioning number for the Hessian of the variational cost function by Haben et al. <cit.> and Tabeart et al. <cit.>.  These works study the information gain from observations, relative to the background uncertainty, due to the assimilation step.  The primary difference between these earlier works and our study here is that the background error covariance is static in these variational formulations, while in the present study the root _k is flow dependent.  In this flow dependent context, the constant α (respectively β) is interpreted as the minimal (respectively maximal) observational precision relative to the maximal (respectively minimal) background forecast uncertainty.  The constant α is nonzero if and only if the principal angles between the column span of _k and the kernel of _k are bounded uniformly below.  Generally, we thus take α=0 unless observations are full dimensional.  A nonzero value for α can be understood as an ideal scenario.

Using <ref> and the matrix shift lemma <cit.> we re-write the forecast Riccati equation <ref> as 

    _k    = _k(_n+ _k-1_k-1)^-1_k-1_k^ + _k 
       = _k _k-1(_n + ^_k-1_k-1_k-1)^-1^_k-1_k^ + _k

from which we infer

    1/1+ β_k _k-1_k^ + _k ≤_k ≤1/1+ α_k _k-1_k^ + _k.

Iterating on the above inequality, we obtain the recursive bound

    1/1+ β^k _k:0_0_k:0^ + ^β_k:1≤_k ≤1/1+ α^k _k:0_0_k:0^ + ^α_k:1.


Equation <ref> holds if there is no filtering step, setting β = α =0.
 

The bounds in equation <ref> explicitly describe the previously introduced uncertainty as dynamically evolved to time k, relative to the constraint of the observations.  We will utilize the BLVs vectors to extract the dynamic information from the sequences of matrices _k:l, _k:l^.



§ LYAPUNOV VECTORS



This section contains a short introduction to Lyapunov vectors and the Multiplicative Ergodic Theorem (MET).  For a more comprehensive introduction, there are many excellent resources at different levels of complexity, see for example <cit.>.  There is inconsistent use of the terminology for Lyapunov vectors in the literature, so we choose to use the nomenclature of Kuptsov & Parlitz <cit.> for its accessibility and self-consistency.

Consider the growth or decay of an arbitrary, norm one vector v_0 ∈ℝ^n to its state at time t_k via the propagator _k:0.  This is written as 

    ł| v_k |̊ =  ł| _k:0 v_0 |̊ =  √(   v_0^_k:0^_k:0 v_0  ),

so that the eigenvectors of the matrix _k:0^_k:0 describe the principal axes of the ellipsoid defined by the unit disk evolved to time t_k.  Using the above relationship for the reverse time model, we see the growth or decay of the unit disk in reverse time as

    ł| u_-k|̊ =  ł| ^-1_0:-k u_0 |̊ = √(   u_0^_0:-k^-^-1_0:-k u_0  ).

The principal axes of the past ellipsoid that evolves to the unit disk at the present time are thus precisely the eigenvectors of the matrix ^-_0:-k^-1_0:-k.
There is no guarantee in general that there is consistency between the asymptotic forward and reverse time growth and decay rates, i.e., in equations 
<ref> and <ref> as k →∞.  Generally, models may have Lyapunov spectrum defined as intervals of lower and upper growth rates, see e.g., Dieci & Van Vleck <cit.>.  However, as we are motivated by the tangent-linear model for a nonlinear system, we may assume some “regularity” in the dynamics.

The anti-symmetry of the forward/reverse time, regular and adjoint models' growth and decay is known as Lyapunov-Perron regularity (LP-regularity) <cit.>, and is equivalent to the classical Oseledec decomposition <cit.>[see Theorem 2.1.1].  LP-regularity guarantees that: (i) the Lyapunov exponents are well defined for the linear model as point-spectrum; (ii) the linear space is decomposable into subspaces that evolve covariantly with the linear propagator; and (iii) each such subspace asymptotically grows or decays according to one of the point-spectrum rates.  We summarize the essential results of Oseledec's theorem for use in our work in the following, see Theorem 2.1.1 of Barreira & Pesin <cit.> for a complete statement.

[Oseledec's Theorem]
The model _k = _k _k-1 is LP-regular if and only if there exists real numbers λ_1 > ⋯> λ_p, for 1 ≤ p ≤ n, and subspaces _k^i ⊂ℝ^n, ^i_k = κ_i, such that for every k,l ∈ℤ

    ⊕_i=1^p ^i_k = ℝ^n       _k ± l:k ^i_k = ^i_k ± l,

and ∈^i_k implies
    
    lim_l→∞1/ l log(ł| _k ± l:k|̊) = ±λ_i .



For p≤ n, the Lyapunov spectrum of the system <ref> is defined as the set 
{λ_i : κ_i }_i=1^p
where λ_1 > ⋯ > λ_p and κ_i corresponds to the multiplicity (degeneracy) of the exponent λ_i.
We separate non-negative and negative exponents, λ_n_0≥ 0 > λ_n_0 +1, such that each index i> n_0 corresponds to a stable exponent.  The subspaces ^i_k are denoted Oseledec spaces, and the decomposition of the model space into the direct sum is denoted Oseledec splitting.



 For arbitrary linear systems LP-regularity is not a generic property — it is the MET that shows that this is a typical scenario for a wide class of nonlinear systems.  A point will be defined to be LP-regular if the tangent-linear model along its evolution is LP-regular.  We state a classical version of the MET <cit.> but note that there are more general formulations of this result and more general forms of the associated covariant-subspace decompositions.  These results go beyond the current work, see e.g., Froyland et al. <cit.> and Dieci et al. <cit.> for a stronger version of the MET and related topics.

[Multiplicative Ergodic Theorem]
If f is a C^1 diffeomorphism of a compact, smooth, Riemannian manifold M, the set of points in M which are LP-regular has measure 1 with respect to any f-invariant Borel probability measure ν on M.  If ν is ergodic, then the Lyapunov spectrum is constant with ν-probability 1.

Loosely, the MET states that, with respect to an ergodic probability measure (that is compatible with the map f and the usual topology), there is probability one of choosing initial conditions for which the Lyapunov exponents are well defined and independent of initial condition.  This form of the MET has a wide range of applications in differentiable dynamical systems, but the MET is not limited to this setting.  The strong version of the MET has been applied in, e.g., hard disk systems, the truncated Fourier expansions of PDEs, and with non-autonomous ODEs and their transfer operators <cit.>.  For the rest of this work, we will take the hypothesis that our model satisfies LP-regularity.

The model defined by the deterministic equation 

    _k = _k _k-1,

is assumed to be LP-regular.

The deterministic evolution in equation <ref> comes naturally in the formulation of the Kalman filter, where the mean state is evolved via the deterministic component in the forecast step.  For Gaussian error distributions, the evolution of the forecast error distribution is interpreted in terms of Oseledec's theorem as the evolution of deviations from the mean, propagated via the equations for perturbations.  While Oseledec's theorem guarantees that a decomposition of the model space exists, constructing such a decomposition is non-trivial.  Motivated by equations <ref> and <ref>, we define the following operators as in equations (13) and (14) of Kuptsov & Parlitz <cit.>.

We define the far-future operator as

    ^+(k) ≜   lim_l→∞
        ^_k+l: k_ k+l: k
    ^1/2l,

and the far-past operator as

    ^-(k) ≜   lim_l→∞
        ^-_k: k-l^-1_ k: k-l
    ^1/2l.



In the classical proof of the MET, the far-future/past operators are shown to be well defined positive definite, symmetric operators <cit.>.  As they are diagonalizable over ℝ, we order the eigenvalues of ^+(k) as μ^+_1(k)>⋯ >μ^+_p(k) and the eigenvalues of ^-(k) as μ^-_p(k) > ⋯ > μ^-_1(k).  
By the MET, the eigenvalues μ^±_i(k) are independent of k and satisfy the relationship

    log(μ^±_i) = ±λ_i.



Let the columns of the matrix _k, respectively _k, be any orthonormal eigenbasis for the far-future operator ^+(k), respectively far past operator ^-(k).  Order the columns block-wise, such that for each i =1,⋯, p and each j =1,⋯, κ_i, ^i_j_k is an eigenvector for μ^+_i and ^i_j_k is an eigenvector for μ^-_i.  We define ^i_j_k to be the i_j-th forward Lyapunov vector (FLV) at time k and _k^i_j to be the i_j-th backward Lyapunov vector (BLV). Let the columns of _k form any basis such that for each i=1,⋯, p and each j=1,⋯, κ_i, ^i_j_k ∈^i_k.  Then we define ^i_j_k to be the i-th covariant Lyapunov vector (CLV) at time k.


The CLVs are defined only by the Oseledec spaces, and therefore, are independent of the choice of a norm — any choice of basis subordinate to the Oseledec splitting is valid.  On the other hand, the FLVs and the BLVs are determined specifically with respect to a choice of a norm and the induced metric.  The choice of basis in each case can be made uniquely (up to a scalar and the choice of a norm) only when p=n.  For the remaining work we will focus on the BLVs; for a general survey on constructing FLVs, BLVs and CLVs, see e.g., Kuptsov & Parlitz <cit.> and Froyland et. al. <cit.>.

The Oseledec spaces and Lyapunov vectors can also be defined in terms of filtrations, i.e., chains of ascending or descending subspaces of ℝ^n.  This forms an axiomatic approach to constructing abstract Lyapunov exponents used by, e.g., Barreira & Pesin <cit.>.  The BLVs describe an orthonormal basis for the ascending chain of Oseledec subspaces, the backward filtration <cit.>.  For all 1 ≤ m≤ p we obtain the equality

    ⊕_i=1^m ^i_k = ⊕_i=1^m {_k^i_j}_j=1^κ_i,

by equation (17) by Kuptsov & Parlitz <cit.>, and the decomposition of the backward filtration in equations (1.5.1) and (1.5.2) of Barreira & Pesin <cit.>.  Note that equation <ref> does not imply ^i_j_k ∈^i_k for i> 1, as the BLVs are not themselves covariant with the model dynamics.  However, the BLVs are covariant with the QR algorithm.     

Outside of a set of Lebesgue measure zero, a choice of i_j linearly independent initial conditions for the recursive QR algorithm converges to some choice for the leading i_j BLVs.  For any k, the BLVs satisfy the relationship 

    _k _k-1 =  _k _k, ⇔_k = _k _k _k-1^

where _k is an upper triangular matrix. Moreover, for any i_j and any k,

    lim_l → -∞1/k-llogł|^_k:l^ i_j_k|̊ = λ_i.



The covariance of the BLVs with respect to the QR algorithm in equation <ref> can be derived from equations <ref> and <ref>.  For all 1≤ m ≤ p,

    _k ⊕_i=1^m {_k-1^i_j}_j=1^κ_i = ⊕_i=1^m {_k^i_j}_j=1^κ_i,

due to the covariance of the Oseledec spaces.  Therefore the transformation _k represented in a moving frame of BLVs is upper triangular.  When the spectrum is degenerate, p<n, there is non-uniqueness in the choice of the BLVs.  However, given an initial choice of the BLVs at some time k-1, the choice of BLVs at time k can be defined directly via the relationship in <ref>.  This is the relationship derived in equation (31) by Kuptsov & Parlitz <cit.>, and is the basis of the recursive QR algorithms of Shimada & Nagashima <cit.> and Benettin et. al. <cit.>.  A choice of BLVs gives a special choice of the classical Perron transformation <cit.>[see Theorems 3.3.1 & 3.3.2], and in particular, it is proven by Ershov & Potapov <cit.> that outside of a set of Lebesgue measure zero, the recursive QR algorithm converges to some choice of BLVs.

Note that the far-future/past operators are  also well defined for the propagator of the adjoint model _k = ^-_k _k-1.
Equation <ref> thus follows from the far-past operator for the adjoint model, defined

    ^∗ -(k) ≜lim_l →∞
        _k:k-l_k:k-l^
    ^1/2l.

It is easy to verify that the BLVs defined by the adjoint model agree with those defined via the regular model — in each case, the left singular vectors of _k:k-l converge to a choice of the BLVs as l→∞.  Notice that the eigenvalues of ^∗ -(k) are reciprocal to those of ^-(k), i.e.,
μ^∗ -_i = 1/μ^-_i.  Thus by equation <ref>, logμ^∗ -_i =  λ_i.


Equation <ref> describes the dynamics in the moving frame of BLVs, where the transition map from the frame at time t_k-1 to time t_k is given by _k.  Applying the change of basis sequentially for the matrix _k:l, we recover

    _k =  _k _k _k-1^⇒_k:l =  _k _k:l_l^,

where we define _k:l≜_k ⋯_l.  We note that _k:k=_n implies _k:k≡_n.  Let _i_j denote the i_j-th standard basis vector, such that

    ł|^_k:l^i_j_k |̊^2  
    = _i_j^_k:l_k:l^_i_j 
    =ł|^_k:l^i_j|̊^2

where ^_k:l^i_j denotes the i_j-th column of ^_k:k-l, i.e., the i_j-th row of _k:k-l.  For any k and any ϵ>0, there exists some N_ϵ, k such that if k-l is taken sufficiently large, equation <ref> guarantees

    e^2(λ_i - ϵ)l≤ł|^_k:k-l^i_j|̊^2 ≤ e^2(λ_i + ϵ)l.
 



For each k>l, i=1,⋯, p and j=1, ⋯, κ_i we define the i_j-th local Lyapunov exponent (LLE) from k to l as 1/k-llog| T^i_j_k:l| where T^i_j_k:l is defined to be the i_j-th diagonal entry of _k:l.


For any fixed l,

    lim_k→∞1/k-llog| T^i_j_k:l|   = λ_i



This is also discussed by Ershov & Potapov <cit.>, in demonstrating the convergence of the recursive QR algorithm.  For a discussion on the numerical stability and convergence see, e.g., Dieci & Van Vleck <cit.>.


Perturbations of model error to the mean equation for the Kalman filter are not governed by the asymptotic rates of growth or decay, but rather, the LLEs.
While the LLE 1/k-llog| T^i_j_k:l| approaches the value λ_i as k-l approaches infinity, 
its behavior on short time scales can be highly variable.  Particularly, for an arbitrary LP-regular system, the rate of convergence in equation <ref> may depend on k.  An important class of such systems is, e.g., non-uniformly hyperbolic systems <cit.>[see chapter 2].  To make the LLEs tractable, we make an additional assumption, compatible with the typical assumptions for partial hyperbolicity <cit.>.  We adapt the definition of partial hyperbolicity from Hasselblatt & Pesin <cit.> to our setting.


Let λ_n_0= 0.  For every k we define the splitting into unstable, neutral and stable subspaces: 

    E^u_k ≜⊕_i= 1^n_0 - 1^i_k,        E^c_k ≜^n_0_k       and        E^s_k ≜⊕_i= n_0 + 1^p^i_k.

Suppose there exists constants C>0 and 

    0 < η_s ≤ν_s < η_c ≤ν_c <  η_u ≤ν_u

independent of k such that ν_s <1 < η_u and for any l>0, ∈ E^m_k, ł| |̊=1, and m∈{s,c,u}

    (η_m)^l/C≤ł| _k+l:k|̊≤ (ν_m)^l C .

Then the model <ref> is (uniformly) partially hyperbolic (in the narrow sense).


Partially hyperbolic systems, as in <ref>, have LLEs which are bounded uniformly with respect to rates defined on the subspaces in equation <ref>.  When C is taken large the definition permits transient growth of stable modes and transient decay of unstable modes.  The neutral subspace encapsulates diverse behaviors which always fall below prescribed rates of exponential growth or decay.  We will make a slightly stronger assumption on these uniform growth and decay rates that is equivalent to fixing a uniform window of transient variability on each Lyapunov exponent.

Let ϵ>0 be given.  We assume that for each i there exists some N_i,ϵ, independent of k and j, such that for any ^i_j_k whenever k-l> N_i,ϵ

    -ϵ < 1/k-llogł| ^_k : l^i_j_k |̊  - λ_i < ϵ ,

i.e., the growth and decay is uniform (translation invariant) in k.

Unless specifically stated otherwise, we assume <ref> for the remaining of this paper.  However, our results may be generalized to all systems satisfying <ref> by using only the uniform rates of growth or decay on the entire unstable, neutral and stable subspaces in equation <ref>.  Our results also apply to systems without neutral exponents, i.e. λ_n_0>0, as a trivial extension.



§ DYNAMICALLY INDUCED BOUNDS FOR THE RICCATI EQUATION




 §.§ Autonomous systems

Consider the classical theorem regarding the existence and uniqueness of solutions to the stable Riccati equation for autonomous dynamics.  This is paraphrased from Theorem 2.38, Chapter 7, of Kumar & Varaiya <cit.> in terms of the forecast error covariance recurrence in equation <ref>. 

The autonomous system is defined such that for every k

    _k ≡ ,     _k ≡ ,     _k ≡ ,     _k ≡   
    and    _k ≡.

Let = ^ for some ∈ℝ^n× n, the stable Riccati equation is defined as
      
    = (_n + ^ )^-1^^ +



Let the autonomous system defined by equations <ref>, <ref> and <ref> satisfy <ref>.  There is a positive semi-definite matrix, ≡^, which is the unique solution to the stable Riccati equation <ref>.  For any initial choice of _0, if _k satisfies the recursion in equation <ref>, then lim_k→∞_k =.
   

Slightly abusing notation, take α and β to be defined by the solution to the stable Riccati equation <ref>,

    α≜σ_n^≥ 0       β≜σ_1^ < ∞.

Then for any k we recover the invariant recursion for the stable limit

    (1+ β)^-k^k^^k + ^β_k:1≤≤ (1+ α)^-k^k^^k + ^α_k:1.



Assume equations <ref> and <ref> satisfy <ref> and define α,β for the stable Riccati equation as in equation <ref>.  For any 1≤ i≤ p, if there exists ϵ > 0 such that 

    e^2(λ_i + ϵ)/1+ α < 1,

choose N_i,ϵ as in <ref>. For the eigenvalue μ_i of ^, where |μ_i | = e^λ_i, choose any eigenvector _i_j.  Then  

    _i_j^_i_j   ≤_i_j^_i_j/1-e^2λ_i/1+α.

  Moreover, if ^i_j is the i_j-th BLV, then

    ^i_j^^ i_j   ≤^ i_j^^α_N_i,ϵ:0^ i_j +  e^2(λ_i +ϵ)/1+ α^N_i,ϵ +1q_sup/1-e^2(λ_i+ϵ)/1+α.

 
For every 1≤ i≤ p, any ϵ>0, and associated N_i,ϵ as in <ref>,

    _i_j^_i_j/1-e^2λ_i/1+β≤_i_j^_i_j

and 

    ^ i_j^^β_N_i,ϵ:0^ i_j +  e^2(λ_i - ϵ)/1+ β^N_i,ϵ +1q_inf/1-e^2(λ_i-ϵ)/1+β≤^ i_j^^ i_j .




Note that time invariant propagators trivially satisfy <ref> and it is easy to verify the relationship |μ_i | = e^λ_i directly from the definition of the Lyapunov exponents.  We begin by proving equations <ref> and <ref> for eigenvectors of ^.  If _i_j is an eigenvector of ^ associated to μ_i, equation <ref> implies

    _i_j^_i_j≤|μ_i|^2 /1+ α^k+1_i_j^_i_j + ∑_l=0^k|μ_i |^2 /1+α^l _i_j^_i_j

for every k.  For λ_i < 0 generally, or for any λ_i such that α > e^2λ_i - 1,

    lim_k→∞
        |μ_i|^2 /1+ α^k+1_i_j^_i_j +∑_l=0^k|μ_i |^2 /1+α^l _i_j^_i_j
      = _i_j^_i_j/1-|μ_i |^2 /1+α
 
 and

    _i_j^_i_j≤_i_j^_i_j/1-e^2λ_i/1+α .


The stable Riccati equation <ref> implies ≤.  Therefore, using the left side of <ref> demonstrates that for any eigenvector _i_j

    ∑^k_l=0|μ_i |^2/1+ β^l _i_j^_i_j≤_i_j^_i_j

for all k.  In particular, for every eigenvector _i_j we obtain

    _i_j^_i_j/1 - |μ_i|^2/1 + β   ≤_i_j^_i_j .

  
The above argument does not have a straightforward extension to the generalized eigenspaces so we coarsen the bound to obtain a closed limiting form in terms of the BLVs which retain the important growth characteristics under ^.  For i> n_0, or for any λ_i such that α > e^2λ_i - 1, there is a choice of ϵ as in equation <ref> and N_i,ϵ as in <ref>.  Let ≤p_sup_n, then from the right side of equation <ref> we derive

    ≤p_sup^k+1^^k+1/(1+ α)^k+1  + ∑_l=0^k^l^^l/(1+α)^l
       ≤p_sup^k+1^^k+1/(1+ α)^k+1  + ^α_N_ϵ,i:1 + q_sup∑_l=N_ϵ,i+1^k^l^^l/(1+α)^l,

which implies ^ i_j^^ i_j can be bounded above by

    p_supł| ^^k+1^ i_j|̊^2/(1+ α)^k+1 + ^ i_j^^α_N_i,ϵ:1^ i_j + q_sup∑_l=N_i,ϵ + 1^kł| ^^l ^ i_j|̊^2/(1+α)^l.

Utilizing equation <ref> we bound ^ i_j^^ i_j by

    p_supe^2(λ_i +ϵ)/1+ α^k+1  + ^ i_j^^α_N_i,ϵ:1^ i_j + q_sup∑_l=N_i,ϵ +1^ke^2(λ_i+ϵ)/1+α^l

for every k > N_i,ϵ.  Taking the limit of equation <ref> as k→∞ yields

    ^ i_j^^ i_j   ≤^i_j^^α_N_i,ϵ:1^ i_j +   e^2(λ_i +ϵ)/1+ α^N_i,ϵ +1q_sup/1-e^2(λ_i+ϵ)/1+α,

The lower bound is demonstrated by similar arguments with the lower bound in equation <ref>, utilizing the property < ∞.


<ref> is similar results in perfect models <cit.>, but with some key differences.  Once again that the estimation errors are dissipated by the dynamics in the span of the stable BLVs, but the recurrent injection of model error prevents the total collapse of the covariance to the unstable-neutral subspace.  In equation <ref>, we see that for very strong decay, when e^2λ_i≈ 0, or high precision observations, i.e., when the system is fully observed and as α→∞, the stable limit of the forecast uncertainty reduces to what is introduced by the recurrent injection of model error.  The SEEK filter of Pham et. al. <cit.> has exploited these properties by neglecting corrections in the stable eigenspaces and only making corrections in the unstable directions.  This is likewise the motivation for AUS of Trevisan et. al. <cit.>, though the work of AUS was concerned with nonlinear, perfect models.  

The upper bounds in equations <ref> and <ref> generally hold for i≤ n_0 only when the system is fully observed.  Therefore, these bounds can be considered an ideal bound for the unstable-neutral modes.  However, the lower bound in equation <ref> hold generally for i<n_0.  By assuming the existence of an invariant solution to the stable Riccati equation <ref>, we will recover a necessary condition for its existence.


Assume there exists a solution  to the stable Riccati equation <ref>.  Choose the smallest index i such that 1≤ i ≤ n_0 and there exists some generalized eigenvector _i_j of ^ for which _i_j∉null.
Then it is necessary that

    e^2λ_i/1+σ_1^2^-1/2 <1.



Let _i_1 be an eigenvector for ^ and _i_1≠ 0.  Then by the definition of β in equation <ref>, the equation <ref> holds for all k if and only if equation <ref> holds.  More generally, suppose {_i_j}_j=1^κ_i are (possibly complex) generalized eigenvectors forming a Jordan block for ^.  Let j be the smallest index for which _i_j≠ 0.  Recall that for each j∈{1,⋯,κ_i} the Jordan basis satisfies

    ^ -μ_i _n_i_j = _i_j-1

where _i_0≡ 0.  Therefore, for any m≥ 1, the vector ^ -μ_i _n^m _i_j is in the span of {_i_1, ⋯, _i_j-1}.   Let us define ≜^ -μ_i _n so that

    ∑_l=0^k+1^^l _i_j   = ∑_l=0^k+1+ μ_i _n^l  _i_j
        = ∑_l=0^k+1∑_m=0^l μ_i^l-ml  m^m    _i_j
        = ∑_l=0^k+1μ_i^l _i_j.

Multiply equation <ref> on the left with _i_j^$̋ (the conjugate transpose) and the right with_i_j.  Combining this with the equality in equation <ref>, proves the result. 
<ref> shows that it is necessary for the existence of the stable Riccati equation that observations are precise enough, relative to the background uncertainty, to counteract the strongest dynamic instability forcing the model error.  The quantity in <ref> thus represents the stabilizing effect of the observations, similar to the bounds on the conditioning number provided by Haben et al. <cit.> and Tabeart et al. <cit.>, but in <ref> expressly in response to the system's dynamic instabilities. 



 §.§ Time varying systems

In the following, we will extend the results of <ref> and <ref> to time-varying systems, and derive a uniform bound on the unfiltered errors in the stable BLVs in <ref>.


Assume equations <ref> and <ref> satisfy <ref> (b).  For any 1 ≤ i ≤ p, if there exists ϵ > 0 such that 

    e^2(λ_i + ϵ)/1+ α < 1,

choose N_i,ϵ as in <ref>.  Then there exists a constant 0≤ C_α,N_i,ϵ such that

    lim sup_k→∞^ i_j_k^_k^ i_j_k    ≤  C_α, N_iϵ +  e^2(λ_i +ϵ)/1+ α^N_i,ϵ +1q_sup/1-e^2(λ_i+ϵ)/1+α.

 
If <ref> (a) is also satisfied, then for every 1 ≤ i ≤ p, any ϵ>0 and associated N_i,ϵ, there exists 0≤ C_β,N_i,ϵ such that

    C_β, N_iϵ   +  e^2(λ_i - ϵ)/1+ β^N_i,ϵ +1q_inf/1-e^2(λ_i-ϵ)/1+β≤lim inf_k→∞^ i_j_k^_k^ i_j_k .

If the system satisfies <ref> (b) then

    ^α_k:k-N_i,ϵ   ≤^0_k:k-N_i,ϵ≡_k:k-N_i,ϵ≤ b_N_i,ϵ_n,

where b_N_i,ϵ is independent of k.  Therefore, there exists a constant depending on α and N_i,ϵ, but independent of k, such that

    ^α_k:k-N_i,ϵ≤ C_α,N_i,ϵ_n.
 
Let _0 ≤ p_0 _n bound the prior covariance.  Equation <ref> implies

    _k ≤ p_0 _k:0^_k:0/(1+α)^k +  q_sup∑_l=1^k_k:l^_k:l/(1+α)^k-l.

From the above, we bound ^ i_j_k^_k ^i_j_k with
 
    p_0 ł| ^_k:0^ i_j_k|̊^2/(1+α)^k +  ^ i_j_k^^α_k:k-N_i,ϵ^ i_j_k + q_sup∑_l=0^k-N_i,ϵ-1ł| ^_k:l^i_j_k|̊^2/(1+α)^k-l,

 thus
 
    ^ i_j_k^_k ^i_j_k ≤ p_0 e^2(λ_i + ϵ)/1+α^k +  C_α,N_i,ϵ + q_sup∑_l=N_i,ϵ+1^ke^2(λ_i + ϵ)/1+α^l .

Taking the lim sup in equation <ref> as k →∞ yields equation <ref>.  

Suppose that <ref> (a) and (b) are both satisfied, then by <ref> there exists a uniform bound on _k such that _k must also be uniformly bounded; together with uniform boundedness of _k and _k, this implies
β < ∞.  Note that

    ^β_k:k-N_i,ϵ   ≥1/1 + β^N_i,ϵ^0_k:k-N_i,ϵ≥1/1 + β^N_i,ϵ a_N_i,ϵ_n

for some constant a_N_i,ϵ independent of k.  This implies

    ^β_k:k-N_i,ϵ≥ C_β,N_i,ϵ_n

for a constant C_β,N_i,ϵ depending on β and N_i,ϵ but independent of k.  Utilizing the recursion in equation <ref>, choosing ϵ and an appropriate N_i,ϵ, and finally bounding the weighted controllability matrix with equation <ref> allows one to recover the lower bound in equation <ref> in a similar manner to the upper bound.
 

The above proposition shows that there is a uniform upper and lower bound on the forecast error for the Kalman filter, in the presence of model error, which can be described in terms of inverse, competing factors: the constantα(respectivelyβ) is interpreted as the minimal (respectively maximal) observational precision relative to the maximal (respectively minimal) background forecast uncertainty, represented in the observation variables.
  AdditionallyC_β, N_iϵ, C_α, N_iϵrepresent the lower and upper bounds on local variability of the evolution of model errors, before perturbations adhere within anϵthreshold to their asymptotic behavior.
 

Assume equations <ref> and <ref> satisfy <ref> (b), and there exists uniform bound to the forecast error Riccati equation <ref> for all k.  Then it is necessary that
    e^2λ_1/1+sup_kσ_1^2_k^-1/2_k _k <1.

If the forecast error Riccati equation <ref> is uniformly bounded, there is a 0<p_sup<∞ such that we have the inequality, _k≤ p_sup_n for all k, and β< ∞.  Using the lower bound in equation <ref>, for all k we have

    1/1 + β^k _k:0_0 ^_k:0 + ^β_k:1   ≤ p_sup_n.

The summands in equation <ref> are positive semi-definite such that for any k> N_Υ+1, truncating _k:1^β verifies

    ∑_l=1^N_Υ + 11/1 + β^k-l_k:l_l _k:l^≤_k:1^β≤ p_sup_n.

Note that by <ref>, if k>N_Υ+2
    _k:N_Υ + 1_N_Υ +1 :1_k:N_Υ + 1^= ∑_l=1^N_Υ +1_k:l_l _k:l^,

and therefore, for every k>N_Υ+2
    1/1+β^k -1 _k:N_Υ + 1_N_Υ +1 :1_k:N_Υ + 1^   ≤ p_sup_n.

Using <ref> (b), for every k>N_Υ+2 we derive

    1/1+β^k-N_Υ - 1_k:N_Υ +1_k:N_Υ +1^   ≤p_sup( 1+β)^N_Υ/b_n,

using the inequality in <ref>.  For any j, multiplying equation <ref> on the left by ^1_j_k^ and on the right by ^1_j_k and taking the limit as k→∞ shows that it is necessary for equation <ref> to hold for the left side to be bounded away from ∞.


In contrast to <ref> for autonomous systems, <ref> uses the <ref> (b) to simplify the arguments — this moreover guarantees the necessary criterion is with respect toλ_1, as the controllability matrix is guaranteed to be positive definite and thus nonvanishing on every Oseledec space.  There is, however, a more direct analogue to the statement of <ref> where the adjoint-covariant Lyapunov vectors will play the role of the eigenvectors of^.  It is easy to demonstrate that the adjoint-covariant Lyapunov vectors have the desired covariance and growth/decay with respect to the reverse time adjoint model,^_k.  There exist, under the condition of integrally separated Oseledec spaces, classical constructions for covariant and adjoint-covariant bases that decompose the model propagator into a block-upper-triangular form <cit.>[see Theorem 5.4.9]. This decomposition makes the derivation of a precise statement like <ref> analogous in time varying models, with respect to the adjoint-covariant Laypunov vectors and adjoint-covariant Oseledec spaces.  However, the above arguments require significant additional exposition which we feel unnecessary, as <ref> is sufficiently general.


Assume equations <ref> and <ref> satisfy <ref> (b) and suppose _k _k ≡ for every k such that α = β = 0.  Let k≥ 1 and choose ∈{^i_j_k : n_0 < i ≤ p,  1≤ j ≤κ_i} such that ł| |̊ =1.  There is a C>0 independent of k such that

    ^_k    ≤ C <∞.

The inequality in equation <ref> is an equality for the unfiltered forecast where β = α=0.  Thus the corollary is clear for any stable BLV directly from <ref> and the conclusion extends to norm one linear combinations. 
<ref> extends the intuition of AUS to the presence of model error: corrections may be targeted along the expanding modes while the uncertainty in the stable modes remains bounded by the system's dynamic stability alone.  Particularly, without filtering uncertainty remains uniformly bounded in the span of the stable BLVs. This is analogous to the results of Bocquet et. al. <cit.>, where in perfect models, the stable dynamics alone are sufficient to dissipate forecast error in the span of the stable BLVs.  Withα=0, the uniform bound in <ref> may be understood by the two components which equation <ref> is composed of, the bound on_k:k-N_i,ϵand 

    q_sup e^2(λ_i +ϵ)N_i,ϵ +1/1-e^2(λ_i+ϵ).

The controllability matrix_k:k-N_i,ϵrepresents the newly introduced uncertainty from model error that is yet to be dominated by the dynamics.  Equation <ref> represents an upper bound on the past model errors that have already been dissipated by the stable dynamics.  Nevertheless, this uniform bound is uninformative about the local variability.  In the following sections, we study the variance of the unfiltered uncertainty in the stable BLVs compared to the uncertainty of the Kalman filter.



§ NUMERICAL EXPERIMENTS


 §.§ Experimental setup


To satisfy <ref>, we construct a discrete, linear model from the nonlinear Lorenz-96 (L96) equations <cit.>, commonly used in data assimilation literature see, e.g., <cit.>[and references therein].
For eachm∈{1,⋯,n}, the L96 equations readd/d t ≜(),

    L^m()    =-x^m-2x^m-1 + x^m-1x^m+1 - x^m + F
  
such that the components of the vectorare given by the variablesx^mwith periodic boundary conditions,x^0=x^n,x^-1=x^n-1andx^n+1=x^1.  The termFin L96 is the forcing parameter.  The tangent-linear model <cit.>  is governed by the equations of the Jacobian matrix,∇(),

    ∇ L^m() = 0, ⋯, -x^m-1, x^m+1 - x^m-2, -1, x^m-1,0, ⋯, 0.
     

We fix the model dimensionn≜10and the forcing parameter as the standardF=8, as the model exhibits chaotic behavior, while the small model dimension makes the robust computation of Lyapunov vectors numerically efficient.  The linear propagator_kis generated by computing the discrete, tangent-linear model <cit.> from the resolvent of the Jacobian equation <ref> along a trajectory of the L96, with an interval of discretization atδ≜0.1.  We numerically integrate the Jacobian equation with a fourth order Runge-Kutta scheme with a fixed time step ofh ≜0.01.

ForF=8, the 10 dimensional nonlinear L96 model has a non-degenerate Lyapunov spectrum and we replace the superscripti_jwithifor the BLVs.  The model has three positive, one neutral and six negative Lyapunov exponents, such thatn_0 = 4.  The Lyapunov spectrum for the discrete, linear model is computed directly via the relationship in <ref>, where the average is taken over10^5iterations of the recursive QR algorithm, after pre-computing the BLVs to convergence.  In our simulations, before our analysis, we pre-compute the BLVs and the FLVs over10^5iterations of the recursive QR algorithm for the forward model, or respectively, for the reverse time adjoint model <cit.>[see section 3].  We note that the computed Lyapunov spectrum for the discrete, linear model as in simulations is related to the spectrum of the nonlinear L96 model by rescaling the linear model's exponents by1/δ.



 §.§ Variability of recurrent perturbations


While <ref> guarantees that the uncertainty in the stable BLVs is uniformly bounded, this bound strongly reflects the scale of the model error and the local variance of the Lyapunov exponents.  If model errors are large, or the stable Lyapunov exponents have high variance, this indicates that the uniform bound can be impractically large for forecasting.  
Assume no observational or filtering constraint, i.e.,_k_k≡.  Suppose that the model error statistics are uniform in time and spatially uncorrelated with respect to a basis of BLVs:_k ≜_k _k^, where∈ℝ^n×nis a fixed diagonal matrix with thei_j-th diagonal entry given byD_i_j.  Denote_0 ≡_0, then equation <ref> becomes

    ^ i_j_k^_k ^i_j_k     = ∑_l=0^k ^ i_j_k^_k:l_l ^_k:l^ i_j_k 
        =  ∑_l=0^k _i_j^_k:l^_k:l_i_j
        = D_i_j∑_l=0^k ł| ^_k:l^i_j|̊^2 ,

whereł| ^_k:l^i_j |̊is the norm of thei_j-th row of_k:l.  In equation <ref>,_krepresents the freely evolved uncertainty at timek, and thus∑_l=0^k ł| ^_k:l^i_j |̊^2describes the variance of the free evolution of perturbations in the direction of^i_j_k.

For each 1 ≤ i ≤ p, each 1 ≤ j ≤κ_i and any k, we define

    Ψ_k^i_j≜∑_l=0^k ł| ^_k:l^i_j|̊^2

to be the free evolution of perturbations in the direction of ^i_j_k.


Assuming the errors are uncorrelated in the basis of BLVs is a strict assumption, but studying the free evolution of perturbations has general applicability:_k^_k _k ≤q_sup_n, and therefore, equation <ref> may be interpreted in terms of an upper bound on the variance of the freely evolved forecast uncertainty in thei_j-th mode.  <ref> describes our recursive approximation of the free evolution, given by equation <ref>, fork∈{1,⋯,m}via the QR algorithm.  We assume that the QR algorithm has been run to numerical convergence for the BLVs at time0.


Equation <ref> implies ł| ^^i_j|̊^2 decays exponentially in k-l and the inner loop of <ref> needs only be computed to the first l such that ł| ^^i_j|̊^2 is numerically zero.


The approximation of <ref> with <ref> is numerically stable for allkand anyi > n_0, precisely due to the upper triangular dynamics in the BLVs.  The upper triangularity of all_kmeans the lower right block of_k:lis given as the product of the lower right blocks of the sequence of matrices{_j}_j= l+1^k.  Therefore, computing the stable block of_k:lis independent of the unstable exponents, and the row norms of_k:lconverge uniformly and exponentially to zero by <ref>.

In <ref> we plotΨ_k^5andΨ_k^6as in <ref> and the LLEs for^5_kand^6_kfork∈{1, ⋯, 10^4}.  Assuming that_k ≤q_sup_n,Ψ^i_kbounds the variance in thei-th stable mode at timek, up to the scaling factor ofq_sup.  Asn_0=4, the exponentλ_5is the stable exponent closest to zero.  The left side of <ref> corresponds to the exponentλ_5 ≈-0.0433while the right side corresponds to the exponentλ_6 ≈-0.0878.  The upper row in <ref> plots the evolution ofΨ^i_kfor^5_kand^6_k, while  the bottom row shows the corresponding time series of LLEs.  The mean of the LLEs are approximately equal to their corresponding Lyapunov exponent, while the standard deviation is given by0.142forλ_5and0.133forλ_6respectively.



WhileΨ_k^5is uniformly bounded, <ref> illustrates that it can be on the order of𝒪10^3, with a mean value of approximately 808 over the10^4iterations.  This is in contrast to perfect models where the projection of the unfiltered forecast error into any stable mode converges to zero at an exponential rate <cit.>.  Moreover, the frequency and scale of positive realizations of LLEs of^5_khas a strong impact the variance of the unfiltered error.  The fewer, and weaker, positive realizations of the LLEs of^6_kcorrespond to the lower overall uncertainty represented byΨ^6_k.  The maximum ofΨ_k^6is on the order ofØ(10^2), with a mean value of approximately 28. 




 §.§ Unfiltered versus filtered uncertainty


In the following, we compare the variance of the unfiltered error in the stable BLVs, represented byΨ_k^i, fori∈{5,⋯,10}, with the uncertainty in the Kalman filter.  Assuming that_k ≜_n, in this caseΨ_k^iis equal to the variance of the unfiltered error along^i_k.  While the error in the Kalman filter depends on the observational configuration, the value ofΨ_k^idepends only on the underlying dynamics.  Therefore, we benchmark the variance of the unfiltered error over a range of observational designs to determine under what conditions the unfiltered error in the stable BLVs will exceed the uncertainty of the full rank Kalman filter.  This analysis allows us to evaluate how many of the stable BLVs can remain unfiltered while achieving an acceptable forecast performance.  This comparison has a special significance when considering reduced rank, sub-optimal filters, which is the subject of the sequel <cit.>.

The recent works of Bocquet et al. <cit.> and Frank & Zhuk <cit.>, weaken <ref> to criteria on the observability, or detectability, of the unstable-neutral subspace to obtain filter stability and boundedness in perfect models.  The results in <ref>, <ref> and <ref> similarly suggest that the sufficient condition for filter boundedness, <ref>, may be weakened in the presence of model errors. For this reason, we will study the variance of the filtered error with respect to observations satisfying the criteria discussed by Bocquet et al. <cit.> and Frank & Zhuk <cit.>.

Given a fixed dimension of the observational spaced < n, consider finding an observational operator,_k, which minimizes the forecast uncertainty.  Suppose the singular value decomposition of an arbitrary choice of_kis given as

    _k = _k _k _k^.

For a given observation error covariance matrix, the size of the singular values of_kcorrespond to the precision of observations relative to the uncertainty in the precision matrix,_k ≜_k^^-1_k_k.
Imposing that all singular values of_kmust be equal to one, then up to an orthogonal transformation of^-1_k, we equate the choice of an observational operator_kwith the selection of an orthogonal matrix_k ∈ℝ^n ×d.

For perfect models,_k ≡0, we write the forecast error Riccati equation in terms of a choice of_k ≜_k^as

    _k+1   = _k+1_k
        _n + _k^-1/2^_k _k^_k^-1/2_k^_k
    ^-1^_k_k+1^.

The Frobenius norm,ł| _k+1 |̊_= √(_k+1^2) 
,is bounded by

    {_k+1_k
        _n + _k^-1/2^_k _k^_k^-1/2_k^_k
    ^-1^_k_k+1^}
    ≤   {
        _n + _k^-1/2^_k _k^_k^-1/2_k^_k
    ^-1}^_k_k+1^_k+1_k.

Equation <ref> attains its smallest values when the eigenvalues of

    _n + _k^-1/2^_k _k^_k^-1/2_k^_k

are as large as possible, similar to maximizing the denominator of equation <ref>.

For a fixed sequence of observation error covariances, finding the largest eigenvalues of equation <ref> can be studied by finding the subspace for which the matrix of orthogonal projection coefficients^_k _khas the largest singular values.  In perfect models, the forecast error covariance for the Kalman filter asymptotically has support confined to the span of the unstable and neutral BLVs <cit.>.  This is likewise, evidenced for the ensemble Kalman filter in weakly-nonlinear models <cit.>, suggesting that the columns of_kshould be taken as the leadingdBLVs. 

Given d ≥ 1, let ^1:d_k ∈ℝ^n× d denote the matrix comprised of the first d columns of _k.  We define the observation operator ^ bd_k ≜_k^1:d^.
<ref> is a formalization of the AUS observational paradigms <cit.> utilizing “bred vectors” as proxies for the BLVs.  The breeding method of Toth & Kalnay <cit.> simulates how the modes of fast growing error are maintained and propagated through the successive use of short range forecasts in weather prediction.  The bred vectors are formed by initializing small perturbations of a control trajectory and forecasting these in parallel along the control. Upon iteration, the span of these perturbations generically converge to the leading BLVs.  For a discussion of variants of this algorithm, and the convergence to the BLVs, see e.g., Balci et al.  <cit.>.

The choice of observation operator in <ref> is also related to the numerical study of targeted observations for the L96 model of Law et al. <cit.>.  Law et al. target observations with the eigenvectors of the operator^_k+1 _k+1, but note that for a small intervalδ≜t_k+1 - t_k, the difference between the linearized equations defining^_k+1 _k+1and_k+1^_k+1is negligible <cit.>[see Remark 5.1].  Law et al. suggest that the eigenvectors of either^_k+1 _k+1and_k+1^_k+1may be sensible depending on whether the filter should take into account the principle axes of growth from the past to the current time or from the present to future time.  It is clear from equations <ref> and <ref> that asδbecomes large, the eigenvectors of_k+1^_k+1approach the BLVs, whereas^_k+1 _k+1approach the FLVs.


Given d ≥ 1, let ^1:d_k ∈ℝ^n× d denote the matrix comprised of the first d columns of _k.  We define the observation operator ^ fd_k ≜_k^1:d^.


Note that the observation operator_k^b4uniformly-completely observes the span of the unstable and neutral BLVs, and thus ford≥4,_k^bdsatisfies the sufficient criterion for filter stability in perfect dynamics discussed by Bocquet et al. <cit.>.  The operator_k^b4likewise satisfies the necessary and sufficient detectability criterion for filter stability perfect dynamics of Frank & Zhuk <cit.>.  On the other hand, the operator^fd_kobserves the span of the leadingdFLVs.  Unlike the BLVs, the FLVs define a QL decomposition of the span of the covariant Lyapunov vectors <cit.>[see equation (53)].  This implies that the columns of the operator_k^f4actually spans the orthogonal complement to the stable Oseledec spaces.  Therefore,_k^f4satisfies the criterion of Frank & Zhuk <cit.>, but will not generally satisfy the condition of Bocquet et al. <cit.>.

We perform parallel experiments, fixing the sequence of linear propagators_k, and the initial prior error covariance_0 ≜_n, while varying the choice of the observation operator and the observational dimensiond.  In each parallel experiment, we study the average forecast uncertainty for the full rank Kalman filter as described by Frobenius norm of the forecast error covariance_k, averaged over10^5assimilations, neglecting a separate filter stabilization period of10^4assimilations.  For eachd∈{ 4, ⋯, 9}, we compare the following choices of observation operators: (i)_k^bd; (ii)_k^fd; (iii)_k ≜_k^for randomly drawn orthogonal matrices,_k ∈ℝ^n ×d; and (iv) a fixed network of observations, given by the leadingdrows of the identity matrix, i.e.,_k ≜_d×n.  We also compute the average Frobenius norm of the forecast error covariance for full dimensional observations, with_k ≜_n.  In each experiment, the observational and model error covariances are fixed as_k ≜_dand_k ≜_n.  For eachi, the value ofΨ^i_kis averaged over the10^5assimilations.

In <ref>, we plot the average Frobenius norm of the Kalman filter forecast error covariance matrix as a function of the number of observations,d.  We consider the observation configurations^bd_k,^_kand_n(plotted horizontally).  The average values ofΨ_k^ifori=7, ⋯,10are also plotted horizontally.  While the observational dimensiond < 7, the average uncertainty for the Kalman filter with random observations, or observations in the BLVs, is greater than the average variance of the unfiltered error along^7_k.  Similarly, in <ref> we consider the configurations with observations defined by_k^bd,^fd_kand_d ×n.  The average values ofΨ_k^ifori=5, ⋯,8are plotted horizontally.  The variance of the unfiltered error in^5_kexceeds the uncertainty of the Kalman filter in every configuration.  The Kalman filter with observations fixed, or in the FLVs, do not obtain comparable performance with the unfiltered error in^6_kuntild≥6.  Only the LLEs of^i_kfori=8,9,10are sufficiently stable to bound the unfiltered errors below the Kalman filter with a fully observed system.



Our results have strong implications for the necessary rank of the gain in ensemble-based Kalman filters.  In perfect, weakly nonlinear models, the ensemble span typically aligns with the leading BLVs <cit.>.  From the above results, we conclude that the effective rank of the ensemble-based gain must be increased to account for weakly stable BLVs of high variance in the presence of model errors.  The perturbations of model errors excited by transient instabilities in these modes can lead to the unfiltered errors becoming unacceptably large compared to the filtered errors.  

In <ref> and <ref>, the choice of observations in the span of the leading BLVs dramatically outperforms the observations in the span of the leading FLVs, or fixed observations.  Likewise^bd_kmakes a slight reduction to the overall forecast error over a choice ofdrandom observations.  As the span of the leadingn_0FLVs is orthogonal to the trailing, stable Oseledec spaces, this choice is can be considered closer to the minimum necessary observational constraint on the forecast errors.  Particularly, the kernel of^f n_0_kis identically equal to the sum of the stable Oseledec spaces.  This suggests that a necessary and sufficient condition for filtered boundedness can be described in terms of the observability of then_0leading FLVs, similar to the criterion of Frank & Zhuk <cit.>.  While it is not necessary, the sufficient condition of Bocquet et al. <cit.> leads to a lower filter uncertainty as the span of the leadingn_0BLVs generally contains the largest projection of the forecast error.  This suggest that observing the leading eigenvectors of_k+1^_k+1may generally outperform observing the leading eigenvectors of^_k+1 _k+1when the time between observationsδ= t_k+1 -t_kleads to significant differences in the linear expansions, as was noted as an alternative design by Law et al. <cit.>.  For operational forecasting, this supports the use of the breeding technique <cit.> to target observations, over using the axes of forward growth. 




§ CONCLUSION


This work formalizes the relationship between the Kalman filter uncertainty and the underlying model dynamics, so far understood in perfect models, now in the presence of model error.  Generically, model error prevents the collapse of the covariance to the unstable-neutral subspace and our <ref> and <ref> characterize the asymptotic window of uncertainty.  We provide a necessary condition for the boundedness of the Kalman filter forecast errors for autonomous and time varying dynamics in <ref> and <ref>: the observational precision, relative to the background uncertainty, must be greater than the leading instability which forces the model error.  Particularly, <ref> proves that forecast errors in the span of the stable BLVs remain uniformly bounded, in the absence of filtering, by the effect of dynamic dissipation alone.  

The uniform bound on the errors in the span of the stable BLVs extends the intuition of AUS to the presence of model error, but with qualifications.  Studying this uniform bound with <ref>, we identify an important mechanism for the growth of forecast uncertainty in sub-optimal filters: variability in the LLEs for asymptotically stable modes can produce transient instabilities, amplifying perturbations of model error.  The impact of stable modes close to zero differs from similar results for nonlinear, perfect models by Ng. et. al. <cit.>, and Bocquet et. al. <cit.>, where the authors demonstrate the need to correct weakly stable modes in the ensemble Kalman filter due to the sampling error induced by nonlinearity.  Likewise, this differs from the EKF-AUS-NL of Palatella & Trevisan <cit.>, that accounts for the truncation errors in the estimate of the forecast uncertainty in perfect, nonlinear models.  Our work instead establishes the fundamental impact of these transient instabilities as a linear effect in the presence of model errors.

In addition to our necessary criterion for filter boundedness, in <ref> we discuss the criteria of Bocquet et al. <cit.> and Frank & Zhuk <cit.> in relation to dynamically targeted observations.  Our numerical results suggest how these sufficient, and respectively necessary and sufficient, criteria can be extended to the presence of model errors.   Moreover, we distinguish between the minimal necessary observational constraints for filter boundedness and more operationally effective, sufficient designs.  Particularly, our results suggest that while it may be necessary that the observations uniformly completely observe the span of the unstable-neutral FLVs, it is sufficient and improves performance to uniformly completely observe the span of unstable-neutral BLVs.  In terms of operational forecasting, this strongly supports the use of bred vectors to target observations to constrain the forecast errors.

  <ref>, <ref>, <ref> and the results of <ref> suggest that as a theoretical framework for the ensemble Kalman filter, AUS may be extended to the presence of model errors.  By uniformly completely observing and correcting for the growth of uncertainty in the span of the unstable, neutral and some number of stable BLVs, reduced rank filters in the presence of model errors may obtain satisfactory performance. In practice, one may compute off-line the typical uncertainty in the stable BLVs via <ref> and determine the necessary observational and ensemble dimension at which unfiltered forecast error has negligible impact on predictions.  However, computational limits on ensemble sizes may make this strategy unattainable in practice — the impact of these unfiltered errors on the performance of a reduced rank, sub-optimal filter is the subject of the direct sequel to this work <cit.>.



§ ACKNOWLEDGMENTS

The authors thank two anonymous referees, and their colleagues Karthik Gurumoorthy, Amit Apte, Erik Van Vleck, Sergiy Zhuk and Nancy Nichols, for their valuable feedback and discussions on this work.  CEREA is a member of the Institut Pierre-Simon Laplace (IPSL).
 
plain