Karhunen–Loève expansion for mean-square integrable random fields

Let $H$ be a separable Hilbert space. A linear operator $C\colon H\mapsto H$ is the covariance of two $H$-valued random variables $X$ and $Y$ if \begin{align} (C\phi, \psi)_H = Cov((\phi,X),(\psi,Y)),\quad \forall \phi,\psi\in H. \end{align} X and Y are said to be uncorrelated if $C$ is the zero operator.

Let $c\in L^2(D\times D)$ be a correlation function. Then, a covariance operator $C$ is defined by \begin{align} (C\phi)(x) = \int\limits_{D}c(x,y)\phi(y)\mathrm{d}y,\quad x\in D. \end{align}


For $D\subset \mathbb{R}^d$ and $\kappa\in L^2(\Omega, L^2(D))$ consider for $\phi, \psi\in L^2(D)$ the covariance operator $C\colon L^2(D)\to L^2(D)$ \begin{align} (C\phi, \psi)_{L^2(D)} = Cov((\phi, u)_{L^2(D)}, (\psi, u)_{L^2(D)}). \end{align} By definition of the real valued covariance, one obtains using Fubinis Theorem \begin{align} Cov((\phi, u)_{L^2(D)}, (\psi, u)_{L^2(D)}) = \int_D\int_D \mathbb{E}[(u(x) - \mathbb{E}[u(x)])(u(y) - \mathbb{E}[u(y)])]\phi(x)\psi(y)\mathrm{d}x\mathrm{d}y, \end{align} which implies, $Cov(u(x), u(y))$ is the correlation function (kernel) of the covariance operator $C$, i.e. \begin{align} (C\phi)(x) = \int_D Cov(u(x), u(y)) \phi(y)\mathrm{d}y. \end{align}

Let $C$ be positive, self-adjoint and compact. Suppose that $\kappa\in L^2(\Omega, L^2(D))$. Then \begin{align} \kappa(x,\omega) = \mu(x) + \sum\limits_{m=1}^\infty\sqrt{\lambda_m}\phi_m(x)\xi_j(\omega) \end{align} with convergence in $L^2(\Omega,L^2(D))$. Here $\{(\lambda_m,\phi_m)\}$ are eigenfunctions of the operator $C$. The random variables $\xi_j$ have mean zero, unit variance and are pairwise uncorrelated and \begin{align} \xi_j(\omega) := \frac{1}{\sqrt{\lambda_j}}(\kappa(\omega)-\mu, \phi_j)_{L^2(D)}. \end{align} If the field $\kappa$ is Gaussian, then $\xi\sim \mathcal{N}(0,1)$ iid.

A truncated KLE and thus field representation of order $M< \infty$ reads \begin{align} \kappa_M(x,\omega):=\mu(x) + \sum\limits_{m=1}^M\sqrt{\lambda_m}\phi_m(x)\xi_j(\omega) \end{align}

The truncation error will depend on the decay of eigenvalues.

remark I

The truncation level $M$ corresponds with a random dimension.

remark II

More regularity of $c$ may imply better convergence situation of the KLE. E.g. if $c\in\mathcal{C}(\overline{D}\times\overline{D})$, then eigenfunctions are continous and \begin{align} \operatorname{sup}\limits_{x\in D} \mathbb{E}[ (\kappa(x) - \kappa_M(x))^2]\rightarrow 0\text{ as } M \to \infty. \end{align}

Matérn covariance function

Let $r = |x-y|$, then the Matérn covariance function is given as \begin{align} c(x,y) = c(r) = \sigma^2\frac{2^{1-\nu}}{\Gamma(\nu)}\left(\sqrt{2\nu}\frac{r}{\rho}\right)^{\nu}K_\nu\left(\sqrt{2\nu}\frac{r}{\rho}\right) \end{align} with gamma function $\Gamma$, $K_\nu$ modified Bessel function of second kind and $\rho,\nu\geq 0$. The parameter $\rho$ accounts for the correlation length and $\nu$ represents the smoothness of the covariance function.

special case I $\nu =1/2$

The exponential covariance function is given as \begin{align} c(r) = \sigma^2 \exp\left(-\frac{r}{\rho}\right). \end{align}

special case II $\nu\to\infty$

The squared exponential covariance function \begin{align} c(r) = \sigma^2\exp\left(-\frac{r^2}{2\rho^2}\right). \end{align}


Solve the eigenvalue problem $C\phi = \lambda \phi$ in $H=L^2(D)$ with an underlying exponential and squared exponential covariance function $c$. Plot and discuss the decay of eigenvalues with varying values of $\rho > 0$. An implementation of the exponential covariance function and the solver is given below.

a simple approach:

Let $V_h\subset H^1(D)\subset L^2(D)$ be a $H^1(D)$ conform FE space with $V_h = \operatorname{span}\{\phi_k\}$. Using the FE approximation $c(x,y) \approx \sum\limits_{ij}c_{ij} \phi_i(x)\phi_j(y)$ in $L^2(D\times D)$. Then, a discretisation yield a generalised eigenvalue problem of the form. \begin{align} M^T[C]M \alpha = \lambda M \alpha, \end{align} with $[C]_{ij}=c_{ij}$. Note that this is a symmetric generalised EVP in opposite to the equivalent EVP $[C]M\alpha = \lambda\alpha$.


Test the slow validation subversion. Increase the refinement level $N$ and discuss the performance (with vectorization only!).

In [1]:
from __future__ import division
from dolfin import *
import scipy.sparse.linalg as spla
import numpy as np

Example implementation of the exponential covariance function.

In [2]:
def cov_exp(r, rho, sigma2=1.0):
    return sigma2 * np.exp(-1/rho * r)
In [3]:
def solve_covariance_EVP(cov, N, k, degree=1):
    def setup_FEM(N):
        mesh = UnitSquareMesh(N,N)
        V = FunctionSpace(mesh, 'CG', degree)
        u = TrialFunction(V)
        v = TestFunction(V)
        return mesh, V, u, v
    # construct FEM space
    mesh, V, u, v = setup_FEM(N)

    # dof to vertex map
    dof2vert = dof_to_vertex_map(V)
    # coords will be used for interpolation of covariance kernel
    coords = mesh.coordinates()
    # but we need degree of freedom ordering of coordinates
    coords = coords[dof2vert]
    # assemble mass matrix and convert to scipy
    M = assemble(u*v*dx)
    M = M.array()

    # evaluate covariance matrix
    L = coords.shape[0]
    if True: # vectorised
            c0 = np.repeat(coords, L, axis=0)
            c1 = np.tile(coords, [L,1])
            r = np.abs(np.linalg.norm(c0-c1, axis=1))
            C = cov(r)
            #C = cov(c0-c1)
            C.shape = [L,L]
    else:   # slow validation
        C = np.zeros([L,L])
        for i in range(L):
            for j in range(L):
                if j <= i:
                    v = cov(np.linalg.norm(coords[i]-coords[j]))
                    C[i,j] = v
                    C[j,i] = v
    # solve eigenvalue problem
    A = np.dot(M, np.dot(C, M))
    w, v = spla.eigsh(A, k, M)
    # return eigenpairs
    print "===== EVP size =====", A.shape, w.shape, v.shape
    #v = np.array([z[dof2vert] for z in v.T])
    return w, v, V

With the method we can now approximate the eigenvalues.

In [4]:
 lambdas, phis, V  = solve_covariance_EVP(lambda r : cov_exp(r, rho=0.1, sigma2=1.0), 
                                            N = 30, k = 30, degree = 1)
 lambdas_2, phis_2, V_2  = solve_covariance_EVP(lambda r : cov_exp(r, rho=1., sigma2=1.0), 
                                            N = 30, k = 30, degree = 1)
===== EVP size ===== (961, 961) (30,) (961, 30)
===== EVP size ===== (961, 961) (30,) (961, 30)

And plot them...

In [5]:
import matplotlib.pyplot as plt
% matplotlib inline
plt.suptitle("Exponential covariance function")
plt.semilogy(range(len(lambdas)), sorted(lambdas, reverse=True))
plt.semilogy(range(len(lambdas_2)), sorted(lambdas_2, reverse=True))
plt.tight_layout(rect=[0, 0.03, 1, 0.95])

Field representation on finite element function space

Affine coefficient representation

As an example of an artificial random field, given e.g. by a KL formulation, we take our usual example field in affine representation of the finite i.i.d. random variables $\boldsymbol \xi = (\xi_1, \ldots,\xi_M)$ $$ \tag 1 \kappa(x, \boldsymbol{\xi}(\omega)) = \langle \kappa \rangle + \gamma \sum_{k=1}^M k^{-\sigma}\cos\left(2\pi\left\lceil \frac{k+1}{2}\right\rceil x_1\right)\cos\left(2\pi\left\lfloor \frac{k+1}{2} \right\rfloor x_2\right) \xi_k(\omega) $$ where $\sigma$ denotes a decay parameter that controls the roughness of a realisation, $\gamma$ is a scaling factor and $\langle \kappa \rangle$ is the mean value of the field.

To ensure uniform boundedness $0 < \check{\kappa} \leq \kappa(x, \boldsymbol\xi(\omega))$ we need to adjust the involved parameter. In the uniform case of $\boldsymbol\xi\sim\mathcal{U}([-1, 1]^M)$, one can easily see the dependence on the choice of $\gamma > 0$, $\sigma > 1$ and the mean value $\langle \kappa \rangle$.


Write a python class/object that represents the diffusion coefficient $\kappa\colon D\times\Omega \to \mathbb{R}$ as in (1). A coefficient must implement the method realisation(y, V).

realisation takes a sample value $y=\boldsymbol\xi(\omega)$ and a FEniCS FunctionSpace $V$ as input and returns the interpolated FEniCS function of the Expression above (1).

Test your implementation using $\boldsymbol\xi$ as uniform random variable on $[-1, 1]^M$ and different values for $\langle\kappa\rangle, \gamma$ and $\sigma$. E.g. $\langle \kappa\rangle = 1.1, \sigma=2, \gamma=1$.

Plot some realisations.

Hint: There exists functions np.ceil and np.floor. Maybe, the parametrized version of FEniCS Expression can be useful.

Log-Normal Diffusion coefficient

In practice, especially when considering the ground-water flow motivation, diffusion fields are often modeled by a more involved structure. The so-called log-normal case is defined by setting $\log \kappa$ as a Gaussian field. We can model this by simply taking the field in (1) and assuming that $\boldsymbol\xi$ is a standard Gaussian random vector, due to the linearity of the affine representation, $\log \kappa$ is again Gaussian.

To ensure positivity, one usually takes $\exp \kappa$ as diffusion coefficient.


Adapt your implementation of the last exercise such that the function $\kappa\colon D\times \Omega\to \mathbb{R}$ can represent a log-normal diffusion field.

Plot some realisations.

Compute statistics of random fields

Consider a random variable $X\colon \Omega\to E$ defined on some probability space $(\Omega, \Sigma, \mathbb{P})$ with values in some vector space $E\subseteq \mathbb{R}^d$. The $E$-valued mean of $X$ can be approximated by the most simple quadrature rule:

Given a set of samples $(x_1, \ldots, x_N)$ of $X$, compute $$ \mathbb{E}[X] = \int_\Omega x \;\mathbb{P}(dx) \approx \frac{1}{N}\sum_{i=1}^N x_i. $$

The same holds for the variance of $X$. Here, we need to compute the second moment of $X$, since $$ \mathbb{V}(X) = \mathbb{E}[X^2] - \mathbb{E}[X]^2 $$ and obviously $$ \mathbb{E}[X^2] = \int_\Omega x^2 \;\mathbb{P}(dx) \approx \frac{1}{N}\sum_{i=1}^N x_i^2. $$


Implement a python method that approximates the mean and variance of a random variable

In [6]:
from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
from dolfin import * # UnitSquareMesh, FunctionSpace, plot, Function, Expression, interpolate
from functools import partial

% matplotlib inline
In [7]:
class CoefField:
    """ artificial M-term KL
    def __init__(self, mean, expfield, M, scale=1.0, decay=2.0):
        # type: (float, bool, int, float, float) -> None
         initialise field with given mean.
         field is log-uniform for expfield=True.
         length M of expansion can optionally be fixed in advance.
         usually it is determined by the length of y when evaluating realisations.
        :param mean: mean value of the field
        :param expfield: Switch to go to lognormal field
        :param M: maximal number of terms in realisation
        :param decay: decay of the field
        :return: void
        #                                           # create a Fenics expression of the affine field
        self.a = Expression('C + S*cos(A*pi*F1*x[0]) * cos(A*pi*F2*x[1])', A=2, C=mean, F1=0, F2=0, S=1, degree=5)
        self.mean = mean
        self.expfield = expfield
        self.M = M
        self.scale = scale
        self.decay = decay

    def realisation(self, y, V):
        # type: (List[float], FunctionSpace) -> Function
          evaluate realisation of random field subject to samples y and return interpolation onto FEM space V.
        :param y: list of samples of the RV in [-1, 1]
        :param V: FunctionSpace
        :return: Fenics Function as field realisation

        def indexer(i):
            m1 = np.floor(i/2)
            m2 = np.ceil(i/2)
            return m1, m2
        assert self.M == len(y)                     # strong assumption but convenient
        a = self.a                                  # store affine field Expression

        a.C, a.S = self.mean, 0                     # get mean function as starting point
        x = interpolate(a, V).vector().array()      # interpolate constant mean on FunctionSpace
        a.C = 0                                     # set mean back to zero. From now on look only at amp_func
        #                                           # add up stochastic modes
        if self.M > 0:
            #                                       # get mean-length ratio
            for m, ym in enumerate(y):              # loop through sample items
                amp = (m+1) ** (-1 * self.decay)    # algebraic decay in front of the sin*sin
                a.F1, a.F2 = indexer(m+2)           # get running values in Expression
                a.S = self.scale                    # multiply a scaling parameter as well
                #                                   # add current Expression value
                x += amp * ym * interpolate(a, V).vector().array()
        f = Function(V)                             # create empty function on FunctionSpace
        #                                           # set function coefficients to realisation coefficients
        f.vector()[:] = x if not self.expfield else np.exp(x)
        return f
In [8]:
def sample_expectation(rv, N):
    mean = 0
    for i in range(N):
        mean += rv.sample()
    mean *= N**(-1)
    return mean

def sample_variance(rv, N):
    var = 0
    mean = sample_expectation(rv, N)
    for i in range(N):
        var += rv.sample()**2
    var *= N**(-1)
    return var - mean**2
In [9]:
mesh = UnitSquareMesh(25, 25)
fs = FunctionSpace(mesh, 'CG', 1)

mean = 1.0
expfield = True
M = 10
scale = 1.0
decay = 2.0
affine_field = CoefField(mean, not expfield, M, scale, decay)
log_field = CoefField(mean, expfield, M, scale, decay)

def sample_field_affine():
    y = np.random.rand(M)*2 - 1
    return affine_field.realisation(y, fs).vector().array()
def sample_field_log():
    y = np.random.randn(M)
    return log_field.realisation(y, fs).vector().array()

affine_field.sample = sample_field_affine
log_field.sample = sample_field_log
In [10]:
def set_fem_fun(vec, fs):
    retval = Function(fs)
    return retval
In [11]:
mean_affine = set_fem_fun(sample_expectation(affine_field, 100), fs)
mean_log = set_fem_fun(sample_expectation(log_field, 100), fs)

mpl.rcParams['figure.figsize'] = [8, 3]
im = plot(mean_affine)
plt.title("Mean of affine field")
im = plot(mean_log)
plt.title("Mean of log normal field")
plt.tight_layout(rect=[0, 0.03, 1, 0.95])
In [12]:
var_affine = set_fem_fun(sample_variance(affine_field, 100), fs)
var_log = set_fem_fun(sample_variance(log_field, 100), fs)

im = plot(var_affine)
plt.title("Variance of affine field")
im = plot(var_log)
plt.title("Variance of log normal field")
plt.tight_layout(rect=[0, 0.03, 1, 0.95])

End of part III. Next: Stochastic Galerkin

© Robert Gruhlke, Manuel Marschall, Phillip Trunschke, 2018-2019