Numerical Treatment of Covariance Stationary Processes in Least Squares Collocation

Abstract

Digital sensors provide long series of equispaced and often strongly correlated measurements. A rigorous treatment of this huge set of correlated measurements in a collocation approach is a big challenge. Standard procedures – applied in a thoughtless brute force approach – fail because these techniques are not suitable to handle such huge systems.

In this article two different strategies, denoted as covariance approach and filter approach, to handle such huge systems are contrasted. In the covariance approach various decorrelation strategies based on different Cholesky approaches to factorize the variance/covariance matrices are reviewed. The focus is on arbitrary distributed data sets with a finite number of data. But also extensions to sparse systems resulting from finite covariance functions and on exploiting the Toeplitz structure which results in the case of equispaced systems are elaborated.

Apart from that filter approaches are discussed to perform a prewhitening strategy for the data and rearrange the whole model to work with this filtered data in a rigorous way. Here, the special focus is on autoregressive processes to model the correlations. Finite, causal, non-recursive filters are constructed as prewhitening filters for the data as well as the model. This approach is extreme efficient, but can only deal with infinite equispaced data sets.

In real data scenarios, finite sequences and data gaps must be handled as well. For the covariance approach this is straightforward but it is a serious problem for the filter approach. Therefore a combination of these approaches is constructed to select the best properties from each. Covariance matrices of equispaced data sets designed by recursively defined covariance sequences are represented by AR processes as well as by Cholesky factorized matrices. It is shown, that it is possible to switch between both strategies to get data gaps and the warm up phase for the filter approach under control.

Zusammenfassung

Digitale Sensoren liefern Zeitreihen von gleichabständigen und oft stark korrelierten Messungen. Eine strenge Auswertung, dieser zumeist umfangreichen Datensätze, in einem Kollokationsansatz stellt eine große Herausforderung dar. Standardverfahren sind nicht in der Lage solch große Systeme zu bewältigen.

In diesem Artikel werden zwei unterschiedliche Verfahren – der Kovarianzansatz und der Filteransatz – einander gegenübergestellt, um deren Potential für den Einsatz bei sehr großen Datenmengen zu diskutieren. Im Kovarianzansatz erfolgt die Modellierung der Korrelationen durch Kovarianzfunktionen. Die Dekorrelation der Messungen erfolgt über die Faktorisierung der Kovarianzmatrizen, wofür unterschiedliche Varianten der Cholesky Faktorisierung untersucht werden. Dünn besetzte Systeme hervorgerufen durch den Einsatz von finiten Kovarianzfunktionen, regelmäßige Toeplitz Systeme resultierend aus der regelmäßigen Abtastung, aber auch die Auswirkungen von lokalen Datenverlusten (Datenlöcher) werden untersucht und maßgeschneiderte Zerlegungsverfahren entwickelt.

Der Filteransatz bietet hingegen kaum Möglichkeiten flexibel auf Eigenheiten der Daten einzugehen. Für regelmäßige Abtastung von sehr langen Zeitreihen, ist die Dekorrelation über Filter extrem effizient. Durch die Modellierung der Kor- relationen durch finite, autoregressive Prozesse (AR-Prozesse) werden kausale, nichtrekursive, finite Filter aufgebaut, die unendlich ausgedehnte Meßreihen effizient dekorrelieren können. Datenlöcher, aber auch die Initialisierung des Filterprozesses bewirken aber einen Datenverlust, der speziell bei großen Korrelationslängen durchaus dramatisch sein kann.

Messserien aus der Praxis bestehen aber aus endlich vielen Messungen, enthalten Unregelmäßigkeiten und Datenlöcher. Während der Kovarianzansatz ohne Probleme auf diese Daten anwendbar ist, sind beim effizienten Filteransatz spezielle Vorkehrungen (Näherungen) notwendig. Durch den Einsatz von rekursiv definierten Kovarianzsequenzen kann jedoch ein Übergang zwischen den beiden Ansätzen hergestellt werden. Die vollbesetzte Kovarianzmatrix kann mit Hilfe einer speziellen Variante der Cholesky Inversion – rekursives rückwärts Rändern – in eine dünnbesetzte Matrix zerlegt werden, die weitgehend eine Band-Toeplitz-Struktur aufweist, die dem kausalen, nichtrekursiven Filter im Filteransatz entspricht. Durch die Kombination der beiden Ansätze ist es somit möglich einen zeitvariablen, kausalen, nicht rekursiven Filter zu entwickeln, der endliche und leicht unregelmäßige Datensätze streng und effizient dekorrelieren kann.

This chapter is part of the series Handbuch der Geodäsie, volume “Mathematical Geodesy/ Mathematische Geodäsieʺ, edited by Willi Freeden, Kaiserslautern.

Appendix

1.1 A Cholesky Approach: Revisited

The Cholesky approach was first published posthumously by [3] in Bulletin Géodésique and based on an unpublished manuscript from [10]. The solution of symmetric linear equation systems is a frequently asked problem in numerical mathematics and data analysis. The efficient and stable solution of the normal equations as well as the factorization and inversion of covariance matrices is an inherent task in least squares adjustment and least squares collocation. Thus, many geodesists use the Cholesky approach and developed a lot of strategies to adopt this method for their special applications [28, 32]. Sparse solvers [37, 45] as well as partial inversion strategies [2, 17] are developed, studies on the numerical stability and the influence of rounding errors in huge systems are performed [12, 27], but also strategies to generate correlated signals for Monte Carlo approaches [1] are developed. Usually Cholesky factorization is connected to symmetric positive definite systems, but in Fig. 7 also other possible applications for infinite, rectangular, asymmetric and ill-posed systems are mentioned. These systems are out of scope of this article, but in many situations the Cholesky approach is extremely efficient and numerically stable applicable.

Fig. 7

Application of the Cholesky approach on different types of systems. This figure is focused on the application on symmetric systems and shows the variety of different strategies only for this special type

1.1.1 A.1 Cholesky Solution

The Cholesky approach was originally developed for the solution of a symmetric linear equation system by the decomposition into a lower and upper triangular matrix. In contrast to the LU decomposition the diagonal elements of both matrices are chosen equal to eliminate the degree of freedom. For symmetric systems the hermitian matrix is decomposed into triangular matrices, where one is the conjugate transposed of the other. In this work we choose the notation

(54)

where R denotes a upper (right) triangular matrix and its conjugate transposed (hermitian) form. If is a positive definite, real valued, symmetric matrix the decomposition is real valued, and thus

(55)

with . This approach denoted a Cholesky forward factorization is often used for least squares problems, especially for Gauss-Markov models and adjustment models with condition equations. The resulting positive definite normal equations can be decomposed very efficiently (roughly twice as efficient as the LU decomposition) and numerically stable [12, 27] by the Cholesky approach.

The solution of the linear symmetric system

(56)

with the symmetric matrix and the right hand side can be efficiently determined after the factorization

(57)

in a two step approach by solving the lower triangular system

(58)

for the auxiliary unknowns in a forward substitution step and then solving the upper triangular system

(59)

for the unknowns x in a backward substitution step. The formulas for each single coefficient of the forward and backward substitution can be derived immediately from (58)

$$\displaystyle \begin{aligned} z_i = \left( n_i - \displaystyle\sum_{k=1}^{i-1} r_{ki} z_k\right) / r_{ii}\ ,\qquad i=1,\ldots,n, \end{aligned} $$

(60)

and (59)

$$\displaystyle \begin{aligned} x_i =\left( z_i - \displaystyle\sum_{k=i+1}^n r_{ik} x_k\right) / r_{ii}\ ,\qquad i=n,\ldots,1 \quad . \end{aligned} $$

(61)

The coefficients of the triangular matrix R can be derived from the identity (55). Written in expanded form

$$\displaystyle \begin{aligned} \begin{bmatrix} n_{11} & n_{12} & \ldots & n_{1n}\\ n_{12} & n_{22} & \ldots & n_{2n}\\ \vdots & \vdots & \ddots & \vdots\\ n_{13} & n_{23} & \ldots &n_{nn} \end{bmatrix} = \begin{bmatrix} r_{11} & & \\ r_{12} & r_{22} & \\ \vdots & \vdots & \ddots\\ r_{13} & r_{23} & \ldots& r_{nn} \end{bmatrix} \begin{bmatrix} r_{11} & r_{12} &\ldots& r_{1n} \\ & r_{22} &\ldots& r_{2n} \\ & & \ddots& \vdots\\ & & & r_{nn} \end{bmatrix}, \end{aligned} $$

(62)

it can be seen immediately that n² – respectively $\frac 12n(n+1)$ linear independent – identities results from this equation. They can be used to derive the Cholesky coefficients by

(63)

It should be mentioned, that the three required loops with respect to i, j and k can be organized in accordance with the used storage scheme of the matrices. These numerical considerations are not in the focus of this article (see for instance [40, p. 138ff]).

Figure 8 summarizes the four steps which are required to determine a solution of a linear symmetric system with the Cholesky approach. It should be mentioned, that the factorization step for off-diagonal elements exactly corresponds to the forward substitution step. The first three steps use a column access only, whereas for the backward substitution a row and a column access is proposed.

Fig. 8

Graphical interpretation of the factorization steps (63), the forward substitution step (60) and the backward substitution step (61). (a) Factorization: diagonal element. (b) Factorization: off-diagonal element. (c) Substitution: forward. (d) Substitution: backward

These figures provide a feeling, what happens during the factorization of sparse input matrices. If an element n_ij is initially zero, it is unchanged during factorization, if the scalar product of the column vectors r_ki and r_kj vanishes. This is the case especially for envelope systems, where the column above a zero element is zero as well. Therefore the envelope structure of a sparse symmetric system is preserved during Cholesky factorization [15, 28]. Special numbering schemes are developed [8, 37, 40, 45] to determine the order of the parameters in the matrix, to optimize the structure to preserve sparsity also during the factorization step. The design of covariance matrices by finite covariance functions [13, 21, 34, 42] yields to sparse systems, which can efficiently solved by the Cholesky approach.

Figure 8 can also be used to get an impression of the case of an indefinite system. In the factorization step of the diagonal element, the value under the square-root becomes negative. Therefore it is necessary to introduce the imaginary unit i. In case this happens, the divisor in all off-diagonal factorization steps is imaginary, therefore the entire row consists of imaginary values. The main operation during the whole factorization and substitution steps are scalar products between two columns. If a row contains imaginary values in both column elements, the result is not complex anymore, but real valued. The product of these two elements is not added during the computation of the scalar product, but must be subtracted. Therefore, indefinite systems can be solved easily by the Cholesky approach by introducing an index vector to indicate which rows involve imaginary numbers. As the normal equations of the commonly used Gauss-Markov model, or the covariance matrices are positive definite systems, this extension of the Cholesky approach is rarely used in geodesy. But the normal equations of the Gauss-Helmert model and also the Gauss-Markov model with restrictions are indefinite systems and the above mentioned extension is crucial to take benefit from the stable and efficient Cholesky approach also for indefinite systems.

1.1.2 A.2 Cholesky Inversion

In this section we are looking for the inverse of the matrix N, where we assume, that the Cholesky factorized matrices R are already computed. A common used strategy is to invert the triangular matrix R and apply

(64)

Here we want to elaborate a direct strategy to switch from to the inverse matrix .

The positive definite matrix N and its inverse can be written in block notation as

(65)

where denotes the upper left block in the inverse . It has to be carefully distinguished from the inverse of the upper left block of the matrix N, which is denoted by . The partitioning of the matrix into blocks is performed in a way, that and are quadratic. The identity

(66)

can be used to find a representation of the blocks of the inverse matrix. They read

(67)

(68)

(69)

Performing a Cholesky factorization of the block partitioned matrix N

(70)

the relation between the matrices and the matrices are given by

(71)

Substituting these quantities in (67–69) the well-known relations between the blocks of the inverse matrix and the Cholesky factorized block matrices can be derived as

(72)

(73)

(74)

To avoid the inversion of large systems the computations are performed by a backward substitution row-by-row starting with the element nnn(−1)

(75)

Figure 9 gives a graphical interpretation of (75). It shows the necessary operations to invert the row i. At this step it is assumed, that for the rows i + 1, …, n the inverse elements are already computed and for the rows 1, … i the Cholesky factorized form is available. By this backward substitution process the whole inverse can be computed.

Fig. 9

Graphical interpretation of the inversion step (75) starting from a Cholesky factorized matrix. By a backward process the inverse matrix is computed row-by-row. In contrast to the upper triangular Cholesky matrix R the inverse matrix is a symmetric matrix. The skeleton-like squares marks this symmetric part, which is typically not stored

The computation of the inverse matrix, with the Cholesky factorization available, can be performed in the order of operations, which is twice expensive as the Cholesky factorization step itself, where only operations are necessary. All in all the inversion (Cholesky and recursive backward edging) process is in the order of . It is well-known, that in general the inversion destroys sparsity, such that the inverse is typically a dense matrix.

In contrast to many other engineering and research areas, where only the solution of the system is requested, in statistics, adjustment theory and geodesy the inverse of the system is of special interest, because the uncertainties and the covariances are deduced from the inverse. The knowledge of this information is crucial in most applications. Sometimes, we are only interested in special elements of the inverse, e.g., the main diagonal elements or elements in a defined bandwidth or off-diagonal elements between special nodes which are characterized by their neighborhood or by points directly connected by measurements. These are typically the entries, which are already non zero in the sparse normal equations. Error propagation within geodetic networks [37], Helmert-Blocking [47] or nested dissection [14] partitioning strategies to solve large finite element mesh systems and kite-structures [39, 40] to perform a spherical harmonic analysis with irregularly fully populated low degree and block-diagonal structured high degree fields are typical applications.

In connection with the new adjustment of the North American horizontal datum [5] a special algorithm to compute the inverse elements for an envelope structured sparse matrices was introduced and denoted as partial inverse. It is shown by Hanson [17], that all the inverse elements within the envelope can be computed in a strict manner without any information about the inverse elements outside of the envelope. This is due to the fact, that the non-computed elements within each column of the inverse corresponds exactly to the zero-elements within the Cholesky reduced row. Therefore, the scalar product in (75) (cf. also the computation of in Fig. 9) for all elements within the envelope is not influenced by inverse elements, that are positioned outside the envelope. This argument holds also for special structures coming from a nested dissection partitioning or for kite-structures [2].

1.1.3 A.3 Cholesky Factorization of a Matrix Starting from Its Inverse

In many statistical applications the inverse of the matrix (covariance matrix) plays a central role. For some special applications, it can be useful to compute the Cholesky factors of the initial matrix N directly from the inverse matrix (covariance matrix). To obtain an algorithm for that purpose, it is necessary to reorder the relations between the Cholesky coefficients and the elements of the inverse given in (75) with respect to the Cholesky coefficients. For the last diagonal element from the inverse nnn(−1), we get immediately the last Cholesky coefficient by

$$\displaystyle \begin{aligned} r_{nn} = 1 / \sqrt{{n^{(-1)}}_{nn}} \ . \end{aligned} $$

(76)

Now in a step-wise approach the Cholesky coefficients can be computed row-by-row by rearranging (75) with respect to the Cholesky coefficients. This yields to

$$\displaystyle \begin{aligned} {n^{(-1)}}_{ij} r_{ii} + {n^{(-1)}}_{i+1,j} r_{i,i+1} + \ldots + {n^{(-1)}}_{n,j} r_{in} = 0 \ ,\quad \quad j=i+1, \ldots,n\end{aligned} $$

(77)

and respectively

$$\displaystyle \begin{aligned} {n^{(-1)}}_{ii} r_{ii} + {n^{(-1)}}_{i+1,i} r_{i,i+1} + \ldots + {n^{(-1)}}_{n,i} r_{in} = \frac 1{r_{ii}} \ . \end{aligned} $$

(78)

After multiplying all equations by r_ii and a rearranging, we end up with a system of linear equations

$$\displaystyle \begin{aligned} \begin{bmatrix} r_{ii} r_{ii}\\ r_{ii} r_{i,i+1}\\ \vdots\\ r_{ii} r_{in}\\ \end{bmatrix} = \begin{bmatrix} {n^{(-1)}}_{ii} & {n^{(-1)}}_{i,i+1} & \ldots & {n^{(-1)}}_{in} \\ {n^{(-1)}}_{i,i+1} & {n^{(-1)}}_{i+1,i+1} & \ldots & {n^{(-1)}}_{i+1,n} \\ \vdots & \vdots & \ddots & \vdots\\ {n^{(-1)}}_{in} & {n^{(-1)}}_{i-1,n} & \ldots & {n^{(-1)}}_{nn} \end{bmatrix}^{-1} \begin{bmatrix} 1 \\ 0 \\ \vdots\\ 0 \end{bmatrix} \qquad \text{for} \qquad i=n-1,\ldots,1 \ \ .\end{aligned} $$

(79)

The Cholesky coefficients of row i multiplied by r_ii are therefore equivalent to the first row of the corresponding sub-matrix of the inverse. The square root of the first element provides r_ii and all the other elements can be computed by dividing the elements by r_ii. On a first glance this factorization process seems to be very inefficient, because in each step the computation of the inverse of the sub-block is necessary. But this successive inversion can be done by recursive backward edging

(80)

assuming that the inverse respectively are already computed and the inverse of the by one row and column extended matrix should be determined [41, p.13f]. The required operations of the overall procedure is in the order of . In applications from signal processing, often deal with covariance matrices resulting from equispaced data. These matrices are Toeplitz structured and the Levinson-Durbin recursion allows for a very efficient solution with operations ( [25]; [11]; [16, p.128–129]; [41, p.291ff]).

1.1.4 A.4 Backward Cholesky Factorization

In the standard Cholesky factorization procedure a symmetric matrix N is split into a lower triangular matrix and an upper triangular matrix, where one is the transposed of the other, , where R denotes the upper (right) triangular matrix. In the backward Cholesky factorization the symmetric matrix is factorized into an upper triangular matrix and a lower triangular matrix. This yields to

(81)

or in extended notation to

(82)

This factorization can be seen as the standard Cholesky approach for a top-down and left-right flipped symmetric system. The factorization starts at the last diagonal element n_nn and processes row by row to the top. The notation with the bar on the top should express, that in general the coefficients r_ij differ from . From the identity (82), the recursion formulas for the computation of the coefficients R can be immediately deduced

(83)

In correspondence to Sect. A.2, the recursion formulas for the coefficients of the inverse of the symmetric matrix can be computed by forward substitution,

(84)

and can be reorganized with respect to the computation of the backward Cholesky factors from the inverse

(85)

(86)

The computation of the backward Cholesky factors of the column j can be organized as the solution of the linear system of equations

(87)

as a recursive forward edging process. If the inverse matrix is Toeplitz structured the Levinson-Durbin algorithm [11, 25, 46] can be applied.

1.1.5 A.5 Resumé on Cholesky Factorization

In the last sections two different approaches of the Cholesky factorization are discussed. The standard or forward Cholesky factorization (63), where the symmetric matrix N is split into a lower triangular matrix and an upper one

(88)

and the backward Cholesky factorization (83), where the symmetric matrix N is factorized into

(89)

It is important to note that in general R is different from R . In a further step, the inverse matrix can be computed directly from the Cholesky factorized triangular matrices R or R by backward substitution (75) or forward substitution (84). It is also presented, how the Cholesky factorized matrices R and R of the initial matrix N can be derived from the given inverse matrix . The matrix R can be computed by recursive backward edging (79) and the matrix R by recursive forward edging (87). The circle can be closed by a simple matrix multiplication (88) or (89) to come back to the initial matrix N from the triangular matrices R and R . Fig. 10a, b summarizes this possible paths.

Fig. 10

Review the Cholesky factorization steps applied to N and . (a) (Forward) Cholesky factorization of N. (b) Backward Cholesky factorization of N. (c) Cholesky approach of . (d) Backward Cholesky approach of

A further interesting perspective is to start with the inverse matrix instead of the matrix N and set up the same sequence of Cholesky factorization steps as done before. Formally the inversion of the matrix N with respect to the forward Cholesky factorization (88) can be written

(90)

and

(91)

with respect to backward Cholesky factorization (89). It is now crucial that in contrast to the common factorized matrices their inverse forms and occur. Changing the view, also the inverse matrix can be factorized into a lower and upper triangular matrix by the forward Cholesky factorization (55) we see, that we get immediately or by the backward Cholesky factorization (83) respectively. This means, starting from the inverses of the triangular matrices can be directly computed. In the next steps the Cholesky inversion by backward substitution (75) or forward substitution (84) can be computed to get the inverse of the inverse which is for sure the matrix N again. But also the other way is now possible, i.e., the computation of the matrix directly from the matrix N by recursive backward edging (79) and by recursive forward edging (87) the matrix . In Fig. 10 all these factorization steps are summarized for direct factorization of the matrix N as well as for .

The Cholesky approach allows a lot of variations to factorize a symmetric matrix and opens a way to directly compute R and R or or . It should be mentioned that some of the factorization steps can be extreme efficiently computed for Toeplitz structures resulting from equispaced data distribution.