observational noises. Due to very large n, it is impossible to apply traditional estimation algorithms for producing the estimate x^k and that is the reason there exist different approximation algorithms for solving this estimation problem. Theoretically, the optimal in mean square error (MSE) estimate x^k=N based on the set of observations Z½1, N� :¼ {z1,…; zN} is a filtered estimate for N = k and smoothed estimate for N > k [1, 2]. For the linear dynamical system Eqs. (1.1) and (1.2), the computation of x^k :¼ x^k=k can be efficiently performed using the Kalman filter (KF) [3] which is a sequential procedure. The KF provides also the equations for computation of estimation errors. If we are interested in obtaining x^k=N—the best estimate for xk based on Z[1, N], the Kalman smoother can serve as an efficient algorithm for its computation. The KF approach, however, is inappropriate for solving estimation problems in highdimensional systems. In this chapter, the high-dimensional system means the system whose

capacity, in both computational power and memory, is still very far to be sufficient to implement the KF in real time to produce the filtered estimate and to make corresponding forecast. For suboptimal schemes for atmospheric data assimilation based on the KF, see Ref. [4].

In this chapter, the emphasis is put principally on comparison of the AF with VM. For the review on the data assimilation methods in meteorology and oceanography, see Ref. [5]. To see more the advantages of the AF, we implement the extended KF (EKF) [6] in Section 6 and will compare its performance with that of the AF (the experiment with Lorenz system). The Cooper-Haines filter (CHF) [7], widely used in data assimilation in oceanography, is also applied in Section 7 to produce the estimate for the ocean state. It serves as a reference to be compared

In the next section, the variational method (VM), which is widely used in data assimilation for high-dimensional systems in meteorology and oceanography, is outlined. Section 3 provides the recently developed AF approach to data assimilation. The main idea of the AF is to take the innovation representation for the input-output system as a departure point to formulate the optimization problem, with the parameters of the filter gain as control variables. Section 4 presents the tools to implement the AF in a simple and efficient way which is adapted for highdimensional setting. This includes the objective function, filter stability, structure of the error covariance matrix (ECM), gain parameterization, algorithm for optimization known as simultaneous perturbation stochastic approximation (SPSA). It is shown how the ECM is estimated

E½vk� ¼ 0, E½vkvT

E½ðx0−x0ÞwT

0 � ¼ M0,

l � ¼ δklQ,

k � ¼ 0:

, and there are uncertainties in statistics of the model and

. At the present and in the near future, the computer

(1.3)

l � ¼ δklR, E½wkvl� ¼ 0,

teristics of the entering random variables are given as

302 Advances in Statistical Methodologies and their Application to Real Problems Advances in Statistical Methodologies and Their Application to Real Problems

dimension n of xk is of order 106 – 108

state dimension is of order 106 – 108

with that produced by the AF in high-dimensional setting.

Consider the problem of estimating {xk} in Eqs. (1.1)–(1.3). The VM consists of minimizing the following objective function

$$f[\mathbf{x}\_0, \dots, \mathbf{x}\_N] = \mathfrak{e}\_0 M\_0^{-1} \mathfrak{e}\_0 + $$

$$\sum\_{k=1}^N (\mathbf{z}\_k - H\_k \mathbf{x}\_k)^T R^{-1} (\mathbf{z}\_k - H\_k \mathbf{x}\_k), \mathfrak{e}\_0 := \mathbf{x}\_0 - \overline{\mathbf{x}}\_0,\tag{2.1}$$

$$J[\mathbf{x}\_0, \dots, \mathbf{x}\_N] \to \min\_{\left[\mathbf{x}\_0, \dots, \mathbf{x}\_N\right]} \tag{2.2}$$

under the constraints ð1:1Þ (2.3)

Thus, in the VM, we seek optimal solutions in the functional space (space of functions {xk}). For systems of high dimension, this task is impossible to perform. The simplification is required. Suppose the system Eq. (1.1) is linear and perfect, that is,

$$\mathbf{x}\_{k+1} = \Phi\_k \mathbf{x}\_k, k = 0, 1, \dots \tag{2.4}$$

Expressing all xk as functions of the initial state x0,

$$\mathbf{x}\_{k} = \Phi(k, \mathbf{0})\mathbf{x}\_{0},$$

$$\Phi(k, l) = \Phi\_{k-1}\dots\Phi\_{l}, (k > l), \Phi(k, k) = I,\tag{2.5}$$

$$I \text{ is the identity matrix of appropriate dimension}$$

and substituting xk, ∀k Eq. (2.5) into Eq. (1.1), at each k th observation instant, the following set of observations is available for x0,

$$\boldsymbol{z}\_k^1 = H\_k^1 \mathbf{x}\_0 + \boldsymbol{v}\_1^k, k = 1, 2, \dots$$

$$H\_k^1 := [(H\_1 \boldsymbol{\Phi}(\mathbf{1}, \mathbf{0}))^T, \dots, (H\_k \boldsymbol{\Phi}(k, \mathbf{0}))^T]^T,\tag{2.6}$$

$$\boldsymbol{v}\_k^1 = [\boldsymbol{v}\_1^T, \dots, \boldsymbol{v}\_k^T]^T.$$

Under the assumption on perfect model, the optimization problem Eqs. (2.1)–(2.3) is simplified as

$$J[\mathbf{x}\_0] \to \min\_{\left[\mathbf{x}\_0\right]},\tag{2.7}$$

$$J[\mathbf{x}\_0] := \mathbf{e}\_0^T \mathbf{M}\_0^{-1} \mathbf{e}\_0 + \sum\_{k=1}^N (\mathbf{z}\_k \mathbf{-} \boldsymbol{H}\_k^\prime \mathbf{x}\_0)^T \mathbf{R}\_k^{-1} (\mathbf{z}\_k \mathbf{-} \boldsymbol{H}\_k^\prime \mathbf{x}\_0), \tag{2.8}$$

$$\boldsymbol{H}\_k^\prime := \boldsymbol{H}\_k \boldsymbol{\Phi}(\mathbf{k}, \mathbf{0}).$$

We have now the unconstrained optimization problem Eqs. (2.7) and (2.8) with the vector of unknown parameters θ :¼ x0—the initial state. This problem can be solved using standard optimization techniques [8].

It is not hard to write out a solution to the problem Eqs. (2.7) and (2.8). For high-dimensional systems, there is no computational and memory resources to handle such implementation. In practice, the solution to the problem Eqs. (2.7) and (2.8) is found by solving iteratively the equation

$$\begin{array}{c} \nabla\_{\theta} \mathsf{J}[\theta] = \mathbf{0}, \\ \nabla\_{\theta} \mathsf{J}[\theta] := \left[ \eth \mathsf{J} / \partial \theta\_1, \dots, \eth \mathsf{J} / \partial \theta\_1 \right]^T. \end{array} \tag{2.9}$$

Comment 2.1. Usually, finding a solution to Eqs. (2.7) and (2.8) is a heavy task: in addition to storing the model solution produced by the direct model, minimization requires 20–30 iterations to reach a relatively good approximate solution.

Comment 2.2. Writing out ∇θJ½θ� shows that solving Eq. (2.9) requires

$$\Phi \left( H\_k' \right)^T y = \Phi\_k^T \Phi\_{k-1}^T \dots \Phi\_1^T H\_k^T y \tag{2.10}$$

for some y. As ΦT k is impossible to store, the approach known as adjoint equation (AE) is used which requires to construct AE code for computing the product ΦT k y. Each iteration in minimization of Eqs. (2.7) and (2.8) thus requires one integration of the model over the assimilation period, followed by one adjoint integration. The cost of one adjoint integration is about twice the cost of one direct integration, so that one minimization requires the equivalent of between 50–100 integrations of the model over the assimilation period (p. 205 [9]). In the next section, we see that the SPSA can also be used to solve this problem at a much lower cost.

Comment 2.3 As θ—initial state—has a physical meaning, it is important to introduce constraints on the appropriate physically realistic structure of the correction for θ during the estimation process. A poor (in the physical sense) structure of the guess for the initial state can lead to large estimation errors.

## 3. Adaptive filtering (AF)

z1 k ¼ H1

H1

32 Advances in Statistical Methodologies and Their Application to Real Problems

J½x0� :¼ eT

tions to reach a relatively good approximate solution.

Comment 2.2. Writing out ∇θJ½θ� shows that solving Eq. (2.9) requires

ðH′ kÞ T y ¼ ΦT k ΦT k−1…ΦT 1 HT

which requires to construct AE code for computing the product ΦT

optimization techniques [8].

0 M−1

0 e0 þ ∑ N k¼1

as

equation

for some y. As ΦT

can lead to large estimation errors.

kx0 þ vk

v1 k ¼ ½vT

k :¼ ½ðH1Φð1; 0ÞÞT,…; ðHkΦðk, 0ÞÞT�

Under the assumption on perfect model, the optimization problem Eqs. (2.1)–(2.3) is simplified

ðzk−H′

We have now the unconstrained optimization problem Eqs. (2.7) and (2.8) with the vector of unknown parameters θ :¼ x0—the initial state. This problem can be solved using standard

It is not hard to write out a solution to the problem Eqs. (2.7) and (2.8). For high-dimensional systems, there is no computational and memory resources to handle such implementation. In practice, the solution to the problem Eqs. (2.7) and (2.8) is found by solving iteratively the

∇θJ½θ� :¼ ½∂J=∂θ1,…; ∂J=∂θ1�

Comment 2.1. Usually, finding a solution to Eqs. (2.7) and (2.8) is a heavy task: in addition to storing the model solution produced by the direct model, minimization requires 20–30 itera-

zation of Eqs. (2.7) and (2.8) thus requires one integration of the model over the assimilation period, followed by one adjoint integration. The cost of one adjoint integration is about twice the cost of one direct integration, so that one minimization requires the equivalent of between 50–100 integrations of the model over the assimilation period (p. 205 [9]). In the next section,

Comment 2.3 As θ—initial state—has a physical meaning, it is important to introduce constraints on the appropriate physically realistic structure of the correction for θ during the estimation process. A poor (in the physical sense) structure of the guess for the initial state

we see that the SPSA can also be used to solve this problem at a much lower cost.

kx0Þ TR−1 k ðzk−H′

H′

∇θJ½θ� ¼ 0,

k is impossible to store, the approach known as adjoint equation (AE) is used

1, k ¼ 1; 2;…

1 ,…; vT k � T:

T,

J½x0� ! min ½x0�, (2.7)

kx0Þ,

T: (2.9)

k y (2.10)

k y. Each iteration in minimi-

k :¼ HkΦðk, 0Þ:

(2.6)

(2.8)

To overcome the difficulties listed in Comments 2.1–2.3, an adaptive filtering (AF) has been proposed in [10]. The main difference of the AF with the VM is lying in the choice of innovation representation for the original input-output system Eqs. (1.1) and (1.2) as a departure point to formulate the optimization problems. It is well known that under standard conditions, the optimal in MSE estimate x^k can be obtained by the KF. As the innovation process for the system output in the KF forms a white sequence, Kailath [11] has developed an innovation approach, in an elegant way, to derive the optimal filter for more general linear systems like nonstationary, filtering problems with Markovian processes for the model and observation errors. The innovation approach to linear least-squares approximation problems is first to "whiten" the output data and then to treat the resulting simpler white-noise observation problem. Consider the observation (output) sequence zk. The innovation process, associated with zk, is written as

$$\zeta\_k = z\_k - \mathbb{E}[z\_k | z\_{k-1}^1],\tag{3.1}$$

Under standard conditions (Gaussianness, uncorrelated noise sequences …), E½zkjz1 k−1� ¼ Hkx^k=k−1 hence

$$
\zeta\_k = z\_k - H\_k \hat{\mathbf{x}}\_{k/k-1}, \hat{\mathbf{x}}\_{k/k-1} = \Phi\_{k-1} \hat{\mathbf{x}}\_{k-1}, \tag{3.2}
$$

where x^k=k−1 is an optimal in MSE one-step ahead prediction for xk given z1 k−1. Using ζk instead of zk, one can write out the formula for the estimate x^k and the KF under standard conditions. The filter has the form

$$\ddot{\mathbf{x}}\_k = \Phi\_k \ddot{\mathbf{x}}\_{k-1} + \mathbf{K}\_k \zeta\_k,$$

$$\mathbf{K}\_k = \mathbf{M}\_k \mathbf{H}\_k^T \left[ \mathbf{H}\_k \mathbf{M}\_k \mathbf{H}\_k^T + \mathbf{R}\_k \right]^{-1} \tag{3.3}$$

where Mk is the ECM for the prediction x^k=k−1. This matrix is found as a solution to the Riccati equation

$$M\_k = \Phi\_k P\_k \Phi\_k^T + Q\_k,\\ P\_k = [I - K\_k H\_k] M\_k. \tag{3.4}$$

Due to the very expensive computational burden in time stepping the ECM Mk in Eq. (3.4) as well as insufficient memory storage, the KF is impractical for solving data assimilation problems in very high-dimensional setting. The idea of the AF is based on the fact that when the filter is optimal, the innovation ζk has a minimum variance. If we assume that the gain Kk belongs to a set of parameterized gains, that is,

$$K\_k = K\_k(\theta), \theta \in \Theta,\tag{3.5}$$

the optimal AF can be considered as that in some class of parameterized filters of a given structure. The following objective function is introduced

$$J(\theta) = E[\Psi(\theta)] \to \min\_{\theta \in \Theta} \Psi(\zeta\_k) = ||\zeta\_k||\_{\Sigma\_k^{-1}}^2 ||\zeta\_k||\_{\Sigma\_k^{-1}}^2 := <\zeta\_k, \Sigma\_k^{-1}\zeta\_k > . \tag{3.6}$$

In Ref. [12], the different classes of parameterized filters are found which belong to the class of stable reduced-order filters (ROF) [10, 13].

As an example for one class of ROFs, consider

$$K\_k = P\_{r,k} K\_{e,k} \tag{3.7}$$

where Ke, k : Rp ! Rne represents the gain, mapping the innovation vector from the observational space to the reduced space Rne of dimension ne≤n; Pr, k is mapping the reduced space Rne to the full space Rn. The choice of a reduced space is of primary importance since it depends on the main characteristics of the filter known as stability. As proved in [12], under detectability condition, stability of the filter is ensured by forming the columns of Pr, k from unstable and stable eigenvectors (or singular vectors, Schur vectors) of the fundamental matrix Φk, and one can choose

$$K\_{e,k} = H\_{e,k}^T \left[ H\_{e,k} H\_{e,k}^T(k) + R\_k \right]^{-1}, \\ H\_{e,k} := H\_k P\_{r,k}, \tag{3.8}$$

One class of parameterized filters is (Section 5.2.2 in Ref. [12])

$$\begin{array}{c} \mathcal{K}\_{k}(\boldsymbol{\theta}) = P\_{r,k} \Lambda \mathcal{K}\_{e,k}(\boldsymbol{\theta}),\\ \Lambda = \text{diag}\left[\boldsymbol{\theta}\_{1}, \dots, \boldsymbol{\theta}\_{n\_{\boldsymbol{\epsilon}}}\right],\\ \mathbf{1} - \mathbf{1}/|\boldsymbol{\phi}\_{i}| < \dot{\boldsymbol{\phi}}\_{1}(i) \boldsymbol{\Leftarrow} \boldsymbol{\theta}\_{2}(i) < \mathbf{1} + \mathbf{1}/|\boldsymbol{\phi}\_{i}|.\end{array} \tag{3.9}$$

if φi is an unstable or neutral eigenvector of Φ. For the stable φi , we have

$$0 < \dot{\rho}\_1(\dot{\mathbf{i}}) \le \theta\_i \le \dot{\theta}\_2(\dot{\mathbf{i}}) < \mathcal{2},\tag{3.10}$$

#### 4. Differences between VM and AF

#### 4.1. Batch data formulation

We list now the main differences between two approaches VM and AF from which it becomes clear what are the advantages of the AF over the VM.

To make easier comparison between two approaches, let us write out the objective function Eq. (3.6) using a representation in a sample space

$$f\_N(\theta) \simeq f\_N(\theta) = \frac{1}{N} \Sigma\_{k=1}^N (z\_k - H\Phi\_k \hat{x}\_{k-1}(\theta))^T \Sigma\_k^{-1} (z\_k - H\Phi\_k \hat{x}\_{k-1}(\theta)) \tag{4.1}$$

Mention that in practical implementation of the AF, the optimization algorithm is not constructed on the basis of Eq. (4.1), but on Eq. (3.6). That is, due to the fact that Eq. (4.1) is written in a batch form which requires to make optimization over the time interval ½1, N� resulting in a very high computational burden. Minimizing Eq. (3.6) allows to apply SPSA method which is much less consuming for both computational and memory requirements. Below the main differences between VM and AF are listed:

JðθÞ ¼ E½ΨðθÞ� ! min θ∈Θ, ΨðζkÞ ¼ jjζkjj2

34 Advances in Statistical Methodologies and Their Application to Real Problems

stable reduced-order filters (ROF) [10, 13].

can choose

As an example for one class of ROFs, consider

Ke, k ¼ HT

One class of parameterized filters is (Section 5.2.2 in Ref. [12])

1−1=jφi

if φi is an unstable or neutral eigenvector of Φ. For the stable φi

4. Differences between VM and AF

clear what are the advantages of the AF over the VM.

Eq. (3.6) using a representation in a sample space

JðθÞ ≈ JNðθÞ ¼ 1

N ∑N

4.1. Batch data formulation

e, k½He, kHT

Σ−1 k , jjζkjj2 Σ−1 k

In Ref. [12], the different classes of parameterized filters are found which belong to the class of

where Ke, k : Rp ! Rne represents the gain, mapping the innovation vector from the observational space to the reduced space Rne of dimension ne≤n; Pr, k is mapping the reduced space Rne to the full space Rn. The choice of a reduced space is of primary importance since it depends on the main characteristics of the filter known as stability. As proved in [12], under detectability condition, stability of the filter is ensured by forming the columns of Pr, k from unstable and stable eigenvectors (or singular vectors, Schur vectors) of the fundamental matrix Φk, and one

e, kðkÞ þ Rk�

−1

KkðθÞ ¼ Pr, kΛKe, kðθÞ, Λ ¼ diag ½θ1,…; θne �,

j:

0 < ò1ðiÞ ≤ θi ≤ ò2ðiÞ < 2, (3.10)

JNðθÞ ! min θ∈Θ,

k ðzk−HΦkx^k−1ðθÞÞ (4.1)

, we have

j < ò1ðiÞ≤θi≤ò2ðiÞ < 1 þ 1=jφi

We list now the main differences between two approaches VM and AF from which it becomes

To make easier comparison between two approaches, let us write out the objective function

k¼1ðzk−HΦkx^k−1ðθÞÞTΣ−1

Mention that in practical implementation of the AF, the optimization algorithm is not constructed on the basis of Eq. (4.1), but on Eq. (3.6). That is, due to the fact that Eq. (4.1) is written in a batch form which requires to make optimization over the time interval ½1, N�

:¼< ζk, Σ−1

Kk ¼ Pr, kKe, k (3.7)

, He, k :¼ HkPr, k, (3.8)

(3.9)

k ζk > : (3.6)

(D1) Dynamical system (DS): if in Eq. (2.1), the DS is the initial system Eq. (1.1); in Eq. (4.1), the DS is the filtering Eq. (3.3). This difference has an interesting consequence: if in practice, there is very little known about statistics of wk, the sequence ζk is observed, and hence, it is possible to estimate the statistics of ζk.

(D2) The system noise wk in Eq. (1.1) is white, while in Eq. (3.3), ζk is a white sequence only if the filter is optimal. That allows us to easily apply different statistical tests for verifying the optimality of the assimilation procedure.

(D3) Control variable x0 in the VM is the initial state, whereas the control variable in Eq. (3.6) is the parameter vector θ.

This difference has an important consequence: as x0 has to be of precise physical meaning (depending, for example, on the ocean domain of interest), the structure for the guess θ0 :¼ x^0 0 (for the initial state) as well as correction δx^ν 0, generated by iterative algorithm, must be chosen carefully so that at each ν iteration, the estimate x^ν 0, x^ν 0 ¼ x^ν−1 0 þ δx^ν 0, must be of physically realistic structure. This is not an easy task. On the other hand, in the AF, the parameters usually are immaterial [see θ Eq. in (3.10)]; hence, the choice of structure for θ is of no importance.

(D4) Suppose the DS Eq. (1.1) is unstable. It implies that the error in estimating x0 will grow during integration of the direct and AE. As for the AF, by its construction, the filtering system Eq. (3.3) remains stable. This can be seen by representing the filtering Eq. (3.3) through its fundamental matrix Lk,

$$\begin{aligned} \hat{\mathfrak{x}}\_{k} &= L\_{k}\hat{\mathfrak{x}}\_{k-1} + K\_{k}\upsilon\_{k} \\ L\_{k} &= (I - \mathcal{K}\_{k}H\_{k})\mathfrak{O}\_{k}. \end{aligned} \tag{4.2}$$

As shown in Ref. [12], the filter Eq. (4.2) is stable under the conditions Eqs. (3.9) and (3.10). It means that the filtering error is bounded during model integration since the parameters θi are lying in the interval guaranteeing a stability of the filter Eq. (4.2).

(D5) Return to the objective function Eqs. (2.1)–(2.3). First taking the derivative of the objective function Eqs. (2.1)–(2.3) wrt (with respective to) x0, we have

$$\begin{split} \frac{1}{2} \nabla\_{\mathbf{x},0} \boldsymbol{I}[\hat{\mathbf{x}}\_{0}^{\boldsymbol{\nu}}] &= M\_{0}^{-1} \boldsymbol{e}\_{0}^{\boldsymbol{\nu}} - \sum\_{k=1}^{N} \boldsymbol{\Phi}^{T}(\boldsymbol{k}, \mathbf{0}) H\_{k}^{T} R\_{k}^{-1} (\mathbf{z}\_{k} - H\_{k} \boldsymbol{\Phi}(\boldsymbol{k}, \mathbf{0}) \mathbf{x}\_{0}^{\boldsymbol{\nu}}) \\ &= M\_{0}^{-1} \boldsymbol{e}\_{0}^{\boldsymbol{\nu}} - \sum\_{k=1}^{N} \boldsymbol{\Phi}^{T}(\boldsymbol{k}, \mathbf{0}) H\_{k}^{T} R\_{k}^{-1} (H\_{k} \boldsymbol{\Phi}(\boldsymbol{k}, \mathbf{0}) \mathbf{e}\_{0}^{\boldsymbol{\nu}} + \boldsymbol{\nu}\_{k}), \mathbf{e}\_{0}^{\boldsymbol{\nu}} := \mathbf{x}\_{0} - \hat{\mathbf{x}}\_{0}^{\boldsymbol{\nu}}. \end{split} \tag{4.3}$$

One sees that for a batch of N observations, Eq. (4.3) requires computation of N terms (without counting for the term M−1 0 eν 0). The k th term is associated with the assimilation instant k, and one needs to compute first μk :¼ Φðk, 0Þeν 0, that is, to integrate k times the direct model Φκ for κ ¼ 1,…; k and next to integrate backward (k times also) the AE ΦT κ from κ ¼ k to κ ¼ 1, that is, to compute ΦTðk, 0ÞHT κR−1 κ ðHkμk þ vkÞ. The larger the k, the bigger the amplification of the initial error eν 0 and the observation error vk. The error eν 0 is amplified doubly since it is integrated by the direct and adjoint models. But the amplification of vk (and wk when wk=¼ 0) is most worrying since it is integrated in the gradient estimate, making the gradient direction to be, possibly, completely erroneous.

#### 4.2. Implementation of AF

#### 4.2.1. The choice of criteria Eq. (3.6)

The choice of Eq. (3.6) is important in many aspects in order to obtain a simple and efficient data assimilation algorithm. The idea lying in the criteria Eq. (3.6) is to select some pertinent parameters as a control variables for minimizing the mean of the cost function ΨðζkÞ. For example, for the class of filters Eqs. (3.3), (3.7)–(3.10) the vector θ ¼ ðθ1,…; θne Þ T can be chosen as control vector for the problem Eq. (3.6).

The solution to the problem Eq. (3.6) can be found iteratively using a stochastic optimization (SA) algorithm

$$
\Theta\_{k+1} = \Theta\_k \mathbf{-} \mathbf{a}\_k \nabla\_{\theta} \Psi(\zeta\_{k+1}) \tag{4.4}
$$

where {ak} is a sequence of positive scalars satisfying some conditions to guarantee a convergence of the estimation procedure. The standard conditions are

$$a\_k \to 0, \sum\_{k=1}^{\infty} a\_k = \Leftrightarrow, \sum\_{k=1}^{\infty} a\_k^2 < \Leftrightarrow \tag{4.5}$$

The algorithm Eq. (4.4) is much more simple [compared to the computation of Eq. (4.3)] since it requires, at the kth assimilation instant, to compute only the gradient of the sample cost function ΨðζkÞ. The gradient ∇θΨðθkÞ of the sample objective function ΨðθkÞ can be computed using the AE approach (in what follows, for simplicity, the subscript k will be omitted to shorten the notations).

$$\frac{1}{2}[\delta\Psi(\zeta\_{k+1})]\_{\theta\_k} = -(H\Phi P\_r \delta\Lambda K\_t \zeta\_k, \zeta\_{k+1}) = -(\delta\Lambda K\_t \zeta\_k, \zeta\_{k+1}'), \zeta\_{k+1}' := P\_r^T \Phi^T H^T \zeta\_{k+1} \tag{4.6}$$

Thus, minimization of Eq. (3.6) by gradient-based SA algorithm requires only one integration of the direct model and one backward integration of the AE code: direct integration of x^k for producing the forecast x^kþ1=k ¼ Φx^k and backward integration ΦTHkζk−1 in computation of ζ′ kþ1. For the structure of the gain Eq. (3.9), the objective function Ψ is quadratic wrt θ; hence, one can find easily the optimal parameters.

A less computational burden can be achieved by measuring the sample objective function (but not based on a gradient formula): instead of computing the gradient by Eq. (4.6) based on AE, one can approximate the gradient using the values of the cost function [on the basis of finite difference scheme (FDSA)]. Traditionally, the ith component of the gradient can be approximated by

A Comparison Study on Performance of an Adaptive Filter with Other Estimation Methods for State... http://dx.doi.org/10.5772/67005 37

$$\nabla\_{\theta\_i} \Psi(\theta\_k) = \mathbf{g}\_i = \left[ \Psi(\theta\_k + \mathfrak{c}\_k \mathbf{e}\_i) \neg \Psi(\theta\_k - \mathfrak{c}\_k \mathbf{e}\_i) \right] / (2\mathfrak{c}\_k) \tag{4.7}$$

where ei is the unit vector with 1 in the ith component, 0 otherwise.

It is seen that FDSA algorithms do not require the formula for the gradient. However, for the high-dimensional systems (n≈Oð106 Þ−Oð107 Þ), this algorithm is inapplicable due to component-wise derivative approximation: for approximation of each partial derivative of the cost function, we need to make two integrations of the direct model.

In order to overcome the difficulties with very high dimension of θ, recently the class of algorithms known as simultaneous perturbation SA (SPSA) receives a great interest [14, 15]. The algorithm SPSA is of the same structure as that of FDSA Eq. (4.7), with the difference residing in the way to perturb stochastically and simultaneously all the components of θ. Concretely, let Δk ¼ ðΔk, 1,…; Δk,nÞ T be a random vector, Δk,i, i ¼ 1,…; n are Bernoulli independent identically distributed (iid). The gradient of the objective function is estimated as

$$\begin{aligned} \nabla\_{\theta} \Psi(\theta\_k) &= \mathbf{g} = (\mathbf{g}\_1, \dots, \mathbf{g}\_n)^T, \\ \mathbf{g} &= [\Psi(\theta\_k + \mathbf{c}\_k \Delta\_k) \mathbf{-} \Psi(\theta\_k \mathbf{-c}\_k \Delta\_k)] \Delta\_k^{-1} / (2\mathbf{c}\_k), \\ \Delta\_k &= (\Delta\_{k,1}, \dots, \Delta\_{k,n})^T, \Delta\_k^{-1} := (1/\Delta\_{k,1}, \dots, 1/\Delta\_{k,n})^T. \end{aligned} \tag{4.8}$$

It is seen that in the SPSA, all the directions are perturbed at the same time (the numerator is identical in all n components). Thus, SPSA uses only two (or three) times integrations of the model, independently on the dimension of θ which makes it possible to apply to high-dimensional optimization problems. Generally, SPSA converges in the same number of iterations as FDSA, and it follows approximately the steepest descent direction, behaving like the gradient method [14]. On the other hand, SPSA, with the random search direction, does not follow exactly the gradient path. On average, though, it tracks the gradient nearly because the gradient approximation is an almost unbiased estimator of the gradient, as shown in Ref. [15].

For the SPSA algorithm, the conditions for {ak} and {ck} are

$$\begin{aligned} \boldsymbol{a}\_{k} > 0, \boldsymbol{c}\_{k} > 0, \quad \boldsymbol{a}\_{k} \to 0, \quad \boldsymbol{c}\_{k} \to 0, \\ \sum\_{k=1}^{\infty} \boldsymbol{a}\_{k} = \rightsquigarrow, \sum\_{k=1}^{\infty} \left( \boldsymbol{a}\_{k} / \mathbf{c}\_{k} \right)^{2} < \rightsquigarrow \end{aligned} \tag{4.9}$$

#### 4.2.2. On the operator Pr

is, to compute ΦTðk, 0ÞHT

4.2. Implementation of AF

be, possibly, completely erroneous.

4.2.1. The choice of criteria Eq. (3.6)

as control vector for the problem Eq. (3.6).

initial error eν

(SA) algorithm

shorten the notations).

1 2

ζ′

mated by

κR−1

36 Advances in Statistical Methodologies and Their Application to Real Problems

0 and the observation error vk. The error eν

grated by the direct and adjoint models. But the amplification of vk (and wk when wk=¼ 0) is most worrying since it is integrated in the gradient estimate, making the gradient direction to

The choice of Eq. (3.6) is important in many aspects in order to obtain a simple and efficient data assimilation algorithm. The idea lying in the criteria Eq. (3.6) is to select some pertinent parameters as a control variables for minimizing the mean of the cost function ΨðζkÞ. For

The solution to the problem Eq. (3.6) can be found iteratively using a stochastic optimization

where {ak} is a sequence of positive scalars satisfying some conditions to guarantee a conver-

ak ¼ ∞, ∑ ∞ k¼1 a2

The algorithm Eq. (4.4) is much more simple [compared to the computation of Eq. (4.3)] since it requires, at the kth assimilation instant, to compute only the gradient of the sample cost function ΨðζkÞ. The gradient ∇θΨðθkÞ of the sample objective function ΨðθkÞ can be computed using the AE approach (in what follows, for simplicity, the subscript k will be omitted to

Thus, minimization of Eq. (3.6) by gradient-based SA algorithm requires only one integration of the direct model and one backward integration of the AE code: direct integration of x^k for producing the forecast x^kþ1=k ¼ Φx^k and backward integration ΦTHkζk−1 in computation of

kþ1. For the structure of the gain Eq. (3.9), the objective function Ψ is quadratic wrt θ; hence,

A less computational burden can be achieved by measuring the sample objective function (but not based on a gradient formula): instead of computing the gradient by Eq. (4.6) based on AE, one can approximate the gradient using the values of the cost function [on the basis of finite difference scheme (FDSA)]. Traditionally, the ith component of the gradient can be approxi-

example, for the class of filters Eqs. (3.3), (3.7)–(3.10) the vector θ ¼ ðθ1,…; θne Þ

ak ! 0, ∑ ∞ k¼1

gence of the estimation procedure. The standard conditions are

½δΨðζkþ1Þ�θk ¼ −ðHΦPrδΛKeζk, ζkþ1Þ ¼ −ðδΛKeζk, ζ′

one can find easily the optimal parameters.

κ ðHkμk þ vkÞ. The larger the k, the bigger the amplification of the

θkþ1 ¼ θk−ak∇θΨðζkþ1Þ (4.4)

kþ1Þ, ζ′

kþ1 :¼ PT

k < ∞ (4.5)

r ΦTHTζkþ1 (4.6)

0 is amplified doubly since it is inte-

T can be chosen

As shown in Ref. [12], span½Pr�–the subspace, spanned by the columns of Pr, must be chosen so that the filter gain K ensures a stability of the filter. Mention that even the KF may suffer from instability. To ensure filter stability, Pr is constructed from all unstable and neutral eigenvectors of the fundamental matrix Φ (or real Schur vectors (ScVs), singular vectors). In practice, we choose Pr to be consisting of the column vectors of S (called S-PE samples)

$$S = \Phi X \tag{4.10}$$

which are results of integration of leading ScVs (columns of X). The columns in S have the meaning of the PE for the system state and are used to approximate the ECM M. As to the ScVs, they are preferred to eigenvectors (or singular vectors) because the ScVs are real, and their computation is numerically stable. Mention that computation of singular vector requires also adjoint code. The ensemble of columns of S plays the same role as an ensemble PE samples in the ensemble-based filtering technique for approximating the background ECM [16, 17].

#### 4.2.3. On separation of vertical and horizontal variables structure in ECM [18]

Let us consider the situation when the DS is described by PDEs. The state vector at the time instant k is xk ¼ xkði, j, lÞ where ði, j, lÞ represents a grid point in three dimensional space. Introduce the stabilizing structure for the filter gain [12]

$$\begin{aligned} \mathbf{K}\_k &= \mathbf{M}\_k \mathbf{H}\_k^T [\mathbf{H}\_k \mathbf{M}\_k \mathbf{H}\_k^T + \mathbf{R}\_k]^{-1}, \\ \mathbf{M} &= \mathbf{M}\_d, \mathbf{M}\_d := \mathbf{P}\_r \boldsymbol{\Lambda} \mathbf{P}\_r^T, \end{aligned} \tag{4.11}$$

where Λ is symmetric positive definitive. As shown in Ref. [12], one can choose Λ to be diagonal with diagonal elements serving to regularize the amplitude of ECM. The matrix M in Eq. (4.11) is ECM, and if it is computed on the basis of the Riccati equation (3.6), Kk is known as the KF gain. For M ¼ Md, computation of M is realizable if the reduced dimension ne is not too large. Actually the number ne of the ensemble size is of order Oð100Þ which is too small for M to be a good approximation for the true ECM. In Ref. [18], it is assumed that the estimated ECM is a member of the class of ECMs with separation of vertical and horizontal variables structure (SeVHS). Mention that this hypothesis is not new and used in modeling the ECM in meteorological data assimilation [19]. The optimal ECM is found as a solution of the minimization problem

$$J(\boldsymbol{\theta}) = E \| |M\_d \mkern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1mu \skern-1.1$$

As the number of vertical layers in the today's numerical models is of order Oð10Þ, all elements of the vertical ECM Mv (included in θ1) can be considered as tuning to be estimated. As to Mh, it is often chosen in analytical form (e.g., the 1st or 2nd order autoregressive models). The parameters like correlation length can be selected as components of the control vector θ2 in Mh.

Using dominant real Schur vectors has advantages that they are real and their computation is stable [12] while the computation of eigenvectors is unstable, and they may be complex.

#### 4.3. Computational comparison between VM and AF

We give a brief comparison (computational burden) between the VM and AF algorithms. Table 1 shows the number of elementary arithmetic operations required for implementation of VM and AF filtering algorithms based on AE tool. A smaller number of operations are required for the AF if the SPSA method is used (no need of ΦT ). Here, for simplicity, n2

operations are accounted for the product ΦTy (the same number as that required for Φy). In the AF, we assume that the full ECM M is used, whereas in the AROF (adaptive ROF), M :¼ PrPT r as shown in Eq. (4.6). These numbers are calculated on the basis of Eqs. (2.1) and (2.3) since they represent the most computational burdens for these two algorithms. These numbers are rounded up to the dimensions of the entry matrices. Here, Nito is a number of iterations required to solve the minimization problem (2.1)–(2.3); Nit is a number of iterations required to solve the equation

$$
\Xi y = \zeta\_k,\\
\Xi := [HMH + R]. \tag{4.13}
$$

In the VM algorithm, the computation of M−1 0 , R−1 is not taken into account. In Table 1, there is also the number of operations required for the AROF, when the ECM M is given in the product decomposition form M ¼ PrPT r , Pr∈Rn · ne . In this situation, instead of n2 operations in the AF, we need to perform 2nne operations. For ne << n, much less computational and memory requirements are needed to perform the AROF.

Table 1. Number of elementary arithmetic operations.

ScVs, they are preferred to eigenvectors (or singular vectors) because the ScVs are real, and their computation is numerically stable. Mention that computation of singular vector requires also adjoint code. The ensemble of columns of S plays the same role as an ensemble PE samples in the ensemble-based filtering technique for approximating the background

Let us consider the situation when the DS is described by PDEs. The state vector at the time instant k is xk ¼ xkði, j, lÞ where ði, j, lÞ represents a grid point in three dimensional space.

k ½HkMkHT

where Λ is symmetric positive definitive. As shown in Ref. [12], one can choose Λ to be diagonal with diagonal elements serving to regularize the amplitude of ECM. The matrix M in Eq. (4.11) is ECM, and if it is computed on the basis of the Riccati equation (3.6), Kk is known as the KF gain. For M ¼ Md, computation of M is realizable if the reduced dimension ne is not too large. Actually the number ne of the ensemble size is of order Oð100Þ which is too small for M to be a good approximation for the true ECM. In Ref. [18], it is assumed that the estimated ECM is a member of the class of ECMs with separation of vertical and horizontal variables structure (SeVHS). Mention that this hypothesis is not new and used in modeling the ECM in meteorological data assimilation [19]. The optimal ECM is found as a solution of the minimi-

M ¼ Md, Md :¼ PrΛPT

JðθÞ ¼ EjjMd−Mvðθ1Þ⊗Mhðθ2Þjj2

jj:jjF denotes matrix Frobenious norm , JðθÞ ! min θ, θ ¼ ðθT

As the number of vertical layers in the today's numerical models is of order Oð10Þ, all elements of the vertical ECM Mv (included in θ1) can be considered as tuning to be estimated. As to Mh, it is often chosen in analytical form (e.g., the 1st or 2nd order autoregressive models). The parameters like correlation length can be selected as components of the control vector θ2 in Mh. Using dominant real Schur vectors has advantages that they are real and their computation is stable [12] while the computation of eigenvectors is unstable, and they may be complex.

We give a brief comparison (computational burden) between the VM and AF algorithms. Table 1 shows the number of elementary arithmetic operations required for implementation of VM and AF filtering algorithms based on AE tool. A smaller number of operations are required for the AF if the SPSA method is used (no need of ΦT ). Here, for simplicity, n2

F,

1 , θT 2 Þ T:

k þ Rk� −1 ,

r , (4.11)

(4.12)

4.2.3. On separation of vertical and horizontal variables structure in ECM [18]

Kk ¼ MkHT

Introduce the stabilizing structure for the filter gain [12]

38 Advances in Statistical Methodologies and Their Application to Real Problems

4.3. Computational comparison between VM and AF

ECM [16, 17].

zation problem

To have the idea on how work in practice the VM and AF, here the examples of experiments with two numerical models MICOM (see Section 7) and HYCOM (Hybrid Coordinate Ocean Model) [20] developed at SHOM, Toulouse, France. The first experiment is performed with the MICOM model. The observations are available each 10 days (ds) during 2 years. For the MICOM (state dimension n ¼ 3 · 105 Þ, the 10 ds forecast requires 45 s (supercomputer Caparmor, IFREMER, France, sequential run). The 2-year integration takes 54.45 min (73 · 45 s). The AF needs 54:45 min · 3 ¼ 164 min (or 2 h 45 min) to perform the assimilation experiment for the 2-year period. In this context, the VM requires between 5.7 ds and 11.4 ds to perform the experiment (hypothesis 50–100 times integration of the MICOM over the 2-year window, Section 2.2). As to the HYCOM (state dimension n ¼ 7 · 107 Þ, the observations are available each 5 ds. The 5 ds forecast requires 1 h (supercomputer Beaufix, Météo, France, parallel run, 62 processors). Two year integration requires 146 · 1 h = 146 h (or 6 ds). The AF needs, hence 18 ds to make the 2-year experiment. As to the VM, the experiment requires between 304 ds and 608 ds. That is one of the reasons why in operational setting one has to choose a short window for assimilating the observations by the VM.

Comment 4.1. Looking at Table 1, one sees that the dominant numbers nd of operations in the VM and AF are ndðVMÞ ¼ n2N2 Nito and ndðNAFÞ ¼ n2NNit. If we assume that Nito≈Nit (in fact, the number of iterations for solving the optimization problem Eqs. (2.7) and (2.8) is often larger than the number of iterations for solving the system of equations (4.13); it is more critical for the VM when the DS is nonlinear); the number ndðVMÞ is N times larger than ndðNAFÞ.

## 5. Simple numerical examples: scalar case

#### 5.1. Estimation problem

Consider the simple scalar dynamical process xðtÞ ¼ sinðtÞ. As x\_ðtÞ ¼ cosðtÞ, for tk ¼ kδt, using the approximation xðtk þ δtÞ−xðtkÞ ¼ cosðtkÞδt, one has the following discrete dynamical system

$$\begin{aligned} \mathbf{x}\_{k+1} &= \phi\_k \mathbf{x}\_k + \boldsymbol{\mu}\_k + \boldsymbol{\sigma} \mathbf{w}\_k, \phi\_k = \phi = 1, \boldsymbol{\mu}\_k = \cos(k \boldsymbol{\delta} t), \\ \boldsymbol{E}(\boldsymbol{w}\_k) &= \mathbf{0}, \boldsymbol{E}(\boldsymbol{w}\_k \boldsymbol{\sigma}\_l) = \sigma\_w^2 \delta\_{kl}, k = \mathbf{0}, 1, \dots, \text{N-1}. \end{aligned} \tag{5.1}$$

In Eq. (5.1), wk represents the Gaussian model error. Suppose at each tk moment, we observe the state xk corrupted by the Gaussian noise vk, hence

$$\begin{aligned} \mathbf{z}\_k &= h\mathbf{x}\_k + \mathbf{v}\_k, h = \mathbf{1}, k = 1, \dots, N, \\ E(\mathbf{v}\_k) &= \mathbf{0}, E(\mathbf{v}\_k \mathbf{v}\_l) = \sigma\_v^2 \delta\_{kl}, k \neq l \end{aligned} \tag{5.2}$$

where δkl is the Kronecker symbol. Suppose that the true initial state x� 0 ¼ 0:5. The problem we study here is to estimate the system state xk based on the set of observations zk, k ¼ 1,…; N.

#### 5.1.1. Experiment: cost functions

In the experiment, δt ¼ 0:01, N ¼ 1000. The two methods, VM and AF, will be implemented to produce the system estimates. We study two situations: (S1) the model is considered as perfect, that is, wk ¼ 0; and (S2) there exists the model error wk with the variance σ2 w ¼ 0:001. As to vk, σ2 v ¼ 0:1 in both cases.

Let wk ¼ 0. To see the advantages of the AF over VM in finding optimal solutions, Figures 1 and 2 display the curves (time averaged variance of distances between the true trajectory and those resulting from varying the control variables) as functions of tuning parameters in these two methods: the initial state θ :¼ x0 and the parameter θ :¼ λ (see Eq. (3.9)). We remark that for the filtering system (5.1), (5.2), as φ ¼ 1, the system has one stable eigenvalue and one stable eigenvector (singular value and Schur decompositions are of the same structure). The filter fundamental matrix L ¼ ð1−KhÞφ ¼ ð1−KÞ is stable if K∈ð0; 2Þ. For the gain structure (4.11), L is stable for any M > 0 since K ¼ M Mþσ2 r for h ¼ 1 and σ2 v ¼ 0:1 > 0. We have then K∈ð0; 1Þ⊂ð0; 2Þ. Thus, one can choose θ :¼ M > 0 as tuning parameter. This structure is of less interest compared to Eq. (3.9) since θ enters in K in a nonlinear way, and in fact, K is allowed to vary only in the interval ð0; 1Þ.  For this structure, the filter is stable for K∈ð0; 2Þ. We will select the last structure as a departure point to optimize the AF performance.

From Figure 1, it is seen that the curve "noise-free" is equal to 0 when x^0 ¼ 0:5 for S1, but the "noisy" attains the minimal value 0.121 at x^0 ¼ 0:6. We note that almost the same picture is obtained for the cost function (2.1) (time averaged variance of the distance between x^k, k ¼ 1;…; 1000 and observations zk, k ¼ 1;…; 1000) subject to x^0∈½−1 : 1�. Despite the fact that the curves in Figure 1 are quadratic, it is impossible to find the true initial state in the noisy situation since almost the same curves are obtained for the cost function Eq. (2.1).

Figure 1. VM: cost functions resulting from perfect model and that with a model error.

5. Simple numerical examples: scalar case

40 Advances in Statistical Methodologies and Their Application to Real Problems

the state xk corrupted by the Gaussian noise vk, hence

Consider the simple scalar dynamical process xðtÞ ¼ sinðtÞ. As x\_ðtÞ ¼ cosðtÞ, for tk ¼ kδt, using the approximation xðtk þ δtÞ−xðtkÞ ¼ cosðtkÞδt, one has the following discrete dynamical system

xkþ1 ¼ φkxk þ uk þ wk, φk ¼ φ ¼ 1, uk ¼ cosðkδtÞ,

In Eq. (5.1), wk represents the Gaussian model error. Suppose at each tk moment, we observe

zk ¼ hxk þ vk, h ¼ 1, k ¼ 1,…; N,

study here is to estimate the system state xk based on the set of observations zk, k ¼ 1,…; N.

In the experiment, δt ¼ 0:01, N ¼ 1000. The two methods, VM and AF, will be implemented to produce the system estimates. We study two situations: (S1) the model is considered as perfect,

Let wk ¼ 0. To see the advantages of the AF over VM in finding optimal solutions, Figures 1 and 2 display the curves (time averaged variance of distances between the true trajectory and those resulting from varying the control variables) as functions of tuning parameters in these two methods: the initial state θ :¼ x0 and the parameter θ :¼ λ (see Eq. (3.9)). We remark that for the filtering system (5.1), (5.2), as φ ¼ 1, the system has one stable eigenvalue and one stable eigenvector (singular value and Schur decompositions are of the same structure). The filter fundamental matrix L ¼ ð1−KhÞφ ¼ ð1−KÞ is stable if K∈ð0; 2Þ. For the gain structure (4.11), L is

for h ¼ 1 and σ2

Thus, one can choose θ :¼ M > 0 as tuning parameter. This structure is of less interest compared to Eq. (3.9) since θ enters in K in a nonlinear way, and in fact, K is allowed to vary only in

satisfies Eq. (3.10). For this structure, the filter is stable for K∈ð0; 2Þ. We will select the last

From Figure 1, it is seen that the curve "noise-free" is equal to 0 when x^0 ¼ 0:5 for S1, but the "noisy" attains the minimal value 0.121 at x^0 ¼ 0:6. We note that almost the same picture is obtained for the cost function (2.1) (time averaged variance of the distance between x^k, k ¼ 1;…; 1000 and observations zk, k ¼ 1;…; 1000) subject to x^0∈½−1 : 1�. Despite the fact

EðvkÞ ¼ 0, EðvkvlÞ ¼ σ2

where δkl is the Kronecker symbol. Suppose that the true initial state x�

that is, wk ¼ 0; and (S2) there exists the model error wk with the variance σ2

Mþσ2 r

the interval ð0; 1Þ. For Eq. (3.9), Pr ¼ 1, He ¼ 1 hence Ke ¼ 1

structure as a departure point to optimize the AF performance.

wδkl, k ¼ 0, 1,…; N−1: (5.1)

vδkl, k ≠ l (5.2)

v ¼ 0:1 > 0. We have then K∈ð0; 1Þ⊂ð0; 2Þ.

< 1 and the filter is stable if θ

1þσ2 r

0 ¼ 0:5. The problem we

w ¼ 0:001. As to vk,

EðwkÞ ¼ 0, EðwkvlÞ ¼ σ2

5.1. Estimation problem

5.1.1. Experiment: cost functions

stable for any M > 0 since K ¼ M

v ¼ 0:1 in both cases.

σ2

Figure 2 presents the same curves as those in Figure 1 resulting from application of the filter by letting the parameter θ in the gain vary in θ∈ð0 : 2Þ. Figure 2 shows that for the both curves "noise-free" and "noisy", the minimal values are attained at θ ¼ 1:1 for both situations S1, S2. Moreover, the two minimal values are identical. The same picture is observed for the cost function Eq. (4.1) as function of θ∈ð0 : 2Þ. It means that independently on whether the model is perfect or not, AF formulation allows optimization algorithms to find the optimal value for θ, hence, to ensure optimality of the filter.

Two curves "noisy" in Figures 1 and 2 show that when the model is noisy, the minimal value of the curve "noisy" (VM) in Figure 1 is much higher (it is equal to 0.121) than that in Figure 2 (0.009, for filtering). This fact is in favour of the choice of a short window for assimilating observations by the VM.

Figure 2. Filtering: cost functions resulting from perfect model and that with model error.

## 6. Numerical experiment: Lorenz system

#### 6.1. Lorenz equations

The Lorenz attractor is a chaotic map, noted for its butterfly shape. The map shows how the state of a dynamical system evolves over time in a complex, non-repeating pattern [21].

The attractor itself and the equations from which it is derived were introduced by Edward Lorenz [21], who derived it from the simplified equations of convection rolls arising in the equations of the atmosphere.

The equations that govern the Lorenz attractor are:

$$
\frac{dy\_1}{dt} = -\sigma(y\_1 - y\_2), \\
\frac{dy\_2}{dt} = \rho y\_1 - y\_2 - y\_1 y\_3, \\
\frac{dy\_3}{dt} = y\_1 y\_2 - \beta y\_3,\tag{6.1}
$$

where σ is called the Prandtl number, and ρ is called the Rayleigh number. All σ, β, ρ > 0, but usually σ = 10, β = 8/3 and ρ is varied. The system exhibits chaotic behavior for ρ = 28 but displays knotted periodic orbits for other values of ρ.

#### 6.2. Numerical model

In the experiments, the parameters σ, ρ, β are chosen to have the values 10, 28, and 8/3 for which the "butterfly" attractor exists.

The numerical model is obtained by applying the Euler method (first-order accurate method) to approximate Eq. (6.1). Symbolically, we have

$$y(t\_{k+1}) = F\left(y(t\_k)\right), \newline y(t\_k) := \left(y\_1(t\_k), y\_2(t\_k), y\_3(t\_k)\right)^T,\tag{6.2}$$

where δt :¼ tkþ1−tk is the model time step. The observations arrive at the moments Tk and ΔTk :¼ Tkþ1−Tk. The experiment setup is similar to that described in Ref. [22].

#### 6.3. Observations: assimilation

The corresponding δt ¼ 0:005, ΔTk ¼ 1, hence the sequence of observations is given by zðkÞ :¼ zðTkÞ, k ¼ 1,…; No. The dynamical system corresponding to the transition of the states between two time instants Tk and Tkþ1 is denoted as

$$\mathbf{x}\_{k+1} = F(\mathbf{x}\_k) + \mathbf{w}\_k \tag{6.3}$$

In Eq. (6.3), wk simulates the model error. The sequence wk is assumed to be a white noise having variance 2, 12.13 and 12.13 respectively. The observation system is then given by

$$\mathbf{z}\_{k} = \mathbf{H}\mathbf{x}\_{k} + \mathbf{v}\_{k} \tag{6.4}$$

where the operator H ¼ ½hT 1 , hT 2 � T, h1 ¼ ð1; 0; 0Þ, h2 ¼ ð0; 0; 1Þ, that is, the first and third components x1, x3 are observed at each time instant k ¼ 1;…; 100. The noise sequence vk is white with zero mean and variance R ¼ 2I2 where In is the unit matrix of dimension n. The initial estimate in all filters is given by the initial condition x^ð0Þ¼ð1, −1,24Þ T.

The true system state x� is modeled as the solution of Eq. (6.3) subject to x� 0 ¼ ð1:508870, −1:531271, 25:46091Þ T.

The problem considered in this experiment is to apply the extended KF (EKF), nonadaptive filter (NAF), and adaptive filter (AF) to estimate the true system state using the observations zk, k ¼ 1, 2; No and to compare their performances.

Here, the NAF is in fact the prediction error filter (PEF). Mention that the PEF is developed in Ref. [23] in which the prediction error ECM is estimated on the basis of an ensemble of PE samples, that is,

$$M = \frac{1}{T - 1} \sum\_{k=1}^{T} B\_s(k), B\_s(k) = \sum\_{k=1}^{L} \delta \mathbf{x}\_k^{(l)} \delta \mathbf{x}\_k^{(l), T},\tag{6.5}$$

where δx ðlÞ k , l ¼ 1,…; L are members of the set of L S-PE samples obtained by L þ 1 integrations of the model from the reference state and L perturbed states which grow in the directions of the L dominant Schur vectors.

The filter gain is taken in the following form

6. Numerical experiment: Lorenz system

42 Advances in Statistical Methodologies and Their Application to Real Problems

The equations that govern the Lorenz attractor are:

displays knotted periodic orbits for other values of ρ.

dt ¼ −σðy1−y2Þ,

dy2

dy1

The Lorenz attractor is a chaotic map, noted for its butterfly shape. The map shows how the state of a dynamical system evolves over time in a complex, non-repeating pattern [21].

The attractor itself and the equations from which it is derived were introduced by Edward Lorenz [21], who derived it from the simplified equations of convection rolls arising in the

dt ¼ ρy1−y2−y1y3,

where σ is called the Prandtl number, and ρ is called the Rayleigh number. All σ, β, ρ > 0, but usually σ = 10, β = 8/3 and ρ is varied. The system exhibits chaotic behavior for ρ = 28 but

In the experiments, the parameters σ, ρ, β are chosen to have the values 10, 28, and 8/3 for

The numerical model is obtained by applying the Euler method (first-order accurate method)

where δt :¼ tkþ1−tk is the model time step. The observations arrive at the moments Tk and

The corresponding δt ¼ 0:005, ΔTk ¼ 1, hence the sequence of observations is given by zðkÞ :¼ zðTkÞ, k ¼ 1,…; No. The dynamical system corresponding to the transition of the states

In Eq. (6.3), wk simulates the model error. The sequence wk is assumed to be a white noise having variance 2, 12.13 and 12.13 respectively. The observation system is then given by

nents x1, x3 are observed at each time instant k ¼ 1;…; 100. The noise sequence vk is white with

ΔTk :¼ Tkþ1−Tk. The experiment setup is similar to that described in Ref. [22].

dy3

ðyðtkÞÞ, yðtkÞ :¼ ðy1ðtkÞ, y2ðtkÞ, y3ðtkÞÞT, (6.2)

xkþ1 ¼ FðxkÞ þ wk (6.3)

zk ¼ Hxk þ vk (6.4)

T, h1 ¼ ð1; 0; 0Þ, h2 ¼ ð0; 0; 1Þ, that is, the first and third compo-

dt ¼ y1y2−βy3, (6.1)

6.1. Lorenz equations

equations of the atmosphere.

6.2. Numerical model

which the "butterfly" attractor exists.

6.3. Observations: assimilation

where the operator H ¼ ½hT

to approximate Eq. (6.1). Symbolically, we have

yðtkþ1Þ ¼ F′

between two time instants Tk and Tkþ1 is denoted as

1 , hT 2 �

$$K = MH^T \Sigma^{-1}, \Sigma = HMH^T + R \tag{6.6}$$

which is time invariant. At the same time, for the comparison purpose, the EKF is also used for assimilating the observations.

Figure 3 shows the evolution of the prediction errors resulting from three filters: NAF, AF, and EKF. It is of no surprise that the NAF has produced the estimates with larger estimation error. By adaptation, however, it is possible to obtain the AF, which improves significantly the PEF and even behaves better than the EKF.

Mention that the VM is much less appropriate for assimilating the observations in the Lorenz model due to the choice of the initial state as control vector. For simplicity, we simulate the situation when all three components of the system state are observed in additive noise, i.e. with H ¼ I3. Figure 4 displays time averaged variances of the difference between the true trajectory and model trajectory (denoted as AVðx�, x^ÞÞ, resulting from varying the third component of the initial state for two situations of noise-free and noisy models. Namely, we initialize the model by the initial state, which is the same as the true one x� 0 ¼ ð1:508870,−1:531271, 25:46091Þ T, with the difference, that the third component x^3ð0Þ is varying in the interval ½24:5 : 26:5�. The global minimum is attained at x� 0ð3Þ ¼ 25:46091 as expected. However, if the system is initialized by the estimate in a vicinity, even not so far from x� 3ð0Þ, there is no guarantee that the VM can approach the true initial condition. For the noisy model, the global minimum is not attained at x� 0ð3Þ. As for the PEF, the function AVðx�, x^Þ is quadratic wrt to the gain parameter, for both situations of noise-free or noisy models, as seen in Figure 5: here, the sample cost function (4.1) is computed over all assimilation period, by varying the third parameter θ3 in the gain (related to the third observed component of the system state).

Figure 3. Prediction errors resulting from three filters: nonadaptive filter (NAF), adaptive filter (AF), and EKF.

Figure 4. Time averaged variance between the true trajectory and model trajectory in the VM as a function of perturbed third component of the initial state. The global minimum is attained at the true initial condition, but there is no guarantee for the VM to approach the true initial state (no-noisy model). For noisy model, the global minimum is not attained at the true initial state. The curve "noisy model" is scaled by the factor C ¼ 1=15.

Figure 5. Cost function in the PEF as a function of perturbed third gain parameter θ3. It is seen that in the PEF, the cost function is quadratic wrt to the gain parameter in both situations of noise-free and noisy models. The curve "noisy model" is scaled by the factor C ¼ 1=50.

## 7. Assimilation in high-dimensional model

#### 7.1. MICOM model and assimilation problem

In this section, we show how the AF can be designed in a simple way to produce the high performance estimates for the ocean state in the high-dimensional ocean model MICOM. For details on the Miami Isopycnal Coordinate Ocean Model (MICOM) used here, see Ref. [24]. The model configuration is a domain situated in the North Atlantic from 30°N to 60°N and 80°W to 44°W; for the exact model domain and some main features of the ocean current (mean, variability of the SSH, velocity) produced by the model, see Ref. [24]. The grid spacing is about 0.2° in longitude and in latitude, requiring Nh = II · JJ = 25200 (II = 140, JJ = 180) horizontal grid points. The number of layers in the model is KK = 4. It is configured in a flat bottom rectangular basin (1860 km · 2380 km · 5 km) driven by a periodic wind forcing. The model relies on one prognostic equation for each component of the horizontal velocity field and one equation for mass conservation per layer. We note that the state of the model is x :¼ ðh, u, vÞ where h ¼ hði, j, lrÞ is the thickness of lrth layer; u ¼ uði, j, lrÞ, v ¼ vði, j, lrÞ are two velocity components. The layer stratification is made in the isopycnal coordinates, that is, the layer is characterized by a constant potential density of water. The model is integrated from the state of rest during 20 years. Averaging the sequence of states over 2 years 17 and 18 gives a so-called climatology. During the period of 2 years 19 and 20, every 10 days (10 ds), we calculate the sea surface height (SSH) from the layer thickness h which will serve as a source for generating observations to be used in the assimilation experiments (in total, there are 72 observations).

## 7.2. Different filters

Figure 4. Time averaged variance between the true trajectory and model trajectory in the VM as a function of perturbed third component of the initial state. The global minimum is attained at the true initial condition, but there is no guarantee for the VM to approach the true initial state (no-noisy model). For noisy model, the global minimum is not attained at the

Figure 3. Prediction errors resulting from three filters: nonadaptive filter (NAF), adaptive filter (AF), and EKF.

Figure 5. Cost function in the PEF as a function of perturbed third gain parameter θ3. It is seen that in the PEF, the cost function is quadratic wrt to the gain parameter in both situations of noise-free and noisy models. The curve "noisy model"

true initial state. The curve "noisy model" is scaled by the factor C ¼ 1=15.

44 Advances in Statistical Methodologies and Their Application to Real Problems

is scaled by the factor C ¼ 1=50.

The filter used for assimilating SSH observations is of the form

$$\hat{\mathbf{x}}\_{k+1} = F[\hat{\mathbf{x}}\_k] + \mathbf{K}\mathbb{zeta}\_{k+1}, \boldsymbol{k} = \mathbf{0}, \mathbf{1}, \dots \tag{7.1}$$

where x^kþ1 is the filtered estimate for xkþ1, xkþ1 ¼ ½hkþ1, ukþ1, vkþ1� is the system state at ðk þ 1Þ assimilation instant, Fð:Þ represents the integration of the nonlinear MICOM model over 10 days, K is the filter gain, ζkþ1 is the innovation vector. The gain K is of the form (4.11) where the ECM M will be estimated from the MICOM model. In the experiment, to be closed to realistic situations, only SSH at the grid points i ¼ 1;…; 140, j ¼ 1; ::; 180 are collected as observations. Thus, the observations are available not at all model grid points. The gain K is symbolically written as K ¼ ðKh, Ku, KvÞ T with Ku, Kv representing the operators which produce the correction for the velocity ðu, vÞ from the layer thickness correction Khζ using the geostrophy hypothesis. The filter thus is a reduced order which has the gain Kh to be estimated from S-PE samples.

#### 7.2.1. PEF: computation of ECM

In the experiment, two assimilation methods will be implemented. First the PEF is designed. To do that, the data ECM Md (see Eqs. (6.5) and (7.4), below) is performed by generating an ensemble of PE samples (as done in the experiment with Lorenz system, see the sampling procedure in Ref. [23] for more detail). As the number of elements of ECM is of order 1010 (for only the layer thickness component h), it is impossible to simulate a sufficient number of PE samples so that Md would be a good estimate for the ECM. The matrix Md will be used only as data to estimate the parameters in a parametrized ECM as follows (see Ref. [18]):

Let M∈Rnh · nh be the ECM for the layer thickness h, that is, M ¼ Mðs,s′ Þ. One useful and efficient way to simplify the filter structure is to assume that the ECM M has a SeVHS. Assuming there exist two covariance matrices, Mv and Mh such that

$$M(\mathbf{s}, \mathbf{s}') = M\_v(\mathbf{s}\_v, \mathbf{s}\_v) \otimes M\_h(\mathbf{s}\_h, \mathbf{s}\_h; \cdot), \mathbf{s}\_v := l, \mathbf{s}\_h := (i, j), \tag{7.2}$$

where ⊗ denotes the Kronecker product between two matrices [25],

$$M\_{\upsilon}(s\_{\upsilon}, s\_{\upsilon}) \otimes M\_{\hbar}(s\_{\imath}, s\_{\imath'}) = M(i, j, l; \mathbf{i}', \mathbf{j}', l') = $$

$$\begin{pmatrix} m\_{\upsilon}(1, 1)M\_{\hbar} & m\_{\upsilon}(1, 2)M\_{\hbar} & \dots & m\_{\upsilon}(1, n\_{\upsilon})M\_{\hbar} \\\\ m\_{\upsilon}(2, 1)M\_{\hbar} & m\_{\upsilon}(2, 2)M\_{\hbar} & \dots & m\_{\upsilon}(2, n\_{\upsilon})M\_{\hbar} \\\\ \dots & \dots & \dots & \dots \\\\ m\_{\upsilon}(n\_{\upsilon}, 1)M\_{\hbar} & m\_{\upsilon}(n\_{\upsilon}, 2)M\_{\hbar} & \dots & m\_{\upsilon}(n\_{\upsilon}, n\_{\upsilon})M\_{\hbar} \end{pmatrix} \tag{7.3}$$

The main advantage of the separability hypothesis is that the number of parameters to be estimated in the covariance matrix is reduced drastically. As a consequence, even an ensemble of PE samples with small size can serve as a large data set for estimating the unknown parameters. This results in a fast convergence of the estimation procedure. In addition, introducing the SeVHS hypothesis allows to avoid the rank deficiency problem in the estimation of the ECM. In fact, as only a few numbers of ScVs can be computed in very high-dimensional systems, approximation of the ECM M by Eq. (6.5) results in rank deficiency for M. With such an ECM, the resulting filter will probably produce worse results, not to say on instability which may occur during the filtering process.

Suppose we are given the ensemble of S-PE samples Sτ½L�¼½δhð1Þ τ ,…; δhðLÞ τ � which are obtained at the τ time instant by applying the sampling procedure in Ref. [23] subject to L perturbations. For τ ¼ 1,…; T, the ECM Md in Eq. (4.12) is estimated as

$$\begin{aligned} M\_d &= \frac{1}{T} \sum\_{\tau=1}^T M\_{\tau}, \\\\ M\_{\tau} &:= \frac{1}{L} \mathcal{S}\_{\overline{\tau}}(L) \mathcal{S}\_{\overline{\tau}}^T(L) \end{aligned} \tag{7.4}$$

For the problem Eqs. (4.3)–(4.4) let us define the vector of unknown parameters in the ECM Mðs,s′ Þ as

A Comparison Study on Performance of an Adaptive Filter with Other Estimation Methods for State... http://dx.doi.org/10.5772/67005 47

$$\Theta := (\mathbf{c}\_{11}, \dots, \mathbf{c}\_{1n\_v}, \mathbf{c}\_{21}, \dots, \mathbf{c}\_{2n\_v}, \mathbf{c}\_{n\_l 1}, \dots, \mathbf{c}\_{n\_v n\_v}, L\_d)^T, \mathbf{c}\_{kl} := m\_v(k, l). \tag{7.5}$$

where ckl :¼ mvðk, lÞ. As to the parameter Ld, it represents the correlation length of the horizontal ECM Mh,

$$M\_h(y, y') = \exp(-d(y, y')/L\_d),\tag{7.6}$$

dðy, y′ Þ is the distance between two horizontal points y :¼ ði, jÞ and y′ :¼ ði ′ , j ′ Þ.

Considering Md as data matrix, optimization problem for determining the vector θ looks like

$$J[\theta] = E[\Psi(M\_t, \theta)] \to \min \theta,\tag{7.7}$$

$$[\Psi(M\_t, \theta)] := ||M\_t - M\_v(\mathbf{s}\_v, \mathbf{s}\_v) \otimes M\_h(\mathbf{s}\_h, \mathbf{s}\_h)||\_F^2,\tag{7.7}$$

Mention that the problem Eq. (7.4) is somewhat closely related to the Nearest Kronecker Product (NKP) problem [25]. In the experiment, the correlation length Ld is not estimated and is taken identical in two filters PEF and CHF, Ld ¼ 25.

#### 7.2.2. PEF: computation of gain

ensemble of PE samples (as done in the experiment with Lorenz system, see the sampling procedure in Ref. [23] for more detail). As the number of elements of ECM is of order 1010 (for only the layer thickness component h), it is impossible to simulate a sufficient number of PE samples so that Md would be a good estimate for the ECM. The matrix Md will be used only as

efficient way to simplify the filter structure is to assume that the ECM M has a SeVHS.

Mvðsv,sv′Þ⊗Mhðsh,sh′Þ ¼ Mði, j, l; i

mvð1; 1ÞMh mvð1; 2ÞMh … mvð1, nvÞMh mvð2; 1ÞMh mvð2; 2ÞMh … mvð2, nvÞMh … …… …

mvðnv, 1ÞMh mvðnv, 2ÞMh … mvðnv, nvÞMh

The main advantage of the separability hypothesis is that the number of parameters to be estimated in the covariance matrix is reduced drastically. As a consequence, even an ensemble of PE samples with small size can serve as a large data set for estimating the unknown parameters. This results in a fast convergence of the estimation procedure. In addition, introducing the SeVHS hypothesis allows to avoid the rank deficiency problem in the estimation of the ECM. In fact, as only a few numbers of ScVs can be computed in very high-dimensional systems, approximation of the ECM M by Eq. (6.5) results in rank deficiency for M. With such an ECM, the resulting filter will probably produce worse results, not to say on instability

at the τ time instant by applying the sampling procedure in Ref. [23] subject to L perturbations.

L SτðLÞST

For the problem Eqs. (4.3)–(4.4) let us define the vector of unknown parameters in the ECM

τ ðLÞ

Md ¼ 1 T ∑ T τ¼1 Mτ,

Mτ :¼ 1

Þ ¼ Mvðsv,sv′Þ⊗Mhðsh,sh′Þ,sv :¼ l,sh :¼ ði, jÞ, (7.2)

′ , j ′ , l ′ Þ ¼

1

CCCCCCA

τ ,…; δhðLÞ

τ � which are obtained

(7.4)

Þ. One useful and

(7.3)

data to estimate the parameters in a parametrized ECM as follows (see Ref. [18]):

Let M∈Rnh · nh be the ECM for the layer thickness h, that is, M ¼ Mðs,s′

Assuming there exist two covariance matrices, Mv and Mh such that

where ⊗ denotes the Kronecker product between two matrices [25],

Mðs,s ′

46 Advances in Statistical Methodologies and Their Application to Real Problems

0

BBBBBB@

which may occur during the filtering process.

Mðs,s′

Þ as

Suppose we are given the ensemble of S-PE samples Sτ½L�¼½δhð1Þ

For τ ¼ 1,…; T, the ECM Md in Eq. (4.12) is estimated as

As to the computation of the gain, introduce the notations: at the instant k, let xði, j, lÞ be the value of the system state defined at the grid points ði, j, lÞ. Let x ! ¼ ðx !T 1 , x !T 2 ,…; x !T nv Þ T be a vector representation for x where x ! l is a vector whose components are the values of x at all the horizontal grid points (ordered in some way) at the l th, vertical layer.

Consider the ECM Eq. (3.10) and the observation equation (1.2). Represent the observation matrix H in a block-matrix form

$$H = [H\_1, \dots, H\_{n\_v}] \tag{7.8}$$

which corresponds to the vector representation x !, that is,

$$H\vec{\mathbf{x}} = \sum\_{\nu=1}^{n\_v} H\_{\nu} \vec{\mathbf{x}}\_{\nu} \tag{7.9}$$

Compute the gain according to Eq. (2.4). We have

$$\begin{aligned} MH^T &= M\_v \otimes M\_h H^T = \begin{bmatrix} \Sigma\_1^T, \dots, \Sigma\_{n\_v}^T \end{bmatrix}^T, \\ \Sigma\_l &= M\_h \sum\_{k=1}^{n\_v} c\_{lk} H\_k^T = M\_h G\_{v,l}, \\ MH^T &= M\_d G\_v, M\_d = \text{block diag} \, \text{diag} \, [M\_h, \dots, M\_h], \\ G\_v &= [G\_{v,1}^T, \dots, G\_{v,n\_v}^T]^T, G\_{v,l} := \sum\_{k=1}^{n\_v} c\_{lk} H\_k^T, \\ \Sigma &:= HM H^T + R = \sum\_{k=1, l=1}^{n\_v} c\_{lk} H\_l M\_h H\_k^T + R, \end{aligned} \tag{7.10}$$

As proved in Ref. [18], in this case, the gain has the following form

$$K = M\_d K\_v,\\ K\_v = \begin{bmatrix} K\_{v,1}^T, \dots, K\_{v,n\_v}^T \end{bmatrix}^T = G\_v \Sigma^{-1},\\ \tag{7.11}$$

$$K\_{v,l} = G\_{v,l} \Sigma^{-1}, l = 1, \dots, n\_v,$$

where Gv,l, Σ are defined in Eq. (7.9).

#### 7.2.3. Cooper-Haines filter: CHF

To see the performance of the AF, we implement also a so-called CHF. The CHF [7] is obtained from (7.9), (7.10) under three hypotheses [24]: (H1) the analysis error of the system output is canceled in the case of noise-free observations; (H2) conservation of the linear potential vorticity (PV); (H3) there is no correction for the velocity at the bottom layer. The AF in Ref. [24] is obtained by relaxing one or several hypotheses (H1)–(H3). From the filtering theory, the difference between the PEF and CHF is lying in the way we estimate the elements of the ECM.

For the choice of the tuning parameters in the PEF, see Ref. [24].

#### 7.3. Numerical results

First, we run the model initialized by the climatology. This run is different from that used for modeling the sequence of true ocean states only by changing the initial state by climatology. This run is denoted as model. Figure 6 shows the (spatial) averaged variance between SSH observations and that produced by the model. We see the error grows as time progresses, meaning instability of the numerical model wrt the perturbed initial condition. That fact signifies that the VM will have difficulties in producing high performance estimates if the assimilation window is long.

Figure 6. SSH prediction errors resulting from "model" and CHF: growing of the prediction error in the "model" signifies instability of the numerical model wrt specification of the initial system state.

Next the two filters, CHF and PEF, are implemented under the same initial condition as those carried out with the experiment model. It is seen from Figure 6 that initialized by the same initial condition, and the CHF is much more efficient than the model in reducing the estimation error. The performance comparison between the CHF and PEF is presented in Figure 7. Here, the (spatially averaged variance) SSH prediction errors, resulting from two filters CHF and PEF, are displayed. The superiority of the PEF over the CHF is undoubted. It is clear from Figure 7 that the PEF is capable of producing the better estimates, with lower error level, along all assimilation period. On the other hand, if the estimation error in CHF decreases continuously at the beginning of the assimilation and is stabilized more or less after (during the interval k∈ð10 : 45Þ), it becomes to increase considerably at the end of the assimilation period. It means that the PEF is more efficient than the CHF and the fact that the ECM in the PEF, constructed on the basis of the S-PE samples, has the effect to stabilize its behavior.

K ¼ MdKv, Kv ¼ ½KT

48 Advances in Statistical Methodologies and Their Application to Real Problems

For the choice of the tuning parameters in the PEF, see Ref. [24].

instability of the numerical model wrt specification of the initial system state.

where Gv,l, Σ are defined in Eq. (7.9).

7.2.3. Cooper-Haines filter: CHF

7.3. Numerical results

assimilation window is long.

v,1,…; KT

Kv,l ¼ Gv,lΣ−1

To see the performance of the AF, we implement also a so-called CHF. The CHF [7] is obtained from (7.9), (7.10) under three hypotheses [24]: (H1) the analysis error of the system output is canceled in the case of noise-free observations; (H2) conservation of the linear potential vorticity (PV); (H3) there is no correction for the velocity at the bottom layer. The AF in Ref. [24] is obtained by relaxing one or several hypotheses (H1)–(H3). From the filtering theory, the difference between the PEF and CHF is lying in the way we estimate the elements of the ECM.

First, we run the model initialized by the climatology. This run is different from that used for modeling the sequence of true ocean states only by changing the initial state by climatology. This run is denoted as model. Figure 6 shows the (spatial) averaged variance between SSH observations and that produced by the model. We see the error grows as time progresses, meaning instability of the numerical model wrt the perturbed initial condition. That fact signifies that the VM will have difficulties in producing high performance estimates if the

Next the two filters, CHF and PEF, are implemented under the same initial condition as those carried out with the experiment model. It is seen from Figure 6 that initialized by the same initial condition, and the CHF is much more efficient than the model in reducing the estimation error. The performance comparison between the CHF and PEF is presented in Figure 7. Here, the (spatially averaged variance) SSH prediction errors, resulting from two

Figure 6. SSH prediction errors resulting from "model" and CHF: growing of the prediction error in the "model" signifies

v,nv �

T ¼ GvΣ−1

, l ¼ 1,…; nv,

,

(7.11)

Figure 7. SSH prediction errors resulting from two filters: CHF and PEF. The PEF is capable of producing the estimates with lower estimation errors and is stable along all assimilation period, whereas the CHF has a difficulty to maintain the same performance at the end of the assimilation period. The PEF is much better than the CHF in providing better estimates for the system states.

To see the effect of adaptation, Figure 8 displays the filtered errors for the u-velocity component estimates at the surface, produced by the PEF and APEF, respectively. The APEF is an AF, which is an adaptive version of the PEF. Here, the tuning parameters are optimized by the SPSA method. From the computational point of view, the SPSA requires much less time

Figure 8. Filtered errors for the u-velocity component estimate, resulting from PEF and APEF (AF based on PEF). Optimization is performed by SPSA. By tuning the parameters in the filter gain, one can improve considerably the performance of the PEF.

integration and memory storage compared with the traditional AE method. At each assimilation instant, we have to make only two integrations of the MICOM for approximating the gradient vector. From Figure 8, one sees that the adaptation allows to reduce significantly the estimation errors produced by the PEF.

## 8. Conclusions

In this chapter, a comparison study on the performance of the AF with other existing methods is presented in the context of its application to data assimilation problems in high-dimensional numerical models. As it is seen, in comparison with the standard VM, the AF is much simpler to implement and produces better estimates. The advantages of the AF over other methods such as EKF, CHF are also demonstrated. The principal reason for high performance of the AF is lying in the choice of innovation representation for the initial input-output system and selection of pertinent gain parameters as control variables to minimize the MSE of the innovation process. If in the VM, the choice of the structure for the initial state is the most important thing to do (but that is insufficient for guaranteeing its high performance), in the AF, however, the initial state has a little impact on the performance of the AF. This happens because the AF is selected as optimal in the class of stable filters, and as a consequence, the error in the initial estimate is attenuated during assimilation process. In contrary, in the VM, the error in the specification of the initial guess is amplified during assimilation if the numerical model is unstable.

In conclusion, it is important to emphasize that the AF approach, presented in this chapter, is consolidated by exploiting the following road map: (i) generate data ECM from S-PE samples which grow in the directions of dominant Schur vectors; (ii) select a parametrized structure for the ECM under the hypothesis SeVHS; (iii) make the choice of tuning parameters in the gain by minimizing the distance between the data ECM and that having the SeHVS structure; (iv) adjust the unknown parameters in the gain in order to minimize the PE error of the system output by applying the SPSA algorithm.

There are a wide variety of engineering problems to which the AF is applicable and that could be worthy of further study. Depending on particular problems, undoubtedly, the other modifications would be helpful to improve the filter performance, to simplify its implementation. But the main features of the AF, presented in this chapter, remain as the key points to follow in order to preserve a high performance of the AF.

## Author details

Hong Son Hoang\* and Remy Baraille

\*Address all correspondence to<EMAIL_ADDRESS>
SHOM/HOM/REC, Toulouse, France

## References

integration and memory storage compared with the traditional AE method. At each assimilation instant, we have to make only two integrations of the MICOM for approximating the gradient vector. From Figure 8, one sees that the adaptation allows to reduce significantly the

In this chapter, a comparison study on the performance of the AF with other existing methods is presented in the context of its application to data assimilation problems in high-dimensional numerical models. As it is seen, in comparison with the standard VM, the AF is much simpler to implement and produces better estimates. The advantages of the AF over other methods such as EKF, CHF are also demonstrated. The principal reason for high performance of the AF is lying in the choice of innovation representation for the initial input-output system and selection of pertinent gain parameters as control variables to minimize the MSE of the innovation process. If in the VM, the choice of the structure for the initial state is the most important thing to do (but that is insufficient for guaranteeing its high performance), in the AF, however, the initial state has a little impact on the performance of the AF. This happens because the AF is selected as optimal in the class of stable filters, and as a consequence, the error in the initial estimate is attenuated during assimilation process. In contrary, in the VM, the error in the specification of the initial guess is amplified during assimilation if the numerical model is

In conclusion, it is important to emphasize that the AF approach, presented in this chapter, is consolidated by exploiting the following road map: (i) generate data ECM from S-PE samples which grow in the directions of dominant Schur vectors; (ii) select a parametrized structure for the ECM under the hypothesis SeVHS; (iii) make the choice of tuning parameters in the gain by minimizing the distance between the data ECM and that having the SeHVS structure; (iv) adjust the unknown parameters in the gain in order to minimize the

There are a wide variety of engineering problems to which the AF is applicable and that could be worthy of further study. Depending on particular problems, undoubtedly, the other modifications would be helpful to improve the filter performance, to simplify its implementation. But the main features of the AF, presented in this chapter, remain as the key points to follow in

PE error of the system output by applying the SPSA algorithm.

order to preserve a high performance of the AF.

\*Address all correspondence to<EMAIL_ADDRESS>
Hong Son Hoang\* and Remy Baraille

SHOM/HOM/REC, Toulouse, France

estimation errors produced by the PEF.

50 Advances in Statistical Methodologies and Their Application to Real Problems

8. Conclusions

unstable.

Author details

#### **Applications of the H-Principle of Mathematical Modelling** Applications of the H-Principle of Mathematical Modelling

#### Agnar Höskuldsson Agnar Höskuldsson

[17] Hamill T.M.: Ensemble-based atmospheric data assimilation. In: Palmer T., editor. Pre-

[18] Hoang H.S. and Baraille R.: A low cost filter design for state and parameter estimation in very high dimensional systems. In: Proceedings of the 19th IFAC Congress; August, 2014,

[19] Daley R.: The effect of serially correlated observation and model error on atmospheric

[20] Baraille R.: Modélisation numérique avec le modèle HYCOM au SHOM. At "http://

[22] Kivman G.A.: Sequential parameter estimation for stochastic systems. Nonlinear Process

[23] Hoang H.S. and Baraille R.: Prediction error sampling procedure based on dominant Schur decomposition. Application to state estimation in high dimensional oceanic model.

[24] Hoang H.S., Baraille R. and Talagrand O.: On an adaptive filter for altimetric data assimilation and its application to a primitive equation model MICOM. Tellus, 57A,

[25] Golub G.H. and Van Loan C.F.: Matrix Computations. Johns Hopkins University Press;

dictability of Weather and Climate. Cambridge Press; 2006; pp. 124–156.

[21] Lorenz E.N.: Deterministic non-periodic flow. J Atmos Sci. 1963; 20: 130–141.

data assimilation. Monthly Weather Rev. 1992; 120: 165–177

mathocean.math.cnrs.fr/presentations/Baraille.pdf"

52 Advances in Statistical Methodologies and Their Application to Real Problems

J Appl Math Comput. 2011; 12: 3689–3709.

Geophys. 2003; 10: 253–259. doi:10.5194/npg-10-253-2003

pp. 3256–3261.

2005, 2: 153–170.

1996.

Additional information is available at the end of the chapter Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/66153

#### Abstract

Traditional statistical test procedures are briefly reviewed. It is pointed out that significance testing may not always be reliable. The author has formulated a modelling procedure, the H-principle, for how mathematical modelling should be carried out in the case of uncertain data. Here it is applied to linear regression. Using this procedure, the author has developed a common framework for carrying out linear regression. Six regression methods are analysed by this framework, two stepwise methods: principal component regression, ridge regression, PLS regression and an H-method. The same algorithm is used for all methods. It is shown how model validation and graphic analysis, which is popular in chemometrics, apply to all the methods. Validation of the methods is carried out by using numerical measures, cross-validation and test sets. Furthermore, the methods are tested by a blind test, where 40 samples have been excluded. It is shown how procedures in applied statistics and in chemometrics both apply to the present framework.

Keywords: linear regression, common framework for linear regression, stepwise regression, ridge regression, PLS regression, H-methods

## 1. Introduction

Regression analysis is the most studied subject within theoretical statistics. Numerous books have been published. Advanced program packages have been developed that make it easy for users to develop and study advanced models and methods.

Advanced program packages such as SAS and SPSS have been used by students since the 1970s and the program packages have become very advanced. The user fills out a 'menu' similar to those at restaurants. The program then carries out the analysis as requested providing the users' possibility of carrying out highly advanced analysis of the data. Users rely upon that they can make interpretation of the output in the way they learned at school. For example,

© The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons

distribution, and eproduction in any medium, provided the original work is properly cited.

if a variable/factor is significant, its presence in the model improves the predictions derived from the modelling results.

Today measurement instruments are becoming more and more advanced, e.g. optical instruments are becoming popular in industry. They may give thousands of values each time a sample is measured. In applied research, the tendency is to include as much information as possible in order to cover possible alternatives of the situation. When adding interactions or non-linear terms, we also see data having hundreds or thousands of variables.

There are three basic challenges in applied data analysis, which are as follows: (1) Data in applied sciences and industry are typically generated by instruments that produce many variables. In these cases, the data have latent structure, which means that the data values are geometrically located in a low-dimensional space. However, mathematical models and methods usually assume that the X-data have full rank. In these cases, the models are often incorrect and may give imprecise results.

(2) Typically, scientists develop solutions that are optimal or unbiased. This is a natural approach, because the derived solutions have important properties, when data satisfy the given model. Simulations based on the model confirm the optimality, but often in practice the data often do not satisfy the assumptions of the model. Forcing an optimal solution on the data may give results that have bad or no prediction ability. A better solution may then be obtained by relaxing on the optimality or unbiasedness of the solution.

(3) It is a tradition to base results of regression analysis on testing the significance of the variables in the model. However, the results may not always be reliable. The influence of a variable can be so small that it may be of no importance, see Section 3.3. At professional organizations, it is considered to be serious problem for the researchers in the field.

H-principle is a prescription of how a solution to a mathematical model should be generated, when data are uncertain. It is proposed that the determination of the solution should be carried out in steps and compute the improvement in the solution and the associated precision at each step. It is suggested that the solution at each step should be an optimal balance between the improvement and the associated precision. Determination of the solution continues as long as it is supported by the given (uncertain) data.

The author has developed a framework linear regression, which is inspired by the H-principle. Using this framework, the same algorithm is applied for carrying out different types of regression analysis. Most regression methods based on linear algebra can be carried out within this framework. Here, we carry out the analysis of six different methods. Numerical and graphic analyses of the results are the same for all the regression methods.

In Section 2, we specify the used notation and briefly describe the used data. In Section 3, we consider some basic issues in modelling data. The background for the H-principle is treated. In Section 4, we consider the latent variable regression model. In Section 5, we present the Hprinciple and show some examples of its usage. In Section 6, we present a common framework for linear regression. In Section 7, we discuss model validation. This is an important topic, but we only consider essential aspects that are used in the analysis of the six methods. In Section 8, we present the results for six different regression methods: (1) stepwise (forward) regression maximising covariance, (2) stepwise (forward) regression maximising R2 ,

(3) principal component regression, (4) ridge regression,

(5) PLS regression and (6) an H-method.

if a variable/factor is significant, its presence in the model improves the predictions derived

Today measurement instruments are becoming more and more advanced, e.g. optical instruments are becoming popular in industry. They may give thousands of values each time a sample is measured. In applied research, the tendency is to include as much information as possible in order to cover possible alternatives of the situation. When adding interactions or

There are three basic challenges in applied data analysis, which are as follows: (1) Data in applied sciences and industry are typically generated by instruments that produce many variables. In these cases, the data have latent structure, which means that the data values are geometrically located in a low-dimensional space. However, mathematical models and methods usually assume that the X-data have full rank. In these cases, the models are often

(2) Typically, scientists develop solutions that are optimal or unbiased. This is a natural approach, because the derived solutions have important properties, when data satisfy the given model. Simulations based on the model confirm the optimality, but often in practice the data often do not satisfy the assumptions of the model. Forcing an optimal solution on the data may give results that have bad or no prediction ability. A better solution may then be obtained

(3) It is a tradition to base results of regression analysis on testing the significance of the variables in the model. However, the results may not always be reliable. The influence of a variable can be so small that it may be of no importance, see Section 3.3. At professional

H-principle is a prescription of how a solution to a mathematical model should be generated, when data are uncertain. It is proposed that the determination of the solution should be carried out in steps and compute the improvement in the solution and the associated precision at each step. It is suggested that the solution at each step should be an optimal balance between the improvement and the associated precision. Determination of the solution continues as long as

The author has developed a framework linear regression, which is inspired by the H-principle. Using this framework, the same algorithm is applied for carrying out different types of regression analysis. Most regression methods based on linear algebra can be carried out within this framework. Here, we carry out the analysis of six different methods. Numerical and

In Section 2, we specify the used notation and briefly describe the used data. In Section 3, we consider some basic issues in modelling data. The background for the H-principle is treated. In Section 4, we consider the latent variable regression model. In Section 5, we present the Hprinciple and show some examples of its usage. In Section 6, we present a common framework for linear regression. In Section 7, we discuss model validation. This is an important topic, but we only consider essential aspects that are used in the analysis of the six methods. In Section 8,

organizations, it is considered to be serious problem for the researchers in the field.

graphic analyses of the results are the same for all the regression methods.

non-linear terms, we also see data having hundreds or thousands of variables.

from the modelling results.

incorrect and may give imprecise results.

it is supported by the given (uncertain) data.

by relaxing on the optimality or unbiasedness of the solution.

54 Advances in Statistical Methodologies and Their Application to Real Problems

Section 9 briefly discusses the results presented. In Section 10, we mention application of the H-principle to multi-block and path modelling, non-linear modelling and extension to multilinear algebra. Section 11 presents conclusions.

## 2. Notation, data and scaling

Matrices and vectors are denoted by bold letters; matrices by upper case and vectors by lower case. In order to facilitate the reading of the equations, different types of indices are used for the steps and for the matrices/vectors. The letters a and b are related to the steps in the algorithm. The letters i, j and k are used for the indices within a matrix/vector. It is assumed that there is given instrumental data X, N +K matrix, and response data Y, N +M matrix. The regression data are denoted by (X, Y). In some cases, only one y-variable, y, is used. This is done to simplify the equations. It is assumed that data are centred, which means that average values are subtracted from each column of X and Y. This also makes it easier to read the equations. xj is the j th column of X and xi is the i th row of X = (xij). A latent variable τ is a linear combination of the original measured variables,

$$
\pi = a\_1 \mathbf{x}\_1 + a\_2 \mathbf{x}\_2 + \dots + a\_K \mathbf{x}\_K \tag{1}
$$

The data used here are from an optical instrument. The instrument gives 1200 values at each measurement. However, technical knowledge of the instrument suggests that only 40 values should be used for determining the substance in question, the y-values. Two hundred samples are measured. Forty samples are put aside for a blind test. Thus, the calibration analysis is based on 160 samples. This gives X as 160 + 40 matrix and y a 160 vector. These data are challenging and represent a common situation, when working with optical instruments (FTIR, NIR, fluorescence etc.).

It is sometimes important to determine if data should be scaled or not. Scaling can be obtained by multiplying X and Y by diagonal matrices. If C1 and C2 are diagonal matrices, the scaling of X is done by the transformation, X←(XC1) and of Y by Y← (YC2). The linear least squares solution for (X, Y) is given by B ¼ ðXTXÞ −1 XTY. This solution can be obtained from the linear least squares solution for scaled data, B1, as follows,

$$\mathbf{B} = \mathbf{C}\_1[(\mathbf{X}\mathbf{C}\_1)^T(\mathbf{X}\mathbf{C}\_1)]^{-1}(\mathbf{X}\mathbf{C}\_1)^T(\mathbf{Y}\mathbf{C}\_2)]\mathbf{C}\_2^{-1} = \mathbf{C}\_1\mathbf{B}\_1\mathbf{C}\_2^{-1} \tag{2}$$

This shows that when computing the linear least squares solution, we can work with the scaled data. The original solution is obtained by scaling 'back' as shown in the equation. We use this property also, when we compute the approximate solution. The effect of scaling is better numerical precision. Scaling is necessary for the present data. If data are not scaled (e.g. to unit variance), numerical results beyond dimension of around 15 may not be reliable. The reason is that Eqs. (20) and (23) are sensitive to small numerical values. Scaling is much debated among researchers. When scaling is used, one must secure that all variables have values above the 'noise' level of the instrument. This is a difficult topic, which is not considered closer here. If original data follow a normal distribution, the scaled ones will not. However, for large sample number like given here, the differences will be negligible.

#### 3. Linear regression

#### 3.1. Traditional regression model

It is assumed that there are given instrumental data X and response data y. A linear regression model is given by

$$y = \beta\_1 \mathbf{x}\_1 + \beta\_2 \mathbf{x}\_2 + \dots + \beta\_K \mathbf{x}\_K + \varepsilon \tag{3}$$

The x-variables are called independent variables and the y-variable is the dependent one. When the parameters β have been estimated as b, the estimated model is now

$$y = b\_1 \mathbf{x}\_1 + b\_2 \mathbf{x}\_2 + \dots + b\_K \mathbf{x}\_K \tag{4}$$

When there is given a new sample x0 = (x10, x20, …, xK0), it gives the estimated or predicted value y0 = b1 x10 + … + bK xK0. There can be many estimates for the regression coefficients. Therefore, Greek letters are used for the theoretical parameters and Roman letters for the estimated values. It is common to use the linear least squares method for estimating the parameters. It is based on minimizing the residuals, (y–Xβ) T (y–Xβ) ! minimum, with respect to β. The solution is given by

$$\mathbf{b} = (\mathbf{X}^{\mathsf{T}}\mathbf{X})^{-1}\mathbf{X}^{\mathsf{T}}\mathbf{y}.\tag{5}$$

It is common to assume that the residuals in Eq. (3) are normally distributed. This is often written as y~N(Xβ, σ2 ). It means that the expected value of y is E(y) = Xβ and the variance is Var(y) = σ2 I, where I is the identity matrix. The linear least squares procedure coincides with the maximum likelihood method in case the data follow a normal distribution. Assuming the normal distribution, the parameter estimate b has the variance given by

$$Var(b) = \sigma^2 (X^T X)^{-1}.\tag{6}$$

Here σ2 is estimated by

$$
\sigma^2 \mathbb{1} \mathbb{E}\_1^\mathrm{N} e\_i^2 / (\text{N-K}) \tag{7}
$$

where the residuals, ei' s, are computed from e = y – Xb.

#### 3.2. Assumptions and properties

numerical precision. Scaling is necessary for the present data. If data are not scaled (e.g. to unit variance), numerical results beyond dimension of around 15 may not be reliable. The reason is that Eqs. (20) and (23) are sensitive to small numerical values. Scaling is much debated among researchers. When scaling is used, one must secure that all variables have values above the 'noise' level of the instrument. This is a difficult topic, which is not considered closer here. If original data follow a normal distribution, the scaled ones will not. However, for large sample

It is assumed that there are given instrumental data X and response data y. A linear regression

The x-variables are called independent variables and the y-variable is the dependent one.

When there is given a new sample x0 = (x10, x20, …, xK0), it gives the estimated or predicted value y0 = b1 x10 + … + bK xK0. There can be many estimates for the regression coefficients. Therefore, Greek letters are used for the theoretical parameters and Roman letters for the estimated values. It is common to use the linear least squares method for estimating the

b ¼ ðXTXÞ

VarðbÞ ¼ σ2

σ2≅ΣN 1 e 2

s, are computed from e = y – Xb.

−1

I, where I is the identity matrix. The linear least squares procedure coincides with

ðXTXÞ −1

It is common to assume that the residuals in Eq. (3) are normally distributed. This is often

the maximum likelihood method in case the data follow a normal distribution. Assuming the

When the parameters β have been estimated as b, the estimated model is now

y ¼ β1x1 þ β2x2 þ … þ βKxK þ ε (3)

y ¼ b1x1 þ b2x2 þ … þ bKxK (4)

(y–Xβ) ! minimum, with respect

XTy: (5)

: (6)

i =ðN−KÞ (7)

T

). It means that the expected value of y is E(y) = Xβ and the variance is

number like given here, the differences will be negligible.

56 Advances in Statistical Methodologies and Their Application to Real Problems

parameters. It is based on minimizing the residuals, (y–Xβ)

normal distribution, the parameter estimate b has the variance given by

3. Linear regression

model is given by

3.1. Traditional regression model

to β. The solution is given by

written as y~N(Xβ, σ2

Here σ2 is estimated by

where the residuals, ei'

Var(y) = σ2

It is assumed that the matrix XT X has an inverse. Furthermore, it is assumed that the samples are independent of each other. In case the model (Eq. (3)) is correct and data follow the normal distribution, the estimate (Eq. (5)) has some important properties compared with other possible estimates:

There is a general agreement that the variance matrix Var(b) should be as small as possible. The matrix (XT X)–1 , the precision matrix, shows how precise the estimates b are. The properties (a) and (b) state that assuming the linear model (Eq. (3)), the solution (Eq. (5)) is, in fact from a theoretical point of view, the best possible one. It often occurs that the model (Eq. (3)) contains many variables. Interpretation of the estimation results (Eq. (4)) is an important part of the statistical analysis. Therefore, it is common to evaluate the parameters with the aid of a significance test. If a variable is not significant, it may be excluded. A parameter associated with a variable is commonly evaluated by a t-test. A t-test of a parameter bi is given by

$$t = \frac{b\_i}{\sqrt{Var(b\_i)}} \text{ = } \frac{b\_i}{\sqrt{(s^2 \times s^{\bar{\imath}})}} \tag{8}$$

Here s 2 is given by Eq. (7) and s ii is the ith diagonal element of (XT X) –1 .

This is the motivation that the program packages such as SAS and SPSS in statistics compute the linear least squares solution in the case of linear regression analysis as the initial solution. Significance of the parameters in Eq. (4) is shown using Eq. (8).

The problem in using the least squares solution appears when the precision matrix becomes close to singular. If the computer program (SAS or SPSS) detects that the computation of b may be imprecise due to close to singularity of the precision matrix, the user is informed that the estimates b may be imprecise. However, the information is related to the judgement of the situation and the numerical precision of the computer. There are practical issues long before the question of the numerical uncertainties of the estimates, b's.

Consider two examples. Suppose that X is N +2 and that x1 = c x2 + δ, for some constant c, where |δ| < 10–5 . We may be able to compute (XT X)–1 . However, inference on the model (Eq. (4)) may be uncertain or incorrect. At the other example, suppose that X is N +40, but that practical rank is 15. Assume that XT X = X1 T X1 + X2 T X2, where X2 = (x2,1 x2,2 … x2,40). If all x2,i i = 1,…,40 are small, say |x2,i| < 10–5 , inference from Eq. (4) may be uncertain or incorrect. The first example may not be realistic. However, the second one is realistic. It is common that data in applied sciences and industry are of reduced practical rank. In these cases, the model (Eq. (3)) is incorrect and leads to uncertain or incorrect results.

#### 3.3. Stepwise linear regression

We shall consider closer the procedure of stepwise linear regression. We do that by using the Cholesky factorization XT X = FFT , where F is lower triangular. Then we can write the columns of the data matrix X as

$$\begin{aligned} \mathbf{x}\_1 &= F\_{11}\mathbf{t}\_1\\ \mathbf{x}\_2 &= F\_{21}\mathbf{t}\_1 + F\_{22}\mathbf{t}\_2\\ &\dots\\ \mathbf{x}\_K &= F\_{K1}\mathbf{t}\_1 + F\_{K2}\mathbf{t}\_2 + \dots + F\_{KK}\mathbf{t}\_K \end{aligned} \tag{9}$$

The vectors (ti) are mutually orthogonal and of length 1, ti T tj = δij. The significance of xK, when x1, x2, …, xK–1 are given, is computed from Eq. (8) with

$$\mathbf{b}\_{K} = \frac{\mathbf{y}^{T}(\mathbf{F\_{KK}}\mathbf{t\_{K}})}{\mathbf{F\_{KK}^{2}}} = \frac{(\mathbf{y}^{T}\mathbf{t\_{K}})}{F\_{KK}} \text{ and } \mathbf{s^{KK}} = \mathbf{1}/F\_{KK}^{2} \tag{10}$$

This gives

$$\mathbf{t} = \frac{(\mathbf{y}^T \mathbf{t} \boldsymbol{\kappa})}{\mathbf{s}} \tag{11}$$

This shows that the t-test for the significance of xK is independent of the size FKK. When the selection of a new variable is carried out among many variables (e.g. 500), there is a considerable risk that FKK is so small that the marginal effect of xK is of no importance although it is being declared as statistically significant. The issue is that the user of program packages is not informed of, if a significant variable improves the prediction of the response variable or not.

#### 3.4. Variance of regression coefficients

It is instructive to study closer the variance of the regression coefficient in the linear least squares model, Eq. (6). It can be written as

$$Var(\mathbf{b}) \mathbb{E}[\mathbf{y}^T \mathbf{y} - \mathbf{y}^T \mathbf{X} (\mathbf{X}^T \mathbf{X})^{-1} \mathbf{X}^T \mathbf{y}] \times [\left(\mathbf{X}^T \mathbf{X}\right)^{-1}]/(\text{N} - \text{K})\tag{12}$$

An important objective of the modelling task is to keep the value of Eq. (12) as small as possible. Equation (12) is a product of two parts. For these two parts it can be shown that, when a new variable is added to the model

(In theory these measures can be unchanged, but in practice these changes always occur). Note that the two terms, (i) and (ii), are equally important and appear in a symmetric way in Eq. (12). If the data are normally distributed, it can be shown that

are stochastically independent. This means that the knowledge of one of them does not give any information on the other one. If, for instance, a certain fit has been obtained, then there is no information on (b). In order to know the quality of predictions, we must be compute (b) and use Eq. (12). Knowing (a), the results concerning the precision matrix can be good or bad. Program packages only use the t-test or some equivalent measure to evaluate to test significance of the variables/components. There is no information on the precision matrix. We must compute the precision matrix together with the significance testing in order to find out how well the model is performing. However, a modelling procedure must include both (a) and (b) in order to secure small values of Eq. (12). The procedure must handle both the decrease of (i) and the increase of (ii) during the model estimation.

## 4. Latent variable model

3.3. Stepwise linear regression

X = FFT

58 Advances in Statistical Methodologies and Their Application to Real Problems

The vectors (ti) are mutually orthogonal and of length 1, ti

bK ¼ yTðFKKtKÞ F2 KK

VarðbÞ≅½yTy−yTXðXTXÞ

−1

Eq. (12). If the data are normally distributed, it can be shown that

−1 , −1

−1 XTy�,

, always increases

x1, x2, …, xK–1 are given, is computed from Eq. (8) with

3.4. Variance of regression coefficients

squares model, Eq. (6). It can be written as

when a new variable is added to the model

i. the residual error, ½yTy−yTXðXTXÞ

a. the residual error, ½yTy−yTXðXTXÞ

b. the precision matrix, ðXTXÞ

ii. the precision matrix, ðXTXÞ

x1 ¼ F11t1 x2 ¼ F21t1 þ F22t2

Cholesky factorization XT

of the data matrix X as

This gives

We shall consider closer the procedure of stepwise linear regression. We do that by using the

… xK ¼ FK1t1 þ FK2t2 þ … þ FKKtK

> ¼ ðyTtKÞ FKK

t ¼ ðyTtKÞ

This shows that the t-test for the significance of xK is independent of the size FKK. When the selection of a new variable is carried out among many variables (e.g. 500), there is a considerable risk that FKK is so small that the marginal effect of xK is of no importance although it is being declared as statistically significant. The issue is that the user of program packages is not informed of, if a significant variable improves the prediction of the response variable or not.

It is instructive to study closer the variance of the regression coefficient in the linear least

−1

An important objective of the modelling task is to keep the value of Eq. (12) as small as possible. Equation (12) is a product of two parts. For these two parts it can be shown that,

(In theory these measures can be unchanged, but in practice these changes always occur). Note that the two terms, (i) and (ii), are equally important and appear in a symmetric way in

XTy� · ½ðXTXÞ

XTy�, always decreases

−1

, where F is lower triangular. Then we can write the columns

tj = δij. The significance of xK, when

s (11)

KK (10)

�=ðN−KÞ (12)

T

and sKK ¼ 1=F2

(9)

The measurement data in applied sciences and industry often represent a 'system', where variables are dependent on each other from chemical equilibrium, action-reaction in physics and physical and technical balance. Geometrically, the samples are located in a low-dimensional space. In these cases, the model (Eq. (3)) is incorrect (apart from simple, designed laboratory experiments) and the nice theory above may not applicable. Why is this the case? When a new sample is available, we cannot automatically insert it into Eq. (4) and compute the y-value. We must secure that the new sample is located geometrically along the samples used for analysis, the rows of X. A correct model for this kind of data is

$$y = a\_1\tau\_1 + a\_2\tau\_2 + \dots + a\_A\tau\_A + \varepsilon.\tag{13}$$

Here (τa) are latent variables and (αa) are the regression coefficients on the latent variables. The regression task is both to determine the latent variables, τ's, and to compute the associated regression coefficients. The resulting Eq. (13) is then converted into Eq. (4) using Eq. (26). The model (Eq. (13)) is called latent structure linear regression. It is sometimes difficult to make 'good' interpretation of the variables, when using latent variables. However, by using appropriate modelling procedure, it is possible to obtain a fairly good interpretation of the variables (see Section 5.4).

## 5. The H-principle of mathematical modelling

#### 5.1. Background

In the 1920s, there were large discussions on the measurement aspects of quantum mechanics. W. Heisenberg pointed out that there are certain magnitudes that cannot be determined exactly at the same time. His famous uncertainty inequality states that there is a lower limit to how well conjugate magnitudes can be determined at the same time. The position and momentum of an elementary particle is an example. For that example, the inequality is

$$
\Delta(\text{position}) \times \Delta(\text{momentum}) \gtrsim \text{constant} \tag{14}
$$

The lower limit of the inequality is related to the Planck's constant of light. These considerations are based on physical theory. In practice, it means that there are some restrictions on the outcome of an experiment. The results may depend on the instrument and the phenomenon being studied. The restrictions are detected by the application of the measurement instrument. The uncertainty inequality and the associated theory are a kind of guidance, when carrying out experiments.

## 5.1.1. The H-principle

When modelling data, there is an analogous situation. Instead of a measurement instrument we have a mathematical method. The conjugate magnitudes are ΔFit and ΔPrecision and they cannot be controlled at the same time. It is recommended to carry out the modelling in steps and at each step evaluate the situation as prescribed by the uncertainty inequality. At each step, it is assumed that there is given a weight vector. The recommendations of the H-principle of mathematical modelling are the following (expressed in the case of regression analysis):

	- a. expression for improvement in fit, Δ(Fit)
	- b. the associated prediction, Δ(Precision)

Note that the H-principle is a recommendation of how to proceed in determining the solution of a mathematical problem, when data are uncertain. It is not like maximum likelihood, where a certain function is to be maximized. However, it is necessary to compute the improvement in the solution at each step and the associated precision. The optimal balance is described in Section 3 and the situation is evaluated to find out if data are following along.

#### 5.1.2. Application to linear regression

Let us consider closer how this applies to linear regression. The task is to determine a weight vector w according to this principle. For the score vector t = Xw we have

If S = XT X the adjustment (Eq. (20)) gives orthogonal score vectors, t's. Therefore, the variance, (b), is used for the precision. Treating Σ as a constant the task is to maximize

Applications of the H-Principle of Mathematical Modelling http://dx.doi.org/10.5772/66153 61

$$
\left[\frac{|\mathbf{Y}^T\mathbf{t}|^2}{\mathbf{t}^T\mathbf{t}}\right] \times \frac{I}{\left[\frac{1}{\mathbf{t}^T\mathbf{t}}\right]} = \mathbf{w}^T\mathbf{X}^T\mathbf{Y}\mathbf{Y}^T\mathbf{X}\mathbf{w} \tag{15}
$$

This is a maximization task because improvement in fit is negative. The maximization is carried out under the restriction that w is of length 1, |w|=1. By using the Lagrange multiplier method, it can be shown that the task is an eigenvalue task,

$$\lambda X^T Y Y^T X \mathfrak{w} = \lambda \mathfrak{w} \tag{16}$$

In case there is only one y-variable, the eigenvalue task has a direct solution

$$w = \frac{X^T y}{|X^T y|} = \frac{\mathbf{C}}{|\mathbf{C}|} \tag{17}$$

These are the solutions used in PLS regression. Thus, we can state that the PLS regression is consistent with the H-principle. H-methods take as a starting point the PLS solution and determine how the prediction aspect or the solution can be improved.

#### 5.2. Interpretation of the modelling results

The fit obtained by the score vector t is

$$\frac{\left(\mathbf{y}^{\mathsf{T}}\mathbf{t}\right)^{2}}{\left(\mathbf{t}^{\mathsf{T}}\mathbf{t}\right)} = \frac{\left(\mathbf{C}^{\mathsf{T}}\mathbf{C}\right)^{2}}{\mathbf{C}^{\mathsf{T}}\mathbf{S}\mathbf{C}} = \left[\frac{c\_{1}^{2}}{f} + \frac{c\_{2}^{2}}{f} + \dots + \frac{c\_{K}^{2}}{f}\right]^{2} \tag{18}$$

where f ¼ ffiffiffiffiffiffiffiffiffiffiffiffi CTSC p

ΔðpositionÞ ·ΔðmomentumÞ ≥ constant (14)

+

Δ(Precision)

The lower limit of the inequality is related to the Planck's constant of light. These considerations are based on physical theory. In practice, it means that there are some restrictions on the outcome of an experiment. The results may depend on the instrument and the phenomenon being studied. The restrictions are detected by the application of the measurement instrument. The uncertainty inequality and the associated theory are a kind of guidance, when carrying out

When modelling data, there is an analogous situation. Instead of a measurement instrument we have a mathematical method. The conjugate magnitudes are ΔFit and ΔPrecision and they cannot be controlled at the same time. It is recommended to carry out the modelling in steps and at each step evaluate the situation as prescribed by the uncertainty inequality. At each step, it is assumed that there is given a weight vector. The recommendations of the H-principle of mathematical modelling are the following (expressed in the case of regression analysis):

1. Carry out determining the solution in steps. You specify how you want to look at the data

4. In case the computed solution improves the prediction abilities of the model, the solution is accepted. If the solution does not provide this improvement, the modelling stops.

Note that the H-principle is a recommendation of how to proceed in determining the solution of a mathematical problem, when data are uncertain. It is not like maximum likelihood, where a certain function is to be maximized. However, it is necessary to compute the improvement in the solution at each step and the associated precision. The optimal balance is described in

Let us consider closer how this applies to linear regression. The task is to determine a weight

X the adjustment (Eq. (20)) gives orthogonal score vectors, t's. Therefore, the variance,

at this step by formulating how the weights are computed.

60 Advances in Statistical Methodologies and Their Application to Real Problems

a. expression for improvement in fit, Δ(Fit)

3. Compute the solution that minimises the product Δ(Fit)

5. The data are adjusted for what has been selected; restart at 1.

Section 3 and the situation is evaluated to find out if data are following along.

vector w according to this principle. For the score vector t = Xw we have

(b), is used for the precision. Treating Σ as a constant the task is to maximize

2 =ðt TtÞ

TtÞ

b. the associated prediction, Δ(Precision)

experiments.

5.1.1. The H-principle

2. At each step compute

5.1.2. Application to linear regression

a. Improvement in fit: jYTtj

b. Associated variance: Σ=ðt

If S = XT

The regression coefficient can be written similarly as

$$\frac{(\mathbf{y}^T \mathbf{t})}{(\mathbf{t}^T \mathbf{t})} = \frac{(\mathbf{C}^T \mathbf{C})^{1.5}}{\mathbf{C}^T \mathbf{S} \mathbf{C}} = \left[ \frac{c\_1^2}{\mathcal{g}} + \frac{c\_2^2}{\mathcal{g}} + \dots + \frac{c\_K^2}{\mathcal{g}} \right]^{1.5} \tag{19}$$

where ¼ ðCTSCÞ 1 1:5.

At each step, we can see how much each variable contributes to the fit and the regression coefficient..

#### 6. A common framework for linear regression

#### 6.1. Background: views on the regression analysis task

There are many ways to carry out a regression analysis. Here we present a general framework for carrying out linear regression that includes most methods based on linear algebra. The basic idea is to separate the computations into two parts: the first part is concerned on how it is preferable to look at data. In practice, there can be different emphasis on the regression analysis. Sometimes it is important to determine important variables. In other cases, it may be the predictions derived that is important. The other part is the numerical algorithm to compute the solution vector and associated measures. At the first part, it is assumed that a weight matrix W = (w1, …,wK) is given. Each weight vector should be of length one, although this is not necessary. The weight vectors (wa) specify how one wants to look at the data. They may be determined by some optimization or significance criteria. They may also be determined by a standard regression method such as PLS regression. In Section 6.4, we show the choices in the case of the six regression analysis. The role of the weight vectors is only to compute the loading vectors. The other part of the computations is a numerical algorithm, which is the same for all choices of W. There is no restriction on the weight vectors except that they may not give zeroloading vector, pa ≠ 0.

#### 6.2. Decomposition of data

The starting point is a variance matrix S and a covariance matrix C. S can be any positive semidefinite matrix, but C is assumed to be the covariance, C = XT Y. The algorithm is formulated for multiple y's.

Initially S0 = S and C0 = C and B = 0. For a = 1, …, K:

$$\text{Loading\\_vector}: \mathbf{p}\_a = \mathbf{S}\_{a-1} \mathbf{w}\_a \tag{20}$$

$$\text{Scaling} \quad \text{constant}: d\_d = \frac{1}{\mathbf{w}\_a^T \mathbf{p}\_a} \tag{21}$$

$$\text{Loading vector}: \mathbf{q}\_a = \mathbf{C}\_{a-1}^T \mathbf{w}\_a \tag{22}$$

$$\text{Loading weight vector}: \mathfrak{v}\_a \tag{23}$$

The loading weight vectors are computed by, see reference [1],

$$\mathbf{w}\_1 = \mathbf{w}\_1, \mathbf{v}\_a = \mathbf{w}\_a - d\_1(\mathbf{p}\_1^T \mathbf{w}\_a)\mathbf{v}\_1 - \dots - d\_{a-1}(\mathbf{p}\_{a-1}^T \mathbf{w}\_a)\mathbf{v}\_{a-1}, a = \text{2, 3, ..., K} \tag{24}$$

S is adjusted by the loading vector pa and similarly for C,

$$\mathbf{S}\_{a} = \mathbf{S}\_{a-1} \mathbf{-} d\_{a} \mathbf{p}\_{a} \mathbf{p}\_{a}^{T} \tag{25}$$

$$\mathbf{C}\_{a} = \mathbf{C}\_{a-1} \mathbf{-} d\_{a} \mathbf{p}\_{a} \mathbf{q}\_{a}^{T} \tag{26}$$

The adjustment of S is of rank one reduction. S also reduces in size by

$$tr(d\_a \mathbf{p}\_a \mathbf{p}\_a^T) = d\_a \mathbf{p}\_a^T \mathbf{p}\_a = \frac{\mathbf{w}\_a^T \mathbf{S}^2 \mathbf{w}\_a}{\mathbf{w}\_a^T \mathbf{S} \mathbf{w}\_a} > 0. \tag{27}$$

The loading weight matrix V satisfies VT P = D–1 , see reference [1].

#### 6.3. Expansion of matrices

analysis. Sometimes it is important to determine important variables. In other cases, it may be the predictions derived that is important. The other part is the numerical algorithm to compute the solution vector and associated measures. At the first part, it is assumed that a weight matrix W = (w1, …,wK) is given. Each weight vector should be of length one, although this is not necessary. The weight vectors (wa) specify how one wants to look at the data. They may be determined by some optimization or significance criteria. They may also be determined by a standard regression method such as PLS regression. In Section 6.4, we show the choices in the case of the six regression analysis. The role of the weight vectors is only to compute the loading vectors. The other part of the computations is a numerical algorithm, which is the same for all choices of W. There is no restriction on the weight vectors except that they may not give zero-

The starting point is a variance matrix S and a covariance matrix C. S can be any positive semi-

Scaling constant : da ¼ 1

Loading vector : qa ¼ CT

1waÞv1−…−da−1ðpT

Sa ¼ Sa−1−dapapT

Ca ¼ Ca−1−dapaqT

a pa ¼ wT

a S2 wa

, see reference [1].

wT a Swa

a Þ ¼ dapT

P = D–1

Y. The algorithm is formulated for

a−1wa (22)

a−1waÞva−1, a ¼ 2, 3, …, K (24)

a (25)

a (26)

> 0: (27)

(21)

Loading vector : pa ¼ Sa−1wa (20)

Loading weight vector : va (23)

wT a pa

definite matrix, but C is assumed to be the covariance, C = XT

62 Advances in Statistical Methodologies and Their Application to Real Problems

The loading weight vectors are computed by, see reference [1],

v1 ¼ w1, va ¼ wa−d1ðpT

S is adjusted by the loading vector pa and similarly for C,

The loading weight matrix V satisfies VT

The adjustment of S is of rank one reduction. S also reduces in size by

trðdapapT

Initially S0 = S and C0 = C and B = 0. For a = 1, …, K:

loading vector, pa ≠ 0.

multiple y's.

6.2. Decomposition of data

At a = K we get SK = 0. Expanding S and C we get

$$\mathbf{S} = d\_1 \mathbf{p}\_1 \mathbf{p}\_1^T + \dots + d\_A \mathbf{p}\_A \mathbf{p}\_A^T + \dots + d\_K \mathbf{p}\_K \mathbf{p}\_K^T = \mathbf{P} \mathbf{D} \mathbf{P}^T \tag{28}$$

$$\mathbf{C} = \mathbf{X}^T \mathbf{Y} = d\_1 \mathbf{p}\_1 \mathbf{q}\_1^T + \dots + d\_A \mathbf{p}\_A \mathbf{q}\_A^T + \dots + d\_K \mathbf{p}\_K \mathbf{q}\_K^T = \mathbf{P} \mathbf{D} \mathbf{Q}^T \tag{29}$$

Inserting appropriate matrices we get

$$\mathbf{S}^{-1} = d\_1 \mathbf{v}\_1 \mathbf{v}\_1^T + \dots + d\_A \mathbf{v}\_A \mathbf{v}\_A^T + \dots + d\_K \mathbf{v}\_K \mathbf{v}\_K^T = \mathbf{V} \mathbf{D} \mathbf{V}^T \tag{30}$$

$$\mathbf{B} = \mathbf{S}^{-1}X^TY = d\_1\mathbf{v}\_1\mathbf{q}\_1^T + \dots + d\_A\mathbf{v}\_A\mathbf{q}\_A^T + \dots + d\_K\mathbf{v}\_K\mathbf{q}\_K^T = \mathbf{V}\mathbf{D}\mathbf{Q}^T\tag{31}$$