# Optimal Modeling and Filtering of Stochastic Time Series for Geoscience Applications

Authors: J. A. Rod Blais
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

Sequences of observations or measurements are often modeled as realizations of stochastic processes with some stationary properties in the first and second moments. However in practice, the noise biases and variances are likely to be different for different epochs in time or regions in space, and hence such stationarity assumptions are often questionable. In the case of strict stationarity with equally spaced data, the Wiener-Hopf equations can readily be solved with fast Fourier transforms (FFTs) with optimal computational efficiency. In more general contexts, covariance matrices can also be diagonalized using the Karhunen-Loève transforms (KLTs), or more generally using empirical orthogonal and biorthogonal expansions, which are unfortunately much more demanding in terms of computational efforts. In cases with increment stationarity, the mathematical modeling can be modified and generalized covariances can be used with some computational advantages. The general nonlinear solution methodology is also briefly overviewed with the practical limitations. These different formulations are discussed with special emphasis on the spectral properties of covariance matrices and illustrated with some numerical examples. General recommendations are included for practical geoscience applications.


---

## Body

## 1. Introduction

Consider a stochastic process or time series{X(t),-∞<t<∞} which has been observed or measured as {y(s),-∞<s<∞}, most often in practice at equispaced discrete time or space intervals Δs=1 for simplicity over some finite domain. In practice, it is also often convenient to assume these processes or time series to be centered, that is, with null mean value or expectation following the appropriate trend modeling in increment stationary situations. When these observations (or measurements) are of Gaussian processes, their second moments or covariances then fully specify the processes. Otherwise, it can only be assumed that the covariance functions are sufficient for the following modeling and filtering operations. Specifically, it is assumed that these covariance functions have positive Fourier transforms or power spectra for realizability.In general, the optimal estimate ofX(t) from observations y(s) is given by the conditional expectation
(1)X^(t)=E[X(t)∣y(s)]
in which the conditional expectation of a random variable U given an observation v is defined by
(2)E[U∣v]=∫-∞∞up(u∣v)du
using the conditional density p(u∣v); that is,
(3)p(u∣v)=p(u,v)pU(u)
for the joint density p(u,v) with positive marginal density pU(u); that is,
(4)pU(u)=∫-∞∞p(u,v)dv.
These quantities can readily be generalized to several variables. The previous expression is also called the conditional mean with corresponding conditional variance
(5)Var[U∣v]=E[(U∣v-E[U∣v])2]≡E[(U∣v-E[U∣v])(U∣v-E[U∣v])*]
which is commonly used in optimization, with the * corresponding to conjugate transpose, or simply transpose in real cases. For the random variable V corresponding to the observations, E[U∣V] and Var[U∣V] are then the conditional expectation for the mean and variance of U given V.Such conditional expectationE[X(t)∣Y(s)] can easily be seen as an optimal estimate X^(t) as
(6)E[X^(t)]=E[E[X(t)∣Y(s)]]=E[X(t)]
and its variance
(7)Var[X^(t)]=Var[E[X(t)∣Y(s)]]=Var[X(t)]-E[Var[X(t)∣Y(s)]]≤Var[X(t)]
as expected. For any other estimator Z(t), that is, any function of Y(t), it is straightforward to verify that
(8)MSE[X(t)-E[X(t)∣Y(s)]]≤MSE[X(t)-Z(t)]
for the mean-square errors (MSEs); that is, the conditional expectation estimator E[X(t)∣Y(s)] is optimal in the mean-square sense. For any unbiased estimator Z(t), the relationship is in the least-squares sense with variances replacing the mean-square errors. More discussion of these relationships can be found, for example, in [1].In the linear context, the conditional expectationX^(t)=E[X(t)∣y(s)] is a linear function of the observations y(s); that is,
(9)X^(t)=∫-∞∞H(t,s)y(s)ds
for some unknown kernel H(t,s), or with discrete data,
(10)X^(t)=∑s=-∞∞H(t,s)y(s),
that is, simply an unknown linear combination of the available data. In practice, only finite data are available so that the previous limits are finite, and furthermore in causal cases, only past observations are available for processing. Furthermore to simplify the presentation in the following, the inner product notation will be used; that is,
(11)X^(t)=〈H(t,s),y*(s)〉
for the continuous and discrete applications whenever appropriate. The superscript * stands for conjugate transpose although only real applications are considered herein.The objective in the following is to estimate optimallyX(t) from the available data y(s) using minimal assumptions. Using the above discussed property of the conditional expectation estimator E[X(t)∣y(s)], one can write the implied orthogonality of the residuals to the observations; that is,
(12)E[(X(t)-X^(t))y*(s)]=E[〈(X(t)-X^(t)),y(s)〉]=0;
that is, using the linear transformation of y(t) for X^(t),
(13)E[〈(X(t)-〈H(t,r),y*(r)〉),y(s)〉]=0
or
(14)E[〈X(t),y(s)〉]=E[〈〈H(t,r),y*(r)〉,y(s)〉]=E[〈H(t,r),〈y(r),y(s)〉〉]=〈H(t,r),E[〈y(r),y(s)〉]〉
which is usually interpreted in terms of cross-covariance between X(t) and y(s) being equal to the autocovariance of y(r) modified by the unknown kernel H(t,r), assuming E[y(r)]=0, or equivalently,
(15)〈Kyy(s,r),H(r,t)〉=KyX(s,t)
using the cross-covariance KyX(s,t) and (symmetric real positive definite) covariance Kyy(s,r), often written as Ky(s,r) in practice, respectively. These normal equations are the well-known Wiener-Hopf equations, which are Fredholm integral equations (of the first kind) in the continuous context. In cases of stationarity, this unknown kernel H(r,t) becomes H(r-t) which will be seen to greatly simplify the situation in practice in terms of computations.It should be mentioned that the previous interchange of the order in the linear expectation operation in the inner products needs to be done with care in the case of continuous applications with infinite limits. However for discrete applications with finite limits, the situation is much simpler and usually carried out without second thoughts. More discussion of those aspects can be found in, for example, the works of Dougherty [1] and Pugachev [2] with references in functional analysis.The general nonlinear solution methodology will also briefly be overviewed in the Gaussian context with some of the practical limitations. These different formulations are discussed with emphasis on the properties of the covariance matrices and illustrated with some numerical examples and references. General recommendations are included for practical geoscience applications.


## 2. Fourier Transform Approaches

The Wiener-Hopf equations written as(16)〈Kyy(s,r),H(r,t)〉=KyX(s,t)
for unknown H(r,t) are trivial for white noise, that is, with Kyy(s,r)≡δ(s,r)≡δ(s-r), the Dirac delta function, implying that H(r,t)≡H(r-t)=KyX(r,t)≡KyX(r-t) in such cases. However, when the white noise has different intensities in different parts of the domain, this simple formulation cannot be used and the subsequent sections will deal with the more general situation.When the covariance kernelKyy(s,r) can be written as Kyy(s-r) meaning that the covariance is only dependent on the spacing between its arguments, then it can be diagonalized using a Fourier transform. The corresponding covariance matrix is Toeplitz; that is, its diagonals are constant in terms of s-r for row index s and column index r with gridded data. Notice that the discrete Fourier transform (DFT) of a finite Toeplitz matrix is not necessarily diagonal unless the matrix is circulant as will be illustrated in Section 6.Explicitly using a Fourier matrixF≡[frs∣frs=e2πirs/N] of order N, the DFT of a vector v is simply ℱ[v]≡v∧≡F*v, and hence in the case of stationarity,
(17)F*Kyy(s-r)FF*H(r)=NF*KyX(r)
or
(18)Kyy∧(τ)H∧(τ)=NKyX∧(τ)
for frequencies τ; that is,
(19)H∧(τ)=N·KyX∧(τ)Kyy∧(τ)
obviously for Kyy∧(τ)≠0 assuming the symmetry and positive definiteness of Kyy, in which F·F*=F*·F=N. This means that the problem reduces to estimating frequency components of the y(s) observations to infer the corresponding frequency components of the X(t) process. Furthermore, in general, the transfer kernel or function H∧(τ) is noncausal, as for each time t, it uses all of the information from -∞<t<∞ in the Fourier transform. Hence in practical cases with only finite or even only past information, complications can be expected (see, e.g., the discussion in [3]).


## 3. Empirical Orthogonal Function Expansions

The covariance kernelKyy(s,r) for the observations is generally assumed to be symmetric and positive definite for physical realizability. Its eigenfunctions are the exponential functions only in the stationary case as F*Kyy(s,r)Fis then diagonal and F*Kyy(s,r)F=diag(λν) or Kyy(s,r)F=N-1Fdiag(λν) for the corresponding eigenvalues λν. The discrete Fourier transform application or orthogonal projection onto the exponential eigenfunctions can be generalized to other systems of orthogonal functions in the sense of Galerkin’s orthogonal projection methodology for such functional operator systems. This will be illustrated in Section 6 using Haar wavelet transforms.Hence for a spectrum of orthogonal eigenfunctionsuk(r) and corresponding eigenvalues λk, usually ordered as λ1≥λ2≥λ3≥⋯>0,
(20)Kyy(s,r)=∑k=1∞λkuk(s)uk*(r)
can be used for an expansion of the random functions. The eigenfunction expansions of random functions or Karhunen-Loève transformed random variables then have a diagonal covariance kernel or matrix. Notice that, in finite dimensions, the expansion of the covariance matrix is simply the eigenvalue expansion which is a special application of the Singular Value Decomposition (SVD) (see, e.g., the works of Klema and Laub [4], Lewis [5], and Parlett [6]). The eigenfunction expansion of the covariance kernel is often referred to as Mercer’s theorem in the mathematical literature. The eigenfunctions derived from the data covariance are also called empirical orthogonal functions (EOFs) and have proven most effective in signal data processing and in modeling dynamical systems such as in oceanography and related fields (see, e.g., Kim and Wu [7], North and Cahalan [8]).The mathematical situation is summarized by the well-knownKarhunen-Loève theorem: Let X(t) be a zero-mean random function on T with a weight function w(t)≥0 such that its covariance kernel satisfies
(21)∫T∫T|KX(s,r)|2w(r)dsdr<∞;
then (i)
for the integral equation(22)∫TKX(s,r)u(r)w(r)dr=λu(s),
there exist eigenvalues λ1≥λ2≥λ3≥⋯>0 with corresponding eigenfunctions u1(s),u2(s),u3(s),… such that
(23)∫TKX(s,r)uk(r)w(r)dr=λkuk(s)
for all k such that {uk(t)} is a deterministic orthonormal system on T relative to the weight function w(t) meaning that
(24)∫Tuk(t)uh*(t)w(t)dt=δkh;(ii)
the corresponding generalized Fourier coefficients ofX(t) relative to these eigenfunctions
(25)ξk=∫TX(t)uk*(t)w(t)dt
are uncorrelated and Var[ξk]=λk;(iii)
X(t) can then be represented as an eigenfunction expansion as
(26)X(t)=∑k=1∞ξkuk(t)
in the mean-square sense.This last orthogonal expansion of the random function X(t) is known in the literature as a generalized Fourier expansion, a Karhunen-Loève expansion, or a canonical expansion. Also, the Fourier and Karhunen-Loève expansions are usually introduced in the functional analysis context of square-integrable functions such as L2(ℝ) while the canonical expansion context is simply a decomposition of the covariance matrix or kernel (see, e.g., the works of Dougherty [1] for more discussion).In practice when only the observational sequence{y(s)} is known, the SVD can be used to obtain the eigenvectors of the corresponding covariance function without explicitly evaluating the covariance matrix. As different data matrices are generally rectangular in shape, different strategies are advisable for the SVD computations (see, e.g., the works of Björck [9] and Blais [10]). Furthermore, especially in multiple dimensions, the SVD and Principal Component Analysis (PCA) terminologies are not always consistent in the applied science and engineering literature largely because of related data array conventions and objectives.


## 4. Empirical Biorthogonal Expansions

Rather than solving the eigenvalue equation for the orthonormal eigenfunctions{uk(t)} with corresponding eigenvalues {λk}, let us assume that an arbitrary sequence {ak(t)} of square-integrable functions are available. Then a biorthogonal sequence {vk(t)} can be generated to satisfy the conditions:(i)
xk(t)=〈KX(t,s),vk(s)〉/〈〈KX(t,s),vk(t)〉,vk(s)〉,(ii)
〈xk(t),vh(t)〉=δkhrecursively in a well-known manner using the following projections: (i)
setv1(t)=a1(t) with x1(t)=〈KX(t,s),v1(s)〉/〈〈KX(t,s),v1(t)〉,v1(s)〉,(ii)
definev2(t)=c21v1(t)+a2(t) with c21 determined from 〈x1(t),v2(t)〉=0; that is,
(27)〈x1(t),c21v1(t)+a2(t)〉=c21〈x1(t),v1(t)〉+〈x1(t),a2(t)〉=c21+〈x1(t),a2(t)〉=0
implying
(28)c21=-〈x1(t),a2(t)〉andhencev2(t)=-〈x1(t),a2(t)〉v1(t)+a2(t)
with
(29)x2(t)=〈KX(t,s),v2(s)〉〈〈KX(t,s),v2(t)〉,v2(s)〉,(iii)
assumingv1(t),v2(t),…,vk-1(t) and x1(t),x2(t),…,xk-1(t) known, define vk(t)=ck1v1(t)+ck2v2(t)+⋯+ck,k-1vk-1(t)+ak(t) with ck1,ck2,…,ck,k-1 determined from
(30)〈x1(t),vk(t)〉=〈x2(t),vk(t)〉=⋯=〈xk-1(t),vk(t)〉=0
with
(31)xk(t)=〈KX(t,s),vk(s)〉〈〈KX(t,s),vk(t)〉,vk(s)〉   (see, e.g., [1] for more details and applications).Notice that when the given sequence{ak(t)} of square-integrable functions are orthogonal, then the generated biorthogonal sequence {vk(t)} is simply the original orthogonal sequence {ak(t)}. Furthermore, if the given sequence {ak(t)} happened to be the eigenfunction sequence {uk(t)} from KX(t,s), then the generated biorthogonal sequence {vk(t)} would simply be the eigensequence {uk(t)} with normalizing coefficients being the corresponding eigenvalues λk and the situation reducing to the previous orthogonal case with unit weight function w(t)≡1.For a random functionX(t) over some domain T with zero-mean value and covariance kernel KX(t,s), let {vk(t),t∈T} denote a sequence of square-integrable functions such that {xk(t)} and {vk(t)} are biorthogonal function sequences on T. Then the following conditions are equivalent:(i)
X(t) is equal in the mean-square sense to its canonical expansion in terms of ξk and xk(t); that is,
(32)X(t)=∑k=1∞〈X(t),vk(t)〉xk(t)=∑k=1∞ξkxk(t);(ii)
KX(t,s) has a canonical expansion in terms of κk, xk(t), and xk(s)(33)KX(t,s)=∑k=1∞κkxk(t)xk*(s);(iii)
Var[X(t)] has a canonical expansion in terms of ξk and xk(t)(34)Var[X(t)]=∑k=1∞ξk|xk(t)|2.More discussion and applications of these results can be found in, for example, [1].Again in some applications where an arbitrary orthogonal function sequence is available such as from previous measurements or observations, the biorthogonal basis function approach is generally advantageous over other approaches that do not make use of the prior information.


## 5. Nonstationary and Nonlinear Stochastic Analysis

Although most nonstationary stochastic processes are not amenable to linear analysis, there are exceptional cases where a simple variable transformation renders linear analysis possible.  For instance, multiplicative noise contamination becomes additive “noise” under a logarithm transformation, and hence linear analysis can be done in the logarithmic space.Some stochastic processes become stationary under differentiation in continuous contexts or differencing in discrete contexts. For example, a (covariance stationary) process with an unknown linear trend as first moment becomes stationary under double differentiation or differencing as commonly done in Global Positioning System (GPS) data processing. Such processes are called increment stationary and their second moments remain stationary under linear operations of differentiation or differencing. Obviously only finitely many such operations are considered corresponding to the (finite) degree of the polynomial trend.Assuming a stochastic processX(t) to have stationary increments of some order k≥0, then ∇k+1X(t)≡X(k+1)(t) is stationary with some generalized covariance G(P,Q)≡G(P-Q). That means that the generalized Fourier transform ℱ[G(P,Q)]>0 for realizability, but using the well-known Khintchine-Wiener theorem, the spectrum of X(k+1)(t) is also given by the generalized Fourier transform of the autocorrelations of X(k+1)(t), that is, |ℱ[X(k+1)(t)]|2. Using the mathematical relation ℱ[X(k+1)(t)]=(iω)k+1ℱ[X(t)] for frequencies ω, then
(35)ℱ[G(P,Q)]=|ℱ[X(k+1)(t)]|2≡ℱ[X(k+1)(t)]·ℱ[X(k+1)(t)]*=ω2k+2|ℱ[X(t)]|2>0
which implies the existence and positive definiteness (for realizability) of a generalized covariance G(P,Q)≡G(P-Q) for X(t) whenever one exists for X(k+1)(t).Assuming thatX(t) is increment stationary, then its expected value can be modeled using an algebraic polynomial MX(t)=a0+a1t+a2t2+⋯+aktk with unknown coefficients a0,a1,a2,…,ak, and the previous Wiener-Hopf equations for unknown H(r,t)(36)〈Kyy(s,r),H(r,t)〉=KyX(s,t)
have to be generalized to the so-called extended normal equations,
(37)〈Kyy(s,r),H(r,t)〉+αMX(t)=KyX(s,t),〈MX(r),H(r,s)〉=MX(s),
for some Lagrange multiplier α for constrained optimization.Explicitly, using the terminology and conventions of regionalized variables in geostatistics, the linear prediction forf(x0)=∑i=1Nλif(xi), with an algebraic polynomial trend or drift model a0+a1x+⋯+akxk for the function f(x), and denoting Cij=E[f(xi)f*(xj)], the normal equations become
(38)[C11⋯C1N1⋯x1k⋮⋯⋮⋮⋯⋮CN1⋯CNN1⋯xNk1⋯10⋯0⋮⋯⋮⋮⋯⋮x1k⋯xNk0⋯0][λ1⋮λNαa0⋮αak]=[C10⋮CN01⋮x0k].
This is the same extended form of normal equations encountered with empirical Radial Basis Functions (RBFs) instead of (generalized) covariance functions. Notice that when the mean is known and constant, the process can be centered with the Lagrange multiplier α set to zero, and these extended normal equations reduce to the standard Wiener-Hopf equations.Such an extended covariance matrix is easily invertible under minimal conditions as(39)(CBBT0)-1=(C-1+C-1BHBTC-1-C-1BH-HBTC-1H)withH=-(BTC-1B)-1,
where H is sometimes called Schur’s Complement of C in the extended matrix (see, e.g., the work of Kailath et al. [3] for more details). Hence under minimal conditions, these extended normal equations can be solved explicitly for the linear coefficients λν and the trend polynomial coefficients. The inverse of a covariance matrix is usually called an information matrix. Also, generalized covariance functions in the sense of generalized functions are often advantageous for estimation in practical applications (see, e.g., the work of Blais [11] for more discussion with simulations).Recalling that given observations{y(s)}, the optimal general estimate (not necessarily linear)
(40)X^(t)=E[X(t)∣y(s)]
with the property that
(41)Var[X(t)-E[X(t)∣y(s)]]≤E[|X(t)-Ψ(y(s))|2]
for any linear or nonlinear estimator Ψ(y(s)) can sometimes be evaluated directly using numerical simulations. In fact, using the previous conventions,
(42)E[X(t)]≡μX(t)≡∫-∞∞xpX(t)dx,
and with the help of Bayes’ theorem,
(43)X^(t)=∫-∞∞xpX∣y(t,s)dx=∫-∞∞x{py∣X(t,s)pX(t)∥py∣X(t,s)pX(t)∥}dx
which can be evaluated directly using Monte Carlo simulations in simple cases.For example with Gaussian variablesZν such that E[Zν∣u]=σν2∑ανηuη and variance σν2, one can write explicitly for the conditional probability density function
(44)p(Zν∣u)=12πσν2exp{-(Zν-σν2∑ανηuη)22σν2},
and these can be multiplied in cases of several independent Gaussian variables. In general, such explicit formulation is obviously unknown. Using Bayes’ Theorem as above and substituting the results into the integral for X^(t), the required conditional expectation can be evaluated in simple applications. More discussion and application examples can be found in the works of Pugachev [2], Dougherty [1], and Zhang [12].


## 6. Computational Examples and Discussion

In many practical applications, solving the preceding normal equations and/or their extended form in the most efficient manner is one primary objective. As mentioned above, in stationary cases with equispaced data, the use of FFTs is the common strategy as the diagonalizing of the covariance matrix is well known to be most efficient. However, this is not always possible in many situations. Computational alternatives are briefly discussed in the following paragraphs.As different software packages use different conventions for discrete Fourier, wavelet Haar, and other transforms, the following is very explicit in the conventions used. First, the Fourier matrixF of order N has elements frs=e2πirs/N, and the simple application with N=8 for the covariance function Cov[d]=1/(1+d2) gives (without any normalization)
(45)F*(112151101171101512)=(2.6591.5070.6590.3750.2590.3750.6591.507),
or using the corresponding circulant Toeplitz matrix (with proper normalization),
(46)F*(112151101171101512121121511011711015151211215110117110110151211215110117117110151211215110110117110151211215151101171101512112121511011711015121)F=8·diag(2.6591.507⋯0.6591.507)
as expected. Notice that formulating the preceding covariance matrix differently for the given covariance function Cov[d]=1/(1+d2) would not imply the above diagonal matrix. Also, the scaling of this diagonal matrix would not be needed with normalized Fourier transforms.Again considering the synthetic covariance function(47)Cov[d]=1α2+d2
for some real scalar α, the plot of Cov[d] for d=0,1,…,512 (using α=1 for simplicity) with corresponding power spectrum or Fourier transform of this covariance function is simply
(48)S(ω)=1αe-α|ω|
for frequencies ω which is shown in Figure 1(a). Adding some noise amplification within -1≤s≤1 in space or time, the covariance function becomes
(49)Cov[s,d]=1α2+d2+H(|s|)δ(d)
for the Heavyside (or step) function
(50)H(|s|)={1for|s|≤1,0otherwise
and scalar α; its spectrum becomes
(51)S(σ,ω)=1αe-α|ω|+σ-1sinσ
which is nondiagonal as shown in Figure 1(b) with ripples away from the diagonal.(a) Exponential spectrum.  (b) Exponential + Sinc spectrum. (c) Haar spectrogram of[1/(1+d2)]. (d)  Haar spectogram of [1/(1+d2) + Nonstationary Noise]. (e) Singular spectrum of extended covariance matrix. (f)  Power spectrum of extended covariance matrix.
(a)
(b)
(c)
(d)
(e)
(f)Considering this covariance functionCov[d] for multiresolution analysis with the Haar wavelet basis, using the discrete orthogonal basis (1/21/2) and (1/2-1/2) with order 512, the corresponding spectrogram is shown in Figure 1(c) which is sparse but not diagonal as the power spectrum. Modifying the level of noise with this covariance function as follows
(52)Cov[s,d]=1α2+d2+{0.5fors=1,2,…,256,0.2fors=257,258,…,512
does not imply much variation in the corresponding spectrogram (see Figure 1(d)). Similar results were obtained by Keller [13] with Haar wavelets and Zhang [12] using a slightly different orthogonal basis. Although the orthogonal Haar wavelet and similar transforms would seem most appropriate for such covariance matrices, they only succeed in making the covariance matrix sparse but not diagonal. The EOF approach does diagonalize any symmetric positive definite covariance matrix by the very nature of the eigenvector decomposition. Such eigenvector and eigenvalue decompositions are computationally demanding, but it is worth noting that the biorthogonal basis construction is computationally advantageous whenever some arbitrary basis is available, as mentioned above.Next, assuming an unknown mean for the simulated stochastic process, the covariance matrixCov[d] can easily be modified such that the last column and row are 1’s with 0 as the last diagonal element; the corresponding singular spectral matrix is shown in Figure 1(e). Using this modified covariance matrix Cov[d], its DFT spectral matrix is shown in Figure 1(f) which shows anomalous effects due to the border modifications of Cov[d]. In other words, the DFT and FFT approaches are not appropriate with extended covariance matrices even with equispaced data. This is obviously an important result for applications of extended covariance and RBF normal matrices in geostatistics (see, e.g., the works of Chilès and Delfiner [14]).


## 7. Concluding Remarks

Stochastic time series are very common in geoscience, and better analytical tools and procedures for their modeling and filtering warrant more research and development. Multiresolution analysis and synthesis of sequences of observations and measurements are needed for all kinds of applications in geoscience and other fields.DFTs and FFTs are normally used in stationary processes with equispaced data. However in practice, nonstationarities are often unavoidable, and even in increment stationary cases with equispaced data, DFTs and FFTs are not appropriate. Galerkin’s projection methods especially with EOFs are more general and unfortunately more computationally demanding. However in some applications, only the leading spectral modes are necessary for analysis and inference thus reducing the necessary computational efforts.For simple prediction applications in stationary cases, the normal equations can be solved using Cholesky’s Square-Root algorithm which has well-known numerical properties; while in increment stationary cases, as well as various RBF formulations, some Cholesky and Givens reduction procedure can be used (see, e.g., the work of Blais [10]). For more in-depth numerical analysis, the SVD and KLT are most recommendable in general and visualization tools are very helpful in the structural analysis of covariance and information matrices.



---

*Source: 895061-2013-05-16.xml*