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
CHAPTER 2. SPATIAL STATISTICS 8
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)
CHAPTER 2. SPATIAL STATISTICS 9
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
CHAPTER 2. SPATIAL STATISTICS 10
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.
CHAPTER 2. SPATIAL STATISTICS 11
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)
CHAPTER 2. SPATIAL STATISTICS 12
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:
CHAPTER 2. SPATIAL STATISTICS 13
• 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)
CHAPTER 2. SPATIAL STATISTICS 14
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
CHAPTER 2. SPATIAL STATISTICS 15
• 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
CHAPTER 2. SPATIAL STATISTICS 16
• 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
CHAPTER 2. SPATIAL STATISTICS 17
• 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]).
CHAPTER 2. SPATIAL STATISTICS 18
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)
CHAPTER 2. SPATIAL STATISTICS 19
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)
CHAPTER 2. SPATIAL STATISTICS 20
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-
CHAPTER 2. SPATIAL STATISTICS 21
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
CHAPTER 2. SPATIAL STATISTICS 22
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
CHAPTER 2. SPATIAL STATISTICS 23
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))
CHAPTER 2. SPATIAL STATISTICS 24
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)
CHAPTER 2. SPATIAL STATISTICS 25
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.
CHAPTER 3. KRIGING ON TILED GRIDS 28
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
CHAPTER 3. KRIGING ON TILED GRIDS 29
S1 S1_2
x1
. x2
.
x1
.x2
.
. .
. .
. .
. .
. .
. .
... .
. ...
. .
S2
Figure 3.1: Overlapping Search Neighbourhoods
CHAPTER 3. KRIGING ON TILED GRIDS 30
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:
CHAPTER 3. KRIGING ON TILED GRIDS 31
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
Figure 3.2: Grid Points in Tiles
CHAPTER 3. KRIGING ON TILED GRIDS 32
• 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)
CHAPTER 3. KRIGING ON TILED GRIDS 33
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
CHAPTER 3. KRIGING ON TILED GRIDS 34
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.
CHAPTER 3. KRIGING ON TILED GRIDS 35
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
Figure 3.3: Grid Points in Overlapping Tiles
CHAPTER 3. KRIGING ON TILED GRIDS 36
• 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
CHAPTER 3. KRIGING ON TILED GRIDS 37
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
CHAPTER 3. KRIGING ON TILED GRIDS 38
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.
CHAPTER 3. KRIGING ON TILED GRIDS 39
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
krige.gridkrige.tileskrige.tilesov
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
krige.gridkrige.tileskrige.tilesov
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
krige.gridkrige.tileskrige.tilesov
Figure 3.4: Computational Burden of krige.grid, krige.tiles and krige.tilesov
CHAPTER 3. KRIGING ON TILED GRIDS 40
178500 179000 179500 180000 180500 181000
3300
0033
1000
3320
0033
3000
Figure 3.5: Plot of maas.kgrid, where nx=ny=100
CHAPTER 3. KRIGING ON TILED GRIDS 41
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
CHAPTER 3. KRIGING ON TILED GRIDS 42
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
CHAPTER 4. PVM (PARALLEL VIRTUAL MACHINE) 45
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.
CHAPTER 4. PVM (PARALLEL VIRTUAL MACHINE) 46
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.
CHAPTER 5. UNIVERSAL PVM KRIGING 49
Rlibrary(pvmkrige)
krige_server
krige_client krige_client krige_client
P MV
Figure 5.1: PVM Computational Model
CHAPTER 5. UNIVERSAL PVM KRIGING 50
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
CHAPTER 5. UNIVERSAL PVM KRIGING 51
. . .
. . .
. . . . . . . .
. . . . .
. . .
. . .
. . .
. . .
. . . . . . . .
. . . . .
. . .
. . .
. . .
. . .
. . .
. . .
. . .
. . .
. . .
. . .
. . .
R SERVER
CLIENT A
CLIENT B
CLIENT C
Figure 5.2: Parallelizing the Tiles
CHAPTER 5. UNIVERSAL PVM KRIGING 52
• 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
CHAPTER 5. UNIVERSAL PVM KRIGING 53
• 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
CHAPTER 5. UNIVERSAL PVM KRIGING 54
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.
CHAPTER 5. UNIVERSAL PVM KRIGING 55
hostname architecture CPU speed memory OS
alpha1-stat Alpha EV4 200 MHz 128 MB Tru64
alpha2-stat Alpha EV4 150 MHz 80 MB Tru64
alpha3-stat Alpha EV4 150 MHz 48 MB Tru64
alpha4-stat Alpha EV6 667 MHz 1 GB Tru64
alpha5-stat Alpha EV5 500 MHz 512 MB 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
CHAPTER 5. UNIVERSAL PVM KRIGING 56
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
tile size:10, overlapping points:6number of grid points
time
in s
econ
ds
8 computers5 computers2 computers
50 100 150 200 250 300
050
100
150
tile size:30, overlapping points:8number of grid points
time
in s
econ
ds
8 computers5 computers2 computers
50 100 150 200 250 300
050
100
150
tile size:30, overlapping points:16number of grid points
time
in s
econ
ds
8 computers5 computers2 computers
Figure 5.3: Computational Effort with different sets of computers
CHAPTER 5. UNIVERSAL PVM KRIGING 57
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
CHAPTER 5. UNIVERSAL PVM KRIGING 58
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
APPENDIX A. FUNCTIONS 62
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.
APPENDIX A. FUNCTIONS 63
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
APPENDIX A. FUNCTIONS 64
# 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,
xne=NULL, yne=NULL, dx=NULL, dy=NULL, nx=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
point.obj spatial data, class ”point ”
at name of variable of interest inpoint.obj
APPENDIX A. FUNCTIONS 65
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
itx number of grid points per tile in x direction
ity number of grid points per tile in y direction
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
APPENDIX A. FUNCTIONS 66
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
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).
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).
APPENDIX A. FUNCTIONS 67
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]
Referencesftp://ftp-stat.uni-klu.ac.at/pub/R/contrib
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),
APPENDIX A. FUNCTIONS 68
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,
xne=NULL, yne=NULL, dx=NULL, dy=NULL, nx=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
point.obj spatial data, class ”point ”
at name of variable of interest inpoint.obj
var.mod.obj variogram model, generated byfit.variogram
APPENDIX A. FUNCTIONS 69
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
itx number of grid points per tile in x direction
ity number of grid points per tile in y direction
nxper number of overlapping grid points per tile in x direction
nyper number of overlapping grid points per tile in y direction
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
APPENDIX A. FUNCTIONS 70
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
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.
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
APPENDIX A. FUNCTIONS 71
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.tiles .
Note
In contrast tokrige.tiles the border effects between tiles are eliminated.
Author N. Samonig, [email protected]
Referencesftp://ftp-stat.uni-klu.ac.at/pub/R/contrib
SeeAlsokrige.cell ,krige
Examples
# compare result and system.time() with example(krige.tiles)
system.time(
maas.ktilesov <-
krige.tilesov(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,
nxper=2, nyper=2, rsearch = 600,
extrap = F, border=maas.bank)
)
# show a picture with tile containing point(50,50) marked :
plot(maas.ktilesov,show.snb=TRUE,snb.i=50,snb.j=50)
APPENDIX A. FUNCTIONS 72
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,
xne=NULL, yne=NULL, dx=NULL, dy=NULL, nx=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
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
APPENDIX A. FUNCTIONS 73
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
itx number of grid points per tile in x direction
ity number of grid points per tile in y direction
nxper number of overlapping grid points per tile in x direction
nyper number of overlapping grid points per tile in y direction
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
APPENDIX A. FUNCTIONS 74
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
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 .
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).
APPENDIX A. FUNCTIONS 75
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]
Referencesftp://ftp-stat.uni-klu.ac.at/pub/R/contrib
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,
min(maas.bank$x), min(maas.bank$y), max(maas.bank$x),
max(maas.bank$y), nx=100,ny=100, itx=10, ity=10,
nxper=2, nyper=2, rsearch = 600,
extrap = F, border=maas.bank)
)
# show a picture with tile containing point(50,50) marked :
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