Efficient Approximation
of Random Fields
for Numerical Applications
Helmut Harbrecht, Michael Peters,
Markus Siebenmorgen
Institute of Mathematics Preprint No. 2014-01
University of Basel January, 2014
Rheinsprung 21
CH - 4051 Basel
Switzerland www.math.unibas.ch
NUMERICAL LINEAR ALGEBRA WITH APPLICATIONS
Numer. Linear Algebra Appl. 2014; 00:1–20
Published online in Wiley InterScience (www.interscience.wiley.com). DOI: 10.1002/nla
Efficient approximation of random fields for numerical
applications
†
Helmut Harbrecht, Michael Peters
∗
, Markus Siebenmorgen
Mathematisches Institut, Rheinsprung 21, 4051 Basel, Switzerland
SUMMARY
This article is dedicated to the rapid computation of separable expansions for the approximation of random
fields. We consider approaches based on techniques from the approximation of non-local operators on the
one hand and based on the pivoted Cholesky decomposition on the other hand. Especially, we provide an
a-posteriori error estimate for the pivoted Cholesky decomposition in terms of the trace norm. Numerical
examples are provided to validate and quantify the presented methods. Copyright
c
2014 John Wiley &
Sons, Ltd.
Received . . .
1. INTRODUCTION
In this article, we present and compare two different approaches for the approximation of random
fields in L
2
P
Ω, H
p
(D)
for a spatial domain D ∈ R
d
and a probability space (Ω, F, P). Stochastic
fields appear for example in the modeling of diffusion problems with random data, see e.g. [1],
and in machine learning, see e.g. [2]. To make a stochastic field a(x, ω) feasible for numerical
computations in a stochastic Galerkin or stochastic collocation method, see e.g. [1, 3, 4, 5, 6, 7] and
the references therein, one has to separate the spatial variable x and the stochastic variable ω. Since
L
2
P
Ω, H
p
(D)
∼
=
L
2
P
(Ω) ⊗H
p
(D), see e.g. [8], this task can be accomplished by computing a basis
representation of a in L
2
P
(Ω) ⊗ H
p
(D). A very common approach to obtain such a representation is
the Karhunen-Lo
`
eve expansion, cf. [1, 9], which can be regarded as linear operator analogue of the
singular value decomposition of matrices.
The main task in the computation of a Karhunen-Lo
`
eve expansion is the solution of a symmetric
and positive semidefinite eigen-problem. In this context, approaches to efficiently compute the
Karhunen-Lo
`
eve expansion have been made by means of the Fast Multipole Method (FMM) based
on interpolation, cf. [10], in [11] and with the aid of H-matrices, cf. [12], in [13]. The idea in
these articles is to provide a data-sparse representation of the covariance operator which is then
used to solve the related eigen-problem numerically by a Krylov subspace method, cf. [14]. Of
course, another algorithm for the efficient approximation of non-local operators, like the Adaptive
Cross Approximation (ACA), cf. [15, 16], or the Wavelet Galerkin Scheme (WGS), cf. [17, 18],
can be considered as well for the representation of the covariance operator. Nevertheless, the major
drawback of these approaches is that the number of eigenvalues to be computed has to be known in
advance which might be a strong assumption in practice.
∗
Correspondence to: E-mail: michael.peters@unibas.ch
†
This research has been supported by the Swiss National Science Foundation (SNSF) through the project “Rapid Solution
of Boundary Value Problems on Stochastic Domains”.
Copyright
c
2014 John Wiley & Sons, Ltd.
Prepared using nlaauth.cls [Version: 2010/05/13 v2.00]
2
To overcome this obstruction, we present here an alternative approach based on the Pivoted
Cholesky Decomposition (PCD). The PCD can be interpreted as a single-block ACA with applicable
total pivoting, cf. [19]. Hence, only the main diagonal of the discretized operator has to be
precomputed, which can be performed in essentially, i.e. up to possible poly-logarithmic terms,
linear complexity, if the quadrature proposed in [20] is applied to discretize the operator. Then, in
each step of the algorithm, the quality of the approximation with respect to the stochastic field is
controllable by means of the trace norm. If the desired accuracy is achieved, the algorithm stops with
an M-term approximation to the operator. If M is substantially smaller than the dimension of the
ansatz space, we end up with a remarkable computational speed-up. The related Karhunen-Lo
`
eve
expansion might then be computed in a post-processing step. Notice that then the PCD yields a full
but relatively small eigen-problem if the operator under consideration exhibits a certain smoothness.
This eigen-problem might be solved numerically exact by e.g. the QR-method, cf. [21].
Now the following question arises: which approach is more efficient? We will try to answer this
question numerically by comparing the PCD with methods lend from the approximation of non-
local operators. We employ here ACA for the data-sparse approximation of the covariance operator
which results in a fast matrix-vector product. Thus, a Krylov subspace method – we use the Implicit
Restarted Arnoldi Method (IRAM), cf. [22, 23, 24] – is feasible to compute the desired eigenvalues
of largest magnitude.
Finally, we would like to emphasize that, although we focus here on the application to random
fields, the presented methods are also applicable in the more general case of approximating bi-
variate functions in L
2
(D
1
) ⊗ L
2
(D
2
) for two domains D
1
⊂ R
d
1
and D
2
⊂ R
d
2
.
The rest of this article is structured as follows. Section 2 is devoted to the Karhunen-Lo
`
eve
expansion. Especially, we discuss here the related error estimates including discretization and
truncation error. To that end, it is crucial to have bounds for the decay of the covariance operator’s
eigenvalues. These bounds are considered here, too. In Section 3, we provide the theoretical
background for the pivoted Cholesky decomposition. Moreover, we establish error estimates for
the approximation of random fields in terms of the trace norm. These estimates are essential for the
a-posteriori control of the approximation error. Section 4 introduces a special class of covariance
functions based on the Mat
´
ern kernel functions. We choose this class of covariance functions for our
numerical tests, since we a-priori know the decay rate of the respective eigenvalues. In particular, we
are also able to analytically compute the eigenfunctions and eigenvalues in the case of the unit sphere
S
2
. Thus, these kernels provide an excellent benchmark to compare both approaches. Section 5 is
dedicated to testing the numerical performance of the methods under consideration. We will solve
the eigenvalue problem for covariance operators related to some of the Mat
´
ern kernels from Section
4 on different geometries. Finally, we sum up the results presented within this article in Section 6.
In the following, in order to avoid the repeated use of generic but unspecified constants, by C . D
we mean that C can be bounded by a multiple of D, independently of parameters which C and D
may depend on. Obviously, C & D is defined as D . C, and C h D as C . D and C & D.
2. THE KARHUNEN-LO
`
EVE EXPANSION
Let (Ω, F, P) be a probability space with σ-field F ⊂ 2
Ω
and a complete probability measure P, i.e.
for all A ⊂ B and B ∈ F with P[B] = 0 it follows A ∈ F. Furthermore, let D ⊂ R
d
for d = 2, 3 be
a sufficiently smooth and bounded domain.
For p ≥ 0, the Lebesgue-Bochner space L
2
P
Ω; H
p
(D)
consists of all maps
a: Ω → H
p
(D)
that satisfy
kvk
L
2
P
(Ω;H
p
(D))
:
=
Z
Ω
kv(·, ω)k
2
H
p
(D)
dP(ω)
1/2
< ∞. (1)
In the following, it will be convenient to identify L
2
P
Ω; H
p
(D)
according to
L
2
P
Ω; H
p
(D)
∼
=
H
p
(D) ⊗ L
2
P
(Ω).
Copyright
c
2014 John Wiley & Sons, Ltd. Numer. Linear Algebra Appl. (2014)
Prepared using nlaauth.cls DOI: 10.1002/nla
3
For further details on Lebesgue-Bochner spaces see e.g. [25].
For the approximation of spatial functions in L
2
(D), we will consider piecewise continuous finite
elements. Therefore, we introduce a family of quasi-uniform triangulations T
h
for D with mesh
width h and define the spaces
V
s
h
:
= {v
h
: D → R : v|
T
is a polynomial of order s for all T ∈ T
h
} ⊂ L
2
(D). (2)
Then, given a function v ∈ H
p
(D) with 0 ≤ p ≤ s, we have due to the Bramble-Hilbert lemma the
approximation estimate
inf
v
h
∈V
s
h
kv − v
h
k
L
2
(D)
. h
p
kvk
H
p
(D)
(3)
uniformly in h, see e.g. [26, 27].
A very common representation of random fields for numerical purposes is given by the Karhunen-
Lo
`
eve expansion. In order to ensure that L
2
P
(Ω) is separable, we have to assume that Ω is a separable
set.
Definition 2.1. Let a ∈ H
p
(D) ⊗ L
2
P
(Ω) for some p ≥ 0 be a random field. The expansion
a(x, ω) = a(x) +
∞
X
m=1
σ
m
ϕ
m
(x)X
m
(ω) (4)
with σ
1
≥ σ
2
≥ ··· ≥ 0, (X
m
, X
n
)
L
2
P
(Ω)
= δ
m,n
and (ϕ
m
, ϕ
n
)
L
2
(D)
= δ
m,n
is called Karhunen-
Lo
`
eve expansion with respect to a. Here, a(x) denotes the mean of a with respect to the stochastic
variable, i.e.
a(x) =
Z
Ω
a(x, ω) dP(ω).
The Karhunen-Lo
`
eve expansion can be regarded as the continuous analogue to the singular value
decomposition of matrices. Especially, it holds σ
m
=
√
λ
m
, where {(λ
m
, ϕ
m
)}
m
are the eigen-pairs
(in decreasing order) of the covariance operator
(Cu)(x)
:
=
Z
D
k(x, y)u(y) dy, (5)
given via the correlation kernel
k(x, y)
:
=
Z
Ω
a(x, ω) − a(x)
a(y, ω) − a(y)
dP(ω).
Notice that, for a ∈ H
p
(D) ⊗ L
2
P
(Ω), it holds k ∈ H
p,p
mix
(D × D)
:
= H
p
(D) ⊗ H
p
(D). Addition-
ally, the random variables {X
m
}
m
are given by
X
m
(ω) =
1
σ
m
Z
D
a(x, ω) − a(x)
ϕ
m
(x) dx.
In the following, we will also make use of the Hilbert-Schmidt operator associated with the
centered random field, i.e. S: L
2
P
(Ω) → H
p
(D) with
(Su)(x) =
Z
Ω
a(x, ω) − a(x)
u(ω) dP(ω) for u ∈ L
2
P
(Ω)
and its adjoint S
?
:
˜
H
−p
(D) → L
2
P
(Ω) with
(S
?
u)(ω) =
Z
D
a(y, ω) − a(y)
u(y) dy for u ∈
˜
H
−p
(D).
Then, we especially find that SS
?
:
˜
H
−p
(D) → H
p
(D) is given by
(SS
?
u)(x) =
Z
Ω
a(x, ω) − a(x)
Z
D
a(y, ω) − a(y)
u(y) dy dP(ω) = (Cu)(x),
Copyright
c
2014 John Wiley & Sons, Ltd. Numer. Linear Algebra Appl. (2014)
Prepared using nlaauth.cls DOI: 10.1002/nla
4
where the last identity holds due to Fubini’s theorem.
For numerical issues, one has to truncate the Karhunen-Lo
`
eve expansion appropriately. Here, the
truncation error depends on the decay of the covariance operator’s eigenvalues. More precisely, for
the decay of the eigenvalues, we have along the lines of [28] the following theorem.
Theorem 2.2. Let a ∈ H
p
(D) ⊗ L
2
P
(Ω). Then, the eigenvalues of the covariance operator
C:
˜
H
−p
(D) → H
p
(D) decay like
λ
m
. m
−2p/d
for m → ∞.
Proof
We consider the operator S
?
S: L
2
P
(D) → L
2
P
(D). Let (λ
m
, X
m
) be an eigen-pair of S
?
S. On the
one hand, we have
C(SX
m
) = λ
m
(SX
m
).
On the other hand, it holds for λ
m
> 0 that
(SX
m
, SX
m
)
L
2
(D)
= (S
?
SX
m
, X
m
)
L
2
P
(Ω)
= λ
m
> 0. (6)
Thus, we conclude that (λ
m
, SX
m
/σ
m
) is an eigen-pair of C. The proof is now based on an
approximation argument.
We consider the approximation space V
dpe
h
⊂ L
2
(D), cf. (2). Let dim(V
dpe
h
) = N. Notice that
h h N
−1/d
, where the constant depends on the polynomial degree dpe. Furthermore, we define the
L
2
(D)-orthogonal projection Q
N
: L
2
(D) → V
dpe
h
. Then, due to the Bramble-Hilbert lemma, we
conclude the estimate
k(I − Q
N
)vk
L
2
(D)
. N
−p/d
kvk
H
p
(D)
for u ∈ H
p
(D).
The min-max principle of Courant-Fisher implies now for arbitrary subspaces V
m
⊂ L
2
P
(Ω) with
dim(V
m
) ≤ m that
λ
m+1
= min
V
m
max
v∈V
⊥
m
,kvk
L
2
P
(Ω)
=1
(S
?
Sv, v)
L
2
P
(Ω)
= min
V
m
max
v∈V
⊥
m
,kvk
L
2
P
(Ω)
=1
(Sv, Sv)
L
2
(D)
.
For the choice V
N
= img(S
?
Q
N
S), the orthogonality of the projection Q
N
yields
λ
N+1
≤ max
v⊥img(S
?
Q
N
S),kvk
L
2
P
(Ω)
=1
(Sv, Sv)
L
2
(D)
= max
v⊥img(S
?
Q
N
S),kvk
L
2
P
(Ω)
=1
Sv, (I − Q
N
)Sv
L
2
(D)
= max
v⊥img(S
?
Q
N
S),kvk
L
2
P
(Ω)
=1
(I − Q
N
)Sv, (I − Q
N
)Sv
L
2
(D)
≤ sup
kvk
L
2
P
(Ω)
=1
(I − Q
N
)Sv, (I − Q
N
)Sv
L
2
(D)
≤ sup
kvk
L
2
P
(Ω)
=1
k(I − Q
N
)Svk
2
L
2
(D)
. N
−2p/d
sup
kvk
L
2
P
(Ω)
=1
kSvk
H
p
(D)
.
This estimate together with the continuity of S, see e.g. [28], yields the assertion.
Now, in accordance with [28], an estimation of the Karhunen-Lo
`
eve expansion’s truncation error
is provided by the following theorem.
Copyright
c
2014 John Wiley & Sons, Ltd. Numer. Linear Algebra Appl. (2014)
Prepared using nlaauth.cls DOI: 10.1002/nla
