Parallel Computing in Spatial StatisticsThis means if we translate a point conﬁguration in a given...

Nadja Samonig

Parallel Computing in SpatialStatistics

Diplomarbeit

zur Erlangung des akademischen Grades

Diplomingenieurin

Studienrichtung: Technische Mathematik

Studienzweig: Angewandte Wirtschaftsmathematik

Universitat Klagenfurt

Fakultat fur Wirtschaftswissenschaften und Informatik

Begutachter: Dipl.-Math. Dr. Albrecht Gebhardt

Institut fur Mathematik

Juli/2001

EHRENWORTLICHE ERKLARUNG

Ich erklare ehrenwortlich, daß ich die vorliegende Schrift verfaßt und alle ihr voraus-

gehenden oder sie begleitenden Arbeiten durchgefuhrt habe. Die in der Schrift verwen-

dete Literatur sowie das Ausmaß der mir im gesamten Arbeitsvorgang gewahrten Un-

terstutzung sind ausnahmslos angegeben. Die Schrift ist noch keiner anderen Prufungsbehorde

vorgelegt worden.

Riegersdorf, am 9. Juli 2001

ΠPOMHΘEYΣ

και µην ε%γω κoυκετι µυθω

χθων σεσαλευται.

β%υχια δ′ ηχω πα%αµυκαται

β%oντ ης, ελικες δ′ εκλαµπoυσι

στε%oπης ζαπυ%oι, στ%oµβoι δε κoνιν

ειλισσoυσιν σκι%τα δ′ ανεµων

πνευµατα παντων εις αλληλα

στασιν αντ ιπνoυν απoδεικνυµενα

ξυντετα%ακται δ′ αιθη% πoντω.

τoιαδ′ επ′ εµoι %ιπη ∆ιoθεν

τενχoυσα ϕoβoν στειχει ϕανε%ως.

ω µητ%oς εµης σεβας, ω παντων

αιθη% κoινoν ϕαoς ειλισσων,

εσo%ας ως εκδικα πασχω.

AIΣXYΛOΣ, ΠPOMHΘEYΣ ∆EΣMΩTHΣ

Contents

1 Introduction 4

2 Spatial Statistics 6

2.1 Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

2.1.1 Spatial Data. . . . . . . . . . . . . . . . . . . . . . . . . . . 7

2.1.2 Regionalized Variable. . . . . . . . . . . . . . . . . . . . . 7

2.2 Covariance Structure. . . . . . . . . . . . . . . . . . . . . . . . . . 8

2.2.1 Stationarity. . . . . . . . . . . . . . . . . . . . . . . . . . . 8

2.2.2 Simple Kriging. . . . . . . . . . . . . . . . . . . . . . . . . 9

2.2.3 Semivariogram and Covariance Function. . . . . . . . . . . 10

2.2.4 Semivariogram Models. . . . . . . . . . . . . . . . . . . . . 12

2.2.5 Anisotropy . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

2.3 Kriging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

2.3.1 Estimation of the Mean. . . . . . . . . . . . . . . . . . . . . 18

2.3.2 Ordinary Kriging. . . . . . . . . . . . . . . . . . . . . . . . 19

2.3.3 Universal Kriging. . . . . . . . . . . . . . . . . . . . . . . . 22

1

CONTENTS 2

3 Kriging on Tiled Grids 26

3.1 Solving Kriging Systems Simultaneously. . . . . . . . . . . . . . . 26

3.2 Tiled Grid Kriging . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

3.3 Tiled Grid Kriging with Overlapping. . . . . . . . . . . . . . . . . . 33

3.4 Example. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

4 PVM (Parallel Virtual Machine) 43

4.1 Introduction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

4.2 PVM’s basics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44

4.3 The PVM System. . . . . . . . . . . . . . . . . . . . . . . . . . . . 45

5 Universal PVM Kriging 47

5.1 Combining R and PVM. . . . . . . . . . . . . . . . . . . . . . . . . 47

5.2 Parallel Tile Kriging . . . . . . . . . . . . . . . . . . . . . . . . . . 50

5.3 Example. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54

6 Conclusion and Prospect 59

A Functions 60

A.1 Kriging Functions. . . . . . . . . . . . . . . . . . . . . . . . . . . . 61

A.1.1 The Functionkrige.grid . . . . . . . . . . . . . . . . . . 61

A.1.2 The Functionkrige.tiles . . . . . . . . . . . . . . . . . 64

A.1.3 The Functionkrige.tilesov . . . . . . . . . . . . . . . 68

A.2 PVM Kriging Function . . . . . . . . . . . . . . . . . . . . . . . . . 72

CONTENTS 3

A.2.1 The Functionkrige.pvm . . . . . . . . . . . . . . . . . . 72

Bibliography 78

Index 79

Chapter 1

Introduction

The object of my thesis is the parallelization of a mathematical computation problem.

Mathematically largescale and wearisome computations can be handled by parallel

processing. This method solves one large problem by dividing it into many small tasks.

My first task was to modify a FORTRAN 77 subroutine for kriging prediction (see sub-

section2.3.3) which was prepared by DR. ALBRECHT GEBHARDT. This subroutine

krgtil.f devides one big krige system into small ones in order to boost computa-

tional speed. This one was enlarged and improved in order to enhance computational

exactness. The new subroutinekrgtilov.f has the modification that the small sys-

tems, the tiles, can overlap. Now some tile points get computed more than once and the

effect is that the prediction is more exact (see section3.4), especially at the tile borders.

The problem of the computation of tile points is very suited for parallelisation. Finally

we yield a big saving of time by distributing small computational tasks to different

computers. In the parallelized version the FORTRAN 77 subroutinekrgtilov.f

gets replaced with the R functionkrige.pvm.R which envokes calls to the C rou-

tine krige pvm.c , in order to send the data prepared in R to a server. The task of

4

CHAPTER 1. INTRODUCTION 5

the server is to send the tiles of data points, the jobs, to the clients which perform the

computation. In order to distribute the data from the server to the clients, PVM (see

chapter 4) is used. When a client is ready with the computation of the job, the results

get sent back to the server. Subsequently the server furnishes the client with further

jobs, as long as non-computed jobs exist. Finally all results get sent to the R applica-

tion and they can be displayed graphically. The implementation of the parallelisation

was prepared by HANNES TSCHOFENIG. Then DR. ALBRECHT GEBHARDT and I

modified and extended this piece of software by adding more and more geostatistical

numerical stuff to the parallel communication framework.

Chapter 2

Spatial Statistics

2.1 Preliminaries

Spatial Statistics has its origin in the mining industry in the early fifties. It helped to

improve ore reserve calculation. The mining engineer DG KRIGE and the statistician

HS SICHEL made the first steps in South Africa. In the late fifties G MATHERON

adopted the concepts of DG KRIGE and developed them into a separate field of statis-

tics. First Spatial Statistics was only used for solving ore reserve estimation problems

but then in the seventies with the advent of high-speed computers it spread into wider

areas of earth sciences. Nowadays it is popular in many fields of science and industry

where there is a need for evaluating spatially or temporally correlated data.

Spatial Statistics takes into account the natural dependency of spatial data. The result-

ing correlation can be computed with avariogram and then it is used for increasing

the exactness of the model and for improving the local prediction of spatial data, called

Kriging .

6

CHAPTER 2. SPATIAL STATISTICS 7

2.1.1 Spatial Data

Spatial data consists of measurements bound to geographical locations. In addition to

the measured values they also contain the locations, resp. the coordinates of these data

values. The data array, in its most general form, may have the following shape

t1 x11 x2

1 x31 z1

1 · · · zi1 · · · zN

1

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

tα x1α x2

α x3α z1

α · · · ziα · · · zN

α

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

tn x1n x2

n x3n z1

n · · · zin · · · zN

n

In this arraytα and (x1α, x2

α, x3α) are the coordinates andzi

α are the variables. The

samples may have been taken at different momentstα at the locations with coordinates

(x1α, x2

α, x3α). Thenzi

α are different measurements ofN variables. The samples are

numbered using an indexα, whereα = 0, 1, 2, . . . , n. The total number of samples is

n. The different variables are labeled using an indexi, wherei = 0, 1, 2, . . . , N . The

number of variables isN .

2.1.2 Regionalized Variable

Now we assume that only one property has been measured (N = 1) and that we have

not recorded the time at which the measurement was made. So the indexi and the time


coordinatetα can be omitted. We describe then observations as

z(xα) with α = 1, 2, . . . , n .

The sampled objects taken from a regionD, whereD is a fixed subset ofRd, are a

fragment of a larger collection of objects. The possibility of more observations of the

same kind leads to theregionalized variableas

z(x) with x ∈ D .

The data setz(xα), α = 1, 2, . . . , n is a family of values of the regionalized variable.

Each measurement in the data set is a regionalized value and a sampled valuez(xα)

represents one draw from the random variableZ(xα)(see [15, page 26]).

z(x) with x ∈ D can be seen as one draw from an infinite family of random vari-

ables constructed at all pointsx of a spatial domainD which is called therandom

functionZ(x). I.e. the regionalized variablez(x) is a realization of a random function

Z(x). Now we can introduce arandom field(or random process)

Z(x) : x ∈ D

with its realizationz(x) : x ∈ D (see [4, p. 8]).

2.2 Covariance Structure

2.2.1 Stationarity

Definition 2.1: A random functionZ(x) is strictly stationaryif for any set ofn points

x1, · · · ,xn and any vectorh holds

Fx1,··· ,xn(z1, · · · , zn) = Fx+h1,··· ,x+hn(z1, · · · , zn) . (2.1)


This means if we translate a point configuration in a given direction, the multiple dis-

tribution does not change.

Definition 2.2: A random functionZ(x) is second-order stationaryif for a vectorh

linking any two pointsx andx + h in the domainD

E [Z(x + h)] = E [Z(x)] and cov [Z(x + h), Z(x)] = C(h) . (2.2)

The expectation and the covariance are translation invariant. The expected value

E [Z(x)] = m is the same at any pointx of the domain and the covariance between

any pair of locations depends only on the vectorh.

Definition 2.3: A random functionZ(x) is intrinsically stationaryif for a vectorh

linking any two pointsx andx + h in the domainD

E [Z(x + h)− Z(x)] = m(h) = 0 and var [Z(x + h)− Z(x)] = 2γ(h) . (2.3)

So we can conclude if a random functionZ(x) is second-order stationary,Z(x) is also

intrinsically stationary but vice versa this is not correct.

The quantity2γ(h) is known as the variogram and is the crucial parameter of spatial

statistics (see [4, page 40]). Thus with the properties of an intrinsically stationary

random function we yield the variogram2γ(h), respectively the semivariogramγ(h).

2.2.2 Simple Kriging

Under the assumption of second-order stationarity we want to build a weighted average

to make an estimation of a value at a pointx0 using information at the pointsxα, α =

1, · · · , n. With a meanm supposed constant over the whole domain, simple kriging


is used to estimate residuals from this valuem given a priori. Simple kringing is

also calledkriging with known mean. The definition of the weighted average for the

estimation is

Z?(x0) = m +n∑

α=0

wα (Z(xα)−m)

wherewα are the weights andZ(xα)−m = Z(xα)− E(Z(xα)) are the residuals.

If the estimation errorZ?(x0)− Z(x0) is nil the estimator is unbiased.

E [Z?(x0)− Z(x0)] = m +n∑

α=0

wα (E[Z(xα)]−m)− E[Z(x0)]

= m +n∑

α=0

wα(m−m)−m

= 0

The variance of the estimation errorσ2E is defined as

σ2E = var(Z?(x0)− Z(x0)) = E

[(Z?(x0)− Z(x0))

2]

.

Now we want to minimize the estimation variance by setting the first derivative to zero

∂σ2E

∂wα

= 0 for α = 1, · · · , n .

Thus the equation system for simple kriging can be written as

α = 1, · · · , n

n∑

β=1

wβC(xα − xβ) = C(xα − x0)

.

The solution of the system yields the optimal kriging weightswα.

The resulting optimal variance for each locationx0 is given by (see [15, page 21])

σ2SK = C(0)−

n∑α=1

wαC(xα − x0) .

2.2.3 Semivariogram and Covariance Function

The dissimilarity between pairs of data values,zα andzβ is measured by

γ?αβ =

(zα − zβ)2

2.


Under assumption of intrinsic stationarity the dissimilarityγ? only depends on the

orientation of the point pair which is described by the vectorhαβ.

γ?(hαβ) =1

2(z(xα + h)− z(xα))2 .

A plot of dissimilaritiesγ? against the spatial separationh is calledsemivariogram

cloud.

Now we form an average of dissimilaritiesγ?(hαβ) for a given class of vectors by

grouping allnc point pairs that can be linked by a vectorh. This class of vectors is

calledH. Generally non overlapping classesHk are chosen. So we can compute the

average dissimilarity as a value which is calledexperimental semivariogram.

γ?(Hk) =1

2nc

nc∑α=1

(z(xα + hαβ)− z(xα))2 with hαβ ∈ Hk

Remind thatγ?(Hk) is unbiased only if the assumption of instrinsic stationarity holds.

The behavior at very small scales, near the origin of the semivariogram is of impor-

tance. If the semivariogram is not differentiable at the origin we have a symptom called

nugget-effect. The values of the variable change abruptly at this very small scale. The

name of this effect has its origin in mining industry where it is believed that small

nuggets are causing this discontinuity at the origin.

Now we replace the experimental semivariogram by thetheoretical semivariogram.

The theoretical semivariogram guarantees that the variance of any linear combination

of sample values is positive. This fact is relevant for a kriging system. With the propo-

sition of intrinsic stationarity we yield the definition for the theoretical semivariogram

γ(h) =1

2E[(Z(x + h)− Z(x))2] . (2.4)


The value of the semivariogram at the origin is zero (γ(0) = 0) and all other values

are positive (γ(h) ≥ 0).

We obtain thecovariance functionC(h) from the hypothesis of stationarity of second

order of the random function

E [Z(x)] = m and E [Z(x) · Z(x + h)]−m2 = C(h) (2.5)

for all x,x + h ∈ D. Thecorrelation functionwe get by

ρ(h) =C(h)

C(0).

So far we can deduce a semivariogram function from a covariance function by

γ(h) = C(0)− C(h) . (2.6)

A covariance functionC(h) is apositive semidefinitefunction. This is important for

solving a Kriging system. For any set of pointsxα and any set of weightswα (or

written as a vectorw) the following equation is valid.

var

(n∑

α=0

wαZ(xα)

)=

n∑α=0

n∑β=0

wαwβC(xα − xβ) = w>Cw ≥ 0

A semivariogram is permissible if it is aconditionally negative definitefunction. For

any matrixΓ of semivariogram values it holds

w>Γw ≤ 0 forn∑

α=0

wα = 0 .

2.2.4 Semivariogram Models

We have mentioned that experimental semivariograms have to be fitted to theoretical

semivariogram models. Therefore we use the following ones:


• Spherical Model- only valid inR1, R2 andR3

γ =

0 , ‖h‖ = 0

c0 + c1

[32‖h‖a− 1

2

(‖h‖a

)3]

, 0 < ‖h‖ ≤ a

c0 + c1 , ‖h‖ ≥ a

(2.7)


0 500 1000 1500

050

000

1000

0015

0000

2000

00

x

gam

ma

(x)

range

sill

nugget

Spherical Semi−Variogram Model

Figure 2.1: Spherical Model


• Exponential Model - valid in all dimensions

γ =

0 , ‖h‖ = 0

c0 + c1

(1− e−3‖h‖/a

), ‖h‖ 6= 0

(2.8)

0 500 1000 1500

050

000

1000

0015

0000

2000

00

x

gam

ma

(x)

asymptotical range

asymptotical sill

nugget

Exponential Semi−Variogram Model

Figure 2.2: Exponential Model


• Gaussian Model- valid in all dimensions

γ =

0 , ‖h‖ = 0

c0 + c1

(1− e−3(‖h‖/a)2

), ‖h‖ 6= 0

(2.9)

0 500 1000 1500

050

000

1000

0015

0000

2000

00

x

gam

ma

(x)

asymptotical range

asymptotical sill

nugget

Gaussian Semi−Variogram Model

Figure 2.3: Gaussian Model


• Power Model- valid in all dimensions

γ =

0 , ‖h‖ = 0

c0 + b‖h‖λ , ‖h‖ 6= 0

(2.10)

with c0 ≥ 0, c1 ≥ 0, a ≥ 0 and 0 ≤ λ < 2, b ≥ 0.

If we setλ = 1 we get thelinear model. The power and the linear model are models

without sill.

All these models are described by the following three parameters:

c0 = nugget-effect

c0 + c1 = sill

a = range

The range correspondes to the distanceh whereγ(h) grows to the maximum. Obser-

vations which lie beyond the range, can be considered uncorrelated (γ(h) = 0.95(c0 +

c1)). The nugget-effect (c0 = limh→0 γ(h)) describes the microvariability of the val-

ues. The assymtotical range for the exponential semi-variogramm model is3a and for

the gaussian semi-variogramm√

3a.

2.2.5 Anisotropy

If the experimental semivariogram has different behavior in different directions, this

behavior is calledanisotropy. In this case we must try to obtain anisotropic random

functions from the isotropically defined semivariogram models by transforming the

coordinates (see [15, pages 46-48]).


In 2-D-space the behavior of the experimental semivariogram can be represented by

drawing a map of iso-semivariogram lines as a function of a vectorh. If the iso-

semivariogram lines are circular around the origin the phenomenon is isotropic. The

semivariogram depends only on the length fo the vectorh. If not, the iso-semivariogram

lines can be approximated by concentric ellipses. This type of anisotropy is calledgeo-

metric anisotropy. We can relate the class of ellipsoidally anisotropic random functions

to a corresponding isotropic random function by a linear transformation of the spatial

coordinates.

It can also happen that experimental semivariograms calculated in different directions

have a different value for the sill. If this behavior occurs we speak ofzonal anisotropy.

For example, in 2-D-space the sill along thex2 coordinate could be much larger than

alongx1.

2.3 Kriging

2.3.1 Estimation of the Mean

The mean value of spatial samples can be computed either with the arithmetic mean

or with the weighted average. The weighted average takes into account the knowledge

of the spatial correlation of the samples. A first approach for estimatingm by m? is to

use thearithmetic mean

m? =1

n

n∑α=1

Z(xα) . (2.11)

The problem is that the samples cannot be considered independent, however they are

correlated. A second approach to estimatem? is to use aweighted average

m? =n∑

α=1

wαZ(xα) (2.12)


wherewα are the weights. For∑n

α=1 wα = 1 we have no bias, i.e.E[m? − m] = 0,

wherem? is the estimated value andm = E[Z(x)] is the true value.

Now we want to define the ”best” weights. This means that we must reduce as much

as possible the variance of the estimation error. Thus we search the

minimum of var(m? −m) subject ton∑

α=1

wα = 1 .

For solving this problem we are using the method of Lagrange.

ϕ(wα,µ) = var(m? −m)− 2µ

(n∑

α=1

wα − 1

)whereµ is the Lagrange multiplier.

2.3.2 Ordinary Kriging

Ordinary Kriging is a method for estimating a value at a point of a region. In this

section the meanm is unknown. This method is often associated with the acronym

B.L.U.E. for ”best linear unbiased estimator” (see [9, page 278]). Ordinary kriging

is ”linear” because its estimates are weighted linear combinations of available data; it

is ”unbiased” since it tries to have the mean residual or error equal to 0; it is ”best”

because it aims at minimizing the variance of the errors. Ordinary kriging can also be

used for estimating a block value. With the assumption of second-order stationarity

we can evaluate the mean in a moving neighbourhood. First the kriging estimate of the

mean is set up and so a kriging estimator is examined.

We want to estimate a value atx0. Therefore we use the data values from then sample

pointsxα with the weightswα

Z?(x0) =n∑

α=1

wαZ(xα) wheren∑

α=1

wα = 1 . (2.13)


We must also assume that the data are a realization of a second-order stationary ran-

dom function with the covariance functionC(h). Because of second-order stationarity

it is guaranteed that the estimation is unbiased, i.e.E [Z?(x0)− Z(x0)] = 0.

The estimation variance is

σ2E = var(Z?(x0)− Z(x0)) (2.14)

= E[(Z?(x0)− Z(x0))

2]= C(x0 − x0) +

n∑α=1

n∑β=1

wαwβC(xα − xβ)− 2n∑

α=1

wαC(x0 − xα)

= C(0) + w>Kw − 2c>0 w

whereK = (C(xα − xβ))α,β=1,··· ,n , c0 = (C(x0 − x1), · · · , C(x0 − xn)) and

w = (w1, · · · , wn)> .

Now we wish to minimize the estimation variance (minC(0) + w>Kw − 2c>0 w)

with the constraint on the weights. We set the first derivative to zero. Thus we obtain

theordinary kriging system(OK) K 1

1> 0

w

µ

=

c0

1

. (2.15)

The vectorw contains the weights andµ is the Lagrange parameter.

We can say thatZ?(x0) =∑n

α=1 wαZ(xα) is the OK-prediction if and only if the

weightsw1, · · · , wn are solutions of the OK-system. Therefore we invert the positive-


definite matrix

K 1

1> 0

. We get

w

µ

=

K 1

1> 0

−1 c0

1

.

The ordinary kriging variance is

σ2OK = C(0)− c>0 w − µ .

If we assume intrinsic stationarity of the random function and not second-order sta-

tionarity, we can compute with the semivariogramγ(·) instead of the covariance ma-

trix C(·). This is allowed because of the equationγ(h) = C(0) − C(h). Then the

OK-system is given by Γ 1

1> 0

w

µ

=

γ0

1

(2.16)

whereΓ = (γ(xα − xβ))α,β=1,··· ,n , γ0 = (γ(x0 − x1), · · · , γ(x0 − xn))> and

w = (w1, · · · , wn)> .

So the ordinary kriging variance is

σ2OK = w>γ0 − µ .

Ordinary kriging has the property of anexact interpolator. If x0 is identical with the

data locationxα then

Z?(x0) = Z(xα) if x0 = xα .

To estimate a block value instead of a point value we can also use ordinary kriging.

We call it thenblock kriging. Instead of a point we use a blockB(x0) ⊂ D with center


x0. In order to krige a block we modify our ordinary kriging model in the following

manner K 1

1> 0

w

µ

=

CB

1

(2.17)

whereCB = (C(x1, B), · · · ,C(xn, B))> with C(xα, B) = 1|B|

∫B

C(xα − x) dx.

CB(xα, B) with xα ∈ D is called the ”point-to-block-covariance”.

For solving the block kriging system we also need the ”block-covariance” between

block B1 and blockB2. C(B1, B2) = 1|B1||B2|

∫B1

∫B2

C(x − y) dx with x ∈

B1 and y ∈ B2.

2.3.3 Universal Kriging

In this section the meanm is unknown and we have non-stationary random functions.

We involve a non-stationary component called thedrift (see [15, page 177]). We con-

sider functionsfl of the coordinatesx deterministic on a domainD if they are known

at any location of the domain.

The random functionZ(x) is composed by

Z(x) = m(x) + Y (x)

wherem(x) is the deterministic drift andY (x) is a residual random function. Assum-

ing E [Y (x)] = 0 we haveE [Z(x)] = m(x).

The driftm(x) depends onfl in a linear way, such as

m(x) =L−1∑l=0

alfl(x) = f(x)>a


where the coefficient vectora = (a0, · · · , aL−1)> is different from zero andf(x) =

(f0(x), · · · , fL−1(x))> is the function vector. The functionf0(x) is defined byf0(x) =

1.

Now we want to set up auniversal kriging(UK) system. We start at the linear combi-

nation

Z?(x0) =n∑

α=1

wαZ(xα) . (2.18)

With no biasE [Z(x0)− Z?(x0)] = 0 we yield

m(x0)−n∑

α=1

wαm(xα) = 0

and setting form(x) =∑L−1

l=0 alfl(x) we get

L−1∑l=0

al

(fl(x0)−

n∑α=1

wαfl(xα)

)= 0 .

Thus we obtainn∑

α=1

wαfl(xα) = fl(x0) for l = 0, · · · , L− 1 .

So we get the system

F>w = f(x0) (2.19)

with F = (fl(xα)) α=1,2,··· ,nl=0,1,··· ,L−1

=

1 f1(x1) · · · fL−1(x1)

......

......

1 f1(xn) · · · fL−1(xn)

.

Further we wish to minimize the estimationσ2E variance with the constraintF>w =

f(x0). Using this constraint we can say thatZ?(x0) is unbiased.

minσ2E(x0) = minvar(Z?(x0)− Z(x0))

= minvar(n∑

α=1

wαZ(xα)− Z(x0))


With Z(x) = m(x) + Y (x) and the constraints on the weights∑n

α=1 wα = 1 we get

minσ2E(x0) = minvar(

n∑α=1

wαY (xα)− wαY (x0))

= minvar(n∑

α=1

wα(Y (xα)− Y (x0))) .

If we compute the variance and set the problem into vector notation we get

min2w>γ0Y −w>ΓY w

with γ0Y = (γY (x0 − x1), · · · , γY (x0 − xn))>, ΓY = (γY (xα − xβ))α,β=1,2,··· ,n,

where γY (·) is the semivariogram fromY (·) = Z(·)−m(·).

By using the Lagrange parametersµl we set the first derivate to null. We must also

take into consideration the constraintF>w = f(x0). Thus we obtain the UK-system Γ F

F> 0

w

µ

=

γ0Y

f(x0)

(2.20)

To guarantee a solution for this system the matrixF must have linearly independent

column vectors andΓ must be conditionally negative definite.

If wandµ are solutions of the UK-system andγY is the semivariogram of the residuals

thenZ?(x0) =∑n

α=1 wαZ(xα) is the best linear unbiased prediction and the minimal

estimation variance isσ2UK(x0) = w>γ0

Y + µ>f(x0).

If we use the covariance matrix of the residualsY (·) andCY (·) instead of the semivar-

iogramγY (·) the UK-system looks like CY F

F> 0

w

µ

=

c0Y

f(x0)

(2.21)


where

c0Y = (cY (x0 − x1), · · · , cY (x0 − x1)) and CY = (Cy(xα − xβ)α,β=1,2,··· ,n) .

(2.22)

Chapter 3

Kriging on Tiled Grids

3.1 Solving Kriging Systems Simultaneously

Usually kriging prediction is performed on regular rectangular grids. The parameters

of these grids (grid spacing inx andy direction) control the quality of the estimated

map. For that reason fine grids are prefered by the user. Obviously computational

burden increases with the number of points in the grid because a system of linear

equations (see equation 2.21) for each grid point has to be solved.

For each grid pointx0 we have a system matrix CY F

F> 0

(3.1)

and a right hand side c0Y

f(x0)

(3.2)

associated withx0 and its kriging search neighbourhood via equation 2.22.

26

CHAPTER 3. KRIGING ON TILED GRIDS 27

With most of the linear equation system (LES) solvers it is possible to solve a whole

family of linear equations systems of the form

Axi = bi for i = 0, · · · , n (3.3)

simultaneously. I.e. they solve the matrix equation

AX = B

with X = (x1, · · · ,xn) and B = (b1, · · · ,bn).

This solutions are realized e.g. in the LAPACK routine DGESV(see [1]), which ap-

plies a LU-factorization toA = PLU . P is a permutation matrix (for pivoting),L and

U are lower resp. upper triangular matrices. These factors are used to determine the

solutionX. Thus, if it is possible to combine the solving steps of all equations (see

equation 3.3) in one single step, the same factorizationA = PLU can be used for all

right hand sides. Compared with solving each of the equations separately, this idea

will reduce computational time.

This approach can be applied to kriging prediction in the following way. Instead of

solving each krige system(see equation 2.21) for a set of pointsx0i(i = 1, · · · , n) by

establishing individual search neighbourhoodsSi = xi | ‖xi − x0i‖ ≤ r and calcu-

lating CY,Si, FSi

, c0Yi

andf(x0i). We take the union of all search neighbourhoodsSi

and form an overall search neighbourhood∪iSi. We use the covariance matrixCY,∪iSi,

the design matrixF∪iSiand the vectorsc0

Yiandf(x0i) to establish the following krige

system CY,∪iSiF∪iSi

F>∪iSi

0

w1 · · · wn

µ1 · · · µn

=

c0Y 1 · · · c0

Y 1

f(x01) · · · f(x01)

. (3.4)

Now we can use e.g. the LAPACK routine DGESV to solve this system in one step.


If we choose the set of pointsx0i in a way such that the individual search neighbour-

hoods overlap to a large extent, the following inequation will hold.

# ∪i Si ∑

#Si (3.5)

This means that the matrixCY,∪iSiis slightly greater than the individualCY,Si

but it is

not as large as a block matrix composed of all theCY,Si. So we will be able to solven

problems in one single step by the cost of a slightly enlarged system matrix.

In order to achieve the goal of overlapping search neighbourhoods we have to choose

x0i which are close together.

3.2 Tiled Grid Kriging

Practically we will generalize the rectangular grid as follows.

Definition 3.1: A tuple

〈xij | i = 1, · · · , n, j = 1, · · · , m, tkl | k = 1, · · · , r, l = 1, · · · , s, dx, dy, tx, ty〉

is called atiled grid if it consists of a set of pointsxij ∈ R2 where

xij =

xij

yij


S1 S1_2

x1

. x2

.

x1

.x2

.

. .

. .

. .

. .

. .

. .

... .

. ...

. .

S2

Figure 3.1: Overlapping Search Neighbourhoods


with

yij+1 − yij = dy

yi+1j = yij

xi+1j − xij = dx

xij+1 = xij

with i = 1, · · · , n and j = 1, · · · , m and a partition of the setxij | i = 1, · · · , n, j =

1, · · · , m in tiles tkl | tkl ⊆ xij | i = 1, · · · , n, j = 1, · · · , m, xij ∈ tkl, k =

1, · · · , r, l = 1, · · · , s with

xij ∈ tkl

if and only if

ktx ≤ i ≤ (k + 1)tx

lty ≤ j ≤ (l + 1)ty .

In this notationdx anddy are the grid spacings inx andy direction. tx andty are the

tiling parameters.

So one tile is a rectangular part oftx × ty grid points (see figure3.2).

Now it is straigthforward to adopt the usual kriging technology to the tiled grids, using

the following algorithm:


. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

Figure 3.2: Grid Points in Tiles


• determinexij andtkl

• for eachtkl:

– for eachxij in tkl:

∗ Sij = searchneighbourhood forxij

∗ determinec0Y ij andf(x0ij)

– Skl = ∪ijSij

– determineCkl andF kl

– determineC0klY =

((c0

Y ij))

ij, F kl

0 = ((f(x0ij)))ij

– determineK =

CklY F kl

F>kl 0

, RHS =

C0klY

F kl0

– call DGESV for

KX = RHS

– extractλij from X =

Λ

M

with Λ = ((λij))ij andM =((µij)

)ij

– calculatezij = λ>ijz (z : zα = z(xα) for xα ∈ Skl = data set restricted to

searchneighbourhood)


This algorithm is implemented in the functionkrige.tiles of the libraryrgeostat .

This function provides an interface to the underlying FORTRAN functionkrgtil.f .

The ”tiling” process fits very natural in the kriging settings because it resembles the

idea of regionalized variables via collecting neighbouring grid points into small areas

of interest. But it has one drawback. At the boundaries of the tiles arise the well known

border effects. This finally leads to non-continous contourlines at these rectangular

borders(see figure3.6).

3.3 Tiled Grid Kriging with Overlapping

In this section I will pick up the border effects(see section3.2). One obvious solution

will be to allow the tiletkl to overlap up to some extent with its neighbours. This leads

to

Definition 3.2: A tuple

〈xij | i = 1, · · · , n, j = 1, · · · , m, tkl | k = 1, · · · , r, l = 1, · · · , s, dx, dy, tx, ty, ox, oy〉

is called anoverlapping tiled gridif it consists of a set of pointsxij ∈ R2 where

xij =

xij

yij


with

yij+1 − yij = dy

yi+1j = yij

xi+1j − xij = dx

xij+1 = xij

with i = 1, · · · , n and j = 1, · · · , m and a partition of the setxij | i = 1, · · · , n, j =

1, · · · , m in tiles tkl | tkl ⊆ xij | i = 1, · · · , n, j = 1, · · · , m,xij ∈ tkl, k =

1, · · · , r, l = 1, · · · , s with

xij ∈ tkl

if and only if

ktx − ox ≤ i ≤ (k + 1)tx

kty − oy ≤ i ≤ (k + 1)ty

ktx ≤ i ≤ (k + 1)tx

lty ≤ j ≤ (l + 1)ty

In addition to definition 3.2 now we have an overlapping ofox or oy points inx andy

direction. Now it is possible that a pointxij belongs to more than one tiletkl.

Now we can apply almost the same algorithm as in section3.2.


. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

Figure 3.3: Grid Points in Overlapping Tiles


• determinexij andtkl

• for eachtkl:

– for eachxij in tkl:

∗ Sij = searchneighbourhood forxij

∗ determinec0Y ij andf(x0ij)

– Skl = ∪ijSij

– determineCkl andF kl

– determineC0klY =

((c0

Y ij))

ij, F kl

0 = ((f(x0ij)))ij

– determineK =

CklY F kl

F>kl 0

, RHS =

C0klY

F kl0

– call DGESV for

KX = RHS

– extractλij from X =

Λ

M

with Λ = ((λij))ij andM =((µij)

)ij

– calculatezij = λ>ijz (z : zα = z(xα) for xα ∈ Skl = data set restricted to

searchneighbourhood)

• for i = 1, · · · , n

– for j = 1, · · · , m

– zij =

∑xi′j′∈tkl

zkli′j′

nijwherenij is the number of tiles containingxij

nij = #tkl|xij ∈ tkl


So we can determine the krige mapzij|i = 1, · · · , n, j = 1, · · · , m by averaging all

computedzklij .

Because of the addition theoremV ar(a(X + Y )) = a2V ar(X) + a2V ar(Y ) we can

calculate the estimated kriging variance in the following manner

σ2ij =

∑xi′j′∈tkl

σkl2

i′j′

n2ij

(3.6)

3.4 Example

Now it is possible to tune the overlapping result by the parametersdx, dy, tx, ty, and

ox, oy. Now we can avoid the ”checker board” output (compare figures3.5, 3.6 and

3.7) caused by border effects by enlarging the amount of overlapping points. On the

other hand, overlapping means that we calculate some estimates twice, four times or

even more often. So we have to find a compromise between output quality(dx, dy will

increase), avoiding border effects(ox, oy will increase) and time saving(tx, ty will in-

crease).

Figure3.4compares the computational speed of the three kriging functions:krige.grid ,

krige.tiles andkrige.tilesov . In order to get a representative sample the

time for each of the three functions was measured five times and then the mean was

taken. In the first graphic the functionkrige.tiles has tile size 5, which means

that every subrectangle contains 5 points. The functionkrige.tilesov has tile

size 5, too. This function has the modification that we can select how many points,

lying in the tiles, should overlap in order to yield exacter results. So we can see that

the functionkrige.tiles is much faster thankrige.grid but the computational

exactness decreases at the tile borders forkrige.tiles (compare figures3.5, 3.6

and3.7). Therefore the functionkrige.tilesov was built which computes more


exactly because of overlapping tiles. But the greater the number of overlapping points

the greater will be the computational effort. We can see this characteristic in all four

graphics of figure3.4.


50 100 150 200 250 300

05

1015

2025

30

tile size:5, overlapping points:2,3number of grid points

time

in s

econ

ds

krige.gridkrige.tileskrige.tilesov

50 100 150 200 250 300

05

1015

2025

30

tile size:10, overlapping points:2,4,6,7number of grid points

time

in s

econ

ds


50 100 150 200 250 300

05

1015

2025

30

tile size:20, overlapping points:2,4,8,12,14number of grid points

time

in s

econ

ds


50 100 150 200 250 300

05

1015

2025

30

tile size:30, overlapping points:4,8,12,16,20number of grid points

time

in s

econ

ds


Figure 3.4: Computational Burden of krige.grid, krige.tiles and krige.tilesov


178500 179000 179500 180000 180500 181000

3300

0033

1000

3320

0033

3000

Figure 3.5: Plot of maas.kgrid, where nx=ny=100


178500 179000 179500 180000 180500 181000

3300

0033

1000

3320

0033

3000

Figure 3.6: Plot of maas.ktiles, where nx=ny=100, itx=ity=10


178500 179000 179500 180000 180500 181000

3300

0033

1000

3320

0033

3000

Figure 3.7: Plot of maas.ktilesov, where nx=ny=100, itx=ity=10, nxper=nyper=2

Chapter 4

PVM (Parallel Virtual Machine)

4.1 Introduction

In this chapter a software package for developing parallel programs is presented. It is

executable on networked Unix computers. The tool is called Parallel Virtual Machine

(PVM) and it allows a heterogeneous collection of workstations to function as a sin-

gle high-performance parallel machine. PVM is portable and runs on a wide range of

modern platforms. PVM is well accepted by the global computing community. It is

used successfully for solving large-scale problems.

PVM is designed to link computing resources and provide users with a parallel plat-

form for running their computer applications. The number of different computers they

use does not matter and also the location of the computers is irrelevant. When PVM

is correctly installed, it is capable of harnessing the combined resources of typically

heterogeneous networked computing platforms in order to deliver high levels of per-

formance and functionality.

43

CHAPTER 4. PVM (PARALLEL VIRTUAL MACHINE) 44

The PVM project began in summer of 1989 at Oak Ridge National Laboratory (see

[6, preface]). The current version 3.4 is described in this chapter and with this version

(called simply PVM) I use to work. The PVM software has been distributed freely and

is being used in computational applications around the world.

Parallel processing is a method where we have many small tasks which solve a large

problem. It was developed because of the demand for high performance and low cost.

Generally there are two opportunities: massively parallel processors (MPPs) and dis-

tributed computing. MPPs (see [6, page 6]) are very powerful computers. These ma-

chines combine a few hundred to a few thousand CPUs in a single large cabinet con-

nected to hundreds of gigabytes of memory. The second development which concerns

scientific problem solving isdistributed computing. A set of cumputers are connected

by a network. So they are used collectively to solve a single large problem. The most

important reason why distributed computing is used are the costs. Large MPPs are

very expensive and in contrast users see very little cost in running their problems on a

local set of existing computers. In an MPP, every processor is exactly like every other

in capability, resources, software and communication speed. Not so on a network.

These computers are heterogeneous. They may be made by different vendors with dif-

ferent processors and memories and they may have different compilers. So the set of

computers can include a wide range of architecture types. PVM is designed to make

programming possible for a heterogeneous collection of machines. A key concept in

PVM is that it works on a collection of computers and they appear as one largevirtual

machine, hence its name.

4.2 PVM’s basics

The user composes his application as a collection of cooperating tasks. Tasks access

PVM resources through a library via standard interface routines. These routines allow


the initialization and termination of tasks across the network as well as communication

and synchronization between tasks. The PVM message-passing primitives are oriented

towards heterogeneous operation, involving strongly typed constructs for buffering and

transmission. Communication constructs include those for sending and receiving data

structures as well as high-level primitives such as broadcast, barrier synchronization,

and global sum(see [6, page 5]).

At any point in the execution of a concurrent application, any task may start or stop

other tasks or add or delete computers from the virtual machine. Any process may

communicate with any other.

4.3 The PVM System

The PVM System is composed of two parts. The first part is a daemon, calledpvmd3.

It runs on all the computers, the hosts. It serves as a message router and makes up

the virtual machine. Any user with a valid login can install the daemon on a machine.

When a user wishes to run a PVM application, he first creates a virtual machine by

starting up PVM. The PVM application can then be started from a Unix prompt on

any of these hosts. The second part of the system is a library of PVM interface rou-

tines. It contains a functionally complete repertory of primitives that are needed to

cooperate between tasks of an application. This library contains user-callable routines

for massage passing, spawning processes, coordinating tasks, and modification of the

virtual machine. The PVM currently supports C, C++, and Fortran languages. The

C and C++ language bindings for the PVM user interface library are implemented as

functions, but Fortran language bindings are implemented as subroutines. All PVM

tasks are identified by an integertask identifier(TID). It is an address. Messages are

sent to and received from those TIDs.


A user writes one or more sequential programs in C, C++, or Fortran 77 that contain

embedded calls to the PVM library. Each program corresponds to a task which makes

up the application. These programs are compiled for each architecture in the host

pool. Then the resulting object files are placed at a location accessible from machines

in the host pool. When a user wants to execute an application, he typically starts one

copy of one task (usually the ”master” task) by hand from a machine within the host

pool1. This process subsequently starts other PVM tasks. In the master-slave model

there is a separate ”control” program. It is called the master and it is responsible for

process spawning, initialization, collection and display of results. In contrast there are

slave programs which perform the actual computation that is involved. They are either

allocated in their workloads by the master or they perform the allocations themselves.

For example, when a user wants to multiplicate two matrices A and B (see [6, page

36, 76]) a group of processes is first spawned which means that a new process or PVM

task is created and it uses the master-slave paradigm. The matrices to be multiplied are

divided into subblocks. Then each subblock of the A and B matrices is placed on the

corresponding process. During computation subblocks need to be exchanged between

processes. At the end of the computation the resulting matrix subblocks are present

in the individual processes. Their position must be in agreement with their respective

positions on the process grid. It also must be consistent with a data partitioned map of

the resulting matrix C.

1usually$HOME/pvm3/bin/$PVMARCH

Chapter 5

Universal PVM Kriging

5.1 Combining R and PVM

We follow an earlier work(see [13]) written by HANNES TSCHOFENIG. He proved

the possibility of a successful combination of R and PVM. He implemented a function

following the lines of thekrige.all fuction of the librarysgeostat which works

as follows.

The function uses the whole data set as search neighbourhood for all points of the

prediction grid. It calculates the (ordinary) kriging system matrix

K =

C 1

1> 0

The function inverts the matrix and subsequently calculates the

λ

µ

solutions for

47

CHAPTER 5. UNIVERSAL PVM KRIGING 48

each grid point by λ

µ

= K−1

CY0

1

. (5.1)

Now the PVM implementationpvmkrige , which was the starting point of the li-

brary pvmkrige , distributes the matrix vector produced in equation 5.1 across the

participants of PVM in fixed size chunks of grid points. Of course this approach is

numerically not very robust. The reason is that it crucially depends on the condition of

the overall krige matrixK. This matrix gets very large in size (it depends on the data

set). Only a subsequent step (see equation 5.1) is involved in the parallel part of the

functionpvmkrige .

But the focus of this first implementation was not the numerical but rather the technical

aspect of connecting R with a PVM server process and exchanging data and results in

both directions. This structure can be seen in the following figure.


Rlibrary(pvmkrige)

krige_server

krige_client krige_client krige_client

P MV

Figure 5.1: PVM Computational Model


5.2 Parallel Tile Kriging

Now we extend the functionpvmkrige and apply it to the overlapping tile kriging

approach which is described in section 3.3. The task was to adopt the data transfer

functions because now much more data was needed by the clients. Additionally to

the variogram parameters, already needed bypvmkrige , the new function uses the

whole data set and needs the complete coordinates of the estimation grid.

Fortunately all these all these big chunks of data have to be transfered only once. Dur-

ing the calculation phase we only need to send the description of the tilestkl which

consist of the coordinates of the tile’s south-west and north-east corner, expressed in

the index numbering ofxij points. This description of the tiles consists of four integer

values which are to be sent. The result consists of2× tx × ty double precision values

(zij, σ2ij).

The following listings give an overview of the whole implementation. More details

can be found in the source code on the CD-ROM accompaning this work.

krige.pvm.R

• preparation of some parameters

• call of krige pvm.c

krige pvm.c

• creation of socket

• handshake with server


. . .

. . .

. . . . . . . .

. . . . .

. . .

. . .

. . .

. . .

. . . . . . . .

. . . . .

. . .

. . .

. . .

. . .

. . .

. . .

. . .

. . .

. . .

. . .

. . .

R SERVER

CLIENT A

CLIENT B

CLIENT C

Figure 5.2: Parallelizing the Tiles


• sends information about the data and static parameters to the server

• sends the datalon, lat, z to the server

• calculates the grid coordinates matricesxgwork, ygwork

• sends the grid coordinate matrices to the server

• receives the resultszg, varg from the server

krige server.c

• starts communication with R

• gets information about parameter settings

• receiveslon, lat, z, xgwork, ygwork

• starts PVM by callingpvm tiles.c

• sends the resultszg, varg

pvm tiles.c

• starts the clients

• sends the jobs to the clients

• receives the tile resultsz, var

• kills the clients

• builds the result matriceszg, varg

krige client.c

• receives the parameter settings


• receive the datalon, lat, z, xgwork, ygwork

• calculates the covariance matrix

• receives the jobs

• calls the FORTRAN functionkrige.f

• sends results back to the server


5.3 Example

In this section I tested how long different sets of computers need to computekrige.pvm .

A set of eight computers (five Alpha, three x86), a set of five computers (three Alpha,

two x86) and a set of two computers (one Alpha, one x86) got compared concerning

their computational speed. I testedkrige.pvm for 50 to 300 grid points. Two differ-

ent tile sizes (10, 30) and different numbers of overlapping points (2, 6 reps. 8, 16).

The result of the tests is given in figure5.3. A great difference can be seen between

two and five resp. eight computers. Whereas there is nearly no difference between

five and eight computers. For this kriging problem we see that it is not profitable to

take more than five computers. The reason is that in this example the size of data is

too small for distributing the jobs to more than five computers. I can suppose that the

receive - send effort is too great when too many computers are taken.


hostname architecture CPU speed memory OS

alpha1-stat Alpha EV4 200 MHz 128 MB Tru64



alpha4-stat Alpha EV6 667 MHz 1 GB Tru64


pc04-stat x86 Pentium 90 MHz 32 MB FreeBSD

pc05-stat x86 Pentium MMX 166 MHz 64 MB SunOS

pc10-stat x86 Dual PIII 500 MHz 512 MB Linux

Table 5.1: Architectures used in the PVM setup


50 100 150 200 250 300

050

100

150

tile size:10, overlapping points:2number of grid points

time

in s

econ

ds

8 computers5 computers2 computers

50 100 150 200 250 300

050

100

150


time

in s

econ

ds


50 100 150 200 250 300

050

100

150


time

in s

econ

ds


50 100 150 200 250 300

050

100

150


time

in s

econ

ds


Figure 5.3: Computational Effort with different sets of computers


The following figure shows the send - receive interaction between server and clients as

displayed by XPVM1.

1a graphical tool written intcl/tk , distributed with the PVM software


Figure 5.4: XPVM: Communication between Server and Clients

Chapter 6

Conclusion and Prospect

The goal of this work was to develop and implement a parallel strategy of performing

kriging prediction. This was achieved by partitioning the estimating grid into rectan-

gular tiles and by using the parallel virtual machine (PVM) to solve the linear equation

systems associated with these tiles. Using the free open source software R and PVM

it was possible to implement this strategy in a very convenient manner. Finally it was

possible to build an add-on package for R which provides an end-user interface to the

above described algorithms.

Now it should be possible to use this work for subsequent tasks, e.g. for high quality

maps with increasing numbers of grid points or for simulation tasks where kriging

maps have to be calculated repeatedly. Especially, in the field of experimental design

whereσ instead ofz is the object of interest, expensive optimization algorithms are

needed which evaluate manyσ maps for a large set of possible monitoring networks

(set of locations wherez is to be measured). Here parallel computing could help to

improve the performance of the optimization routines.

59

Appendix A

Functions

All three kriging functions,krige.grid , krige.tiles andkrige.tilesov

were developed using the freely available statistical program R. R can be found in the

internet on the CRAN Servers (Comprehensive R Archive Network), for example on

http://www.r-project.org .

The following pages show the online documentation of the functions for kriging on

grids, tiled grids and tiled grids with overlapping points. These functions are imple-

mented in the libraryrgeostat . Finally we will give some advice for using the func-

tion krige.pvm of the librarypvmkrige . All these libraries will be available at

ftp://ftp-stat.uni-klu.ac.at/pub/R/contrib and they are included

on the accompaning CD-ROM.

60

APPENDIX A. FUNCTIONS 61

A.1 Kriging Functions

A.1.1 The Functionkrige.grid

Description

Performs (2D) universal kriging on a given rectangular grid using data and variogram

objects prepared withpoint , pair , est.variogram andfit.variogram from

library sgeostat .

Usage

krige.grid(point.obj, at, var.mod.obj, xsw=NULL, ysw=NULL,

xne=NULL, yne=NULL, dx=NULL, dy=NULL, nx=NULL,

ny=NULL, angle=NULL, maxdist=NULL, extrap=F,

trend=0, rsearch=0, nsearch=0, nsmin=-1,

nsmax=-1, mode=3, duplicate = "error",

dupfun = NULL)

Arguments

point.obj spatial data, class ”point ”

at name of variable of interest inpoint.obj

var.mod.obj variogram model, generated byfit.variogram

xsw x coordinate of south west corner of the grid

ysw y coordinate of south west corner of the grid

xne x coordinate of north east corner of the grid


yne y coordinate of north east corner of the grid

dx grid spacing in x direction

dy grid spacing in y direction

nx grid dimension in x direction, to be used in absence ofdx

ny grid dimension in y direction, to be used in absence ofdy

angle rotation of grid – not yet implemented –

maxdist synonym forrsearch , for compatibility withsgeostat

extrap logical, indicates whether to extrapolate beyond convex hull of data points,

defaultFALSE

border optional polygon (list with two componentsx andy of same length) rep-

resenting a (possibly non convex) region of interest to be used instead of the

convex hull. Needsextrap=FALSE .

trend order of trend function

rsearch fixed search radius

nsearch fixed nearest neighbour search

nsmin minimum no of points in search neighbourhood, default 0

nsmax maximum no of points in search neighbourhood, defaultlength(point.obj[at])

mode 1 - calculate only krige prediction, 2 - calculate only krige variance (useful for

network design purposes), 3 - both (default)

duplicate indicates how to handle duplicate data points. Possible values are"error"

- produces an error message,"strip" - remove duplicate z values,"mean" ,"median" ,"user"

- calculate mean , median or user defined function of duplicate z values.


dupfun this function is applied to duplicate points ifduplicate="user"

Details

Trend models are coded via thetrend parameter (only up to order 2 supported):

trend trend function

0 z ˜ 1

1 z ˜ x+y+1

2 z ˜ x+y+I(x ˆ 2)+I(y ˆ 2)+x*y

The krige system is solved via the LAPACK routine DGESV.

If mode is set to 2 (default is 3) only kriging variances are calculated. In this case the

values inat (the measurements) are not used. This may be useful for experimental

design applications where only new locations are given and measurements do not yet

exist. Of coursevar.mod has to be given (computed from an existing dataset).

Value

An object of typekrige.map .

Note

In contrast to the interpretedkrige function from librarysgeostat this function is

written in FORTRAN, so it should be reasonably faster.

Author A. Gebhardt, [email protected]

Referencesftp://ftp-stat.uni-klu.ac.at/pub/R/contrib

SeeAlsokrige.points ,krige

Examples


# compare result and system.time() with example(krige) from library(sgeostat)

system.time(

maas.kgrid <- krige.grid(maas.point, "zinc", maas.vmod,

min(maas.bank$x), min(maas.bank$y),max(maas.bank$x), max(maas.bank$y),

nx=100,ny=100, rsearch = 600, extrap = F, border=maas.bank)

)

# to get a picture type "example(plot.krige.map)"

A.1.2 The Functionkrige.tiles

Description

Performs (2D) universal kriging on a given grid (divided in subrectangles (tiles)) us-

ing data and variogram objects prepared withpoint , pair , est.variogram and

fit.variogram from librarysgeostat .

Usage

krige.tiles(point.obj, at, var.mod.obj, xsw=NULL, ysw=NULL,


ny=NULL, itx=NULL, ity=NULL, angle=NULL, maxdist=NULL,

extrap=FALSE, border=NULL, trend=0, rsearch=0,

nsearch=0, nsmin=-1, nsmax=-1, mode=3,

duplicate = "error", dupfun = NULL)

Arguments













itx number of grid points per tile in x direction

ity number of grid points per tile in y direction




defaultFALSE
















Details



0 z ˜ 1

1 z ˜ x+y+1

2 z ˜ x+y+I(x ˆ 2)+I(y ˆ 2)+x*y

In contrast tokrige.grid the grid specified forkrige.tiles is built by rectan-

gular subregions (tiles). Kriging predictions for points belonging to the same tile are

calculated in a single step. The intention was to reduce computational burden by col-

lecting some neighbouring krige systems and forming a combined krige system (using

multiple right hand sides in the krige equations, see LAPACK routine DGESV).






Value

An object of typekrige.map . A list (identical tokrige.grid ):

x vector of x coordinates of grid

y vector of y coordinates of grid

z matrix of predicted values, see alsomode

var matrix of prediction variances

tiles

snb matrix with 0-1 search neighbourhood codes, columns refer to grid poiuts, rows

to data points

Note

Drawback: border effects between tiles, not only at the outer regions of the grid. This

tile approach can be viewed as intermediate step to a block kriging predictor, it is only

left to average predictions across tiles (identifying tiles with blocks).

Author A. Gebhardt, [email protected]


SeeAlsokrige.cell ,krige

Examples

# compare result and system.time() with example(krige.grid)

system.time(

maas.ktiles <-

krige.tiles(maas.point, "zinc", maas.vmod,

min(maas.bank$x), min(maas.bank$y), max(maas.bank$x),


max(maas.bank$y), nx=100,ny=100, itx=10, ity=10,

rsearch = 600, extrap = F, border=maas.bank)

)

# show a picture with tile containing point(50,50) marked :

plot(maas.ktiles,show.snb=TRUE,snb.i=50,snb.j=50)

A.1.3 The Functionkrige.tilesov

Description

Performs (2D) universal kriging on a given grid (divided in subrectangles (tiles) which

can overlap) using data and variogram objects prepared withpoint , pair , est.variogram

andfit.variogram from librarysgeostat .

Usage

krige.tilesov(point.obj, at, var.mod.obj, xsw=NULL, ysw=NULL,


ny=NULL, itx=NULL, ity=NULL, nxper=NULL,nyper=NULL,

angle=NULL, maxdist=NULL, extrap=FALSE, border=NULL,

trend=0, rsearch=0, nsearch=0, nsmin=-1, nsmax=-1,

mode=3, duplicate = "error", dupfun = NULL)

Arguments















nxper number of overlapping grid points per tile in x direction

nyper number of overlapping grid points per tile in y direction




defaultFALSE
















Details



0 z ˜ 1

1 z ˜ x+y+1

2 z ˜ x+y+I(x ˆ 2)+I(y ˆ 2)+x*y

In contrast tokrige.tiles the tiles specified forkrige.tilesov can overlap.

The dimension of the overlap is defined bynxper andnxper . Kriging predictions

for points belonging to the same tile are calculated in a single step. Then the mean

of the Kriging predictions for the tile points which are computed more than once be-

cause of overlapping tiles, is taken. The intention was to enhance the computational

exactness and to eliminate the border effects between tiles by computing the marginal

points of a tile more often.






Value

An object of typekrige.tiles .

Note

In contrast tokrige.tiles the border effects between tiles are eliminated.

Author N. Samonig, [email protected]


SeeAlsokrige.cell ,krige

Examples

# compare result and system.time() with example(krige.tiles)

system.time(

maas.ktilesov <-

krige.tilesov(maas.point, "zinc", maas.vmod,



nxper=2, nyper=2, rsearch = 600,

extrap = F, border=maas.bank)

)


plot(maas.ktilesov,show.snb=TRUE,snb.i=50,snb.j=50)


A.2 PVM Kriging Function

A.2.1 The Functionkrige.pvm

Description

Performs (2D) parallel universal kriging on a given grid (divided in subrectangles

(tiles) which can overlap) using data and variogram objects prepared withpoint ,

pair , est.variogram and fit.variogram from library pvmkrige . Small

computational tasks can be distributed over different computers which work in a par-

allel way.

Usage

krige.pvm(point.obj, at, var.mod.obj, xsw=NULL, ysw=NULL,


ny=NULL, itx=NULL, ity=NULL, nxper=NULL,nyper=NULL,

angle=NULL, maxdist=NULL, extrap=FALSE, border=NULL,

trend=0, rsearch=0, nsearch=0, nsmin=-1, nsmax=-1,

mode=3, duplicate = "error", dupfun = NULL)

Arguments















nxper number of overlapping grid points per tile in x direction

nyper number of overlapping grid points per tile in y direction




defaultFALSE
















Details



0 z ˜ 1

1 z ˜ x+y+1

2 z ˜ x+y+I(x ˆ 2)+I(y ˆ 2)+x*y

In contrast tokrige.tilesov the tiles specified forkrige.pvm get parallelized.

The tile jobs (Kriging predictions for points belonging to a tile) are distributed to dif-

ferent computers. The intention was to enhance the computational speed by paral-

lelization.

In the current implementation you have to start the server processkrige server

manually. Don’t forget to add the linepvmkrige 50500 to your file/etc/services .






Value

An object of typekrige.tiles .

Note

In contrast tokrige.tilesov the result of the function gets computed by more

than one computer.

Authors

A. Gebhardt, [email protected]

N. Samonig, [email protected]

H. Tschofenig, [email protected]


SeeAlsokrige.tilesov ,krige

Examples

# compare result and system.time() with example(krige.pvm)

system.time(

maas.ktilesov.pvm <-

krige.pvm.tiles(maas.point, "zinc", maas.vmod,



nxper=2, nyper=2, rsearch = 600,

extrap = F, border=maas.bank)

)


plot(maas.ktilesov.pvm,show.snb=TRUE,snb.i=50,snb.j=50)

List of Figures

2.1 Spherical Model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

2.2 Exponential Model . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

2.3 Gaussian Model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

3.1 Overlapping Search Neighbourhoods. . . . . . . . . . . . . . . . . . 29

3.2 Grid Points in Tiles. . . . . . . . . . . . . . . . . . . . . . . . . . . 31

3.3 Grid Points in Overlapping Tiles. . . . . . . . . . . . . . . . . . . . 35

3.4 Computational Burden of krige.grid, krige.tiles and krige.tilesov. . . 39

3.5 Plot of maas.kgrid, where nx=ny=100. . . . . . . . . . . . . . . . . 40

3.6 Plot of maas.ktiles, where nx=ny=100, itx=ity=10. . . . . . . . . . . 41

3.7 Plot of maas.ktilesov, where nx=ny=100, itx=ity=10, nxper=nyper=2. 42

5.1 PVM Computational Model . . . . . . . . . . . . . . . . . . . . . . 49

5.2 Parallelizing the Tiles. . . . . . . . . . . . . . . . . . . . . . . . . . 51

5.3 Computational Effort with different sets of computers. . . . . . . . . 56

5.4 XPVM: Communication between Server and Clients. . . . . . . . . 58

76

Bibliography

[1] E. Anderson and Others.LAPACK Users’ Guide. SIAM Philadelphia, third

edition, 1999.

[2] Ricardo A.Olea.Geostatistics for Engineers and Earth Scientists. Kluwer Aca-

demic Publishers, 1999.

[3] William S. Cleveland.Visualizing Data. AT & T Bell Laboratories, 1993.

[4] Noel A.C. Cressie.Statistics for Spatial Data. Wiley, 1993.

[5] Peter J. Diggle.Statistical Analysis of Spatial Point Patterns. Academic Press,

1983.

[6] Al Geist and Others.PVM: Parallel Virtual Machine - A Users Guide and Tu-

torial for Networked Parallel Computing. The MIT Press, Cambridge, Mas-

sachusetts, second edition, 1995.

[7] N.L. Hjort and H. Omre. Topics in spatial statistics.Scandinavian Journal of

Statistics, 1994.

[8] Ross Ihaka and Robert Gentleman. R: A language for data analysis and graphics.

Journal of Computational and Graphical Statistics, 1996.

[9] Edward H. Isaaks and R. Mohan Srivastava.Applied Geostatistics. Oxford Uni-

versity Press, 1989.

77

BIBLIOGRAPHY 78

[10] Piere Delfiner J.-P. Chiles and Paul Delfiner.Geostatistics: Modeling Spatial

Uncertainty. Wiley, 1999.

[11] Andreas Krause.Einfuhrung in S und S-Plus. Springer, 1997.

[12] Georges Matheron.Geostatistical Case Studies. Reidel, 1987.

[13] Hannes Tschofenig. Anwendung der parallelen, virtuellen Maschine (PVM) fur

eine Geostatistikanwendung. Universitat Klagenfurt, 2001.

[14] William N. Venables and Brian D. Ripley.Modern Applied Statistics with S-Plus.

Springer, 1994.

[15] Hans Wackernagel.Multivariate Geostatistics - An Introduction with Applica-

tions. Springer, 1995.

Index

anisotropy,16

geometric,16

zonal,16

average

weighted,10, 17

bias,21

block kriging,20

correlation function,12

covariance function,12

covariance matrix,19, 23

deterministic,21

distributed computing,39

drift, 21

estimation error,17

estimation variance,10, 18, 22

exact interpolator,20

grid, 24

overlapping tiles,29

tiled, 26

independent

linearly,23

intrinsically stationary,9, 19

kriging estimator,18

Kriging functions

krige.grid,53

krige.tiles,56

krige.tilesov,60

Lagrange parameter,17, 19, 23

linear equation system,25

LU factorization,25

mean

arithmetic,17

weighted,17

MPP,39

negative definite

conditionally,12

non-stationary,21

nugget-effect,11, 15

ordinary kriging,18, 20

prediction,19

system,19, 20

variance,19, 20

79

INDEX 80

parallel tile kriging,44

positive semidefinite,12, 19

PVM

Introduction,38

library, 40

master-slave,41

pvmd3,40

system,40

task identifier,41

tasks,40

PVM Kriging functions

krige.tilesov,65

random field,8

random function,8

random variable,8

range,16

regionalized data,8

regionalized variable,8

residual,10, 21, 23

search neighbourhood,25

second-order stationary,9, 18

semivariogram,9, 19

cloud,11

experimental,11

exponential model,14

gaussian model,14

linear model,15

models,13

power model,15

spherical model,13

theoretical,12

sill, 15

simple kriging,10

equation system,10

Spatial Data,7

strictly stationary,8

universal kriging,21

system,23

universal pvm kriging,42

variogram,9

virtual machine,39

Parallel Computing in Spatial StatisticsThis means if we translate a point conﬁguration in a given...

Documents

Transcript of Parallel Computing in Spatial StatisticsThis means if we translate a point conﬁguration in a given...