Numerically finding a Mercer expansion for a given covariance kernel

Question

Let $c(r)$ be a nice, continuous function with compact support. For example, $c(r) = \tfrac 1 5 (1-r)^{11} \big( 5 + 55r + 239 r^2 + 429 r^3 \big)$ for $r \in [0,1]$, and $c(r) = 0$ otherwise.

On a rectangular domain $D$ in $\mathbb R^2$, we define $c(x,y) = c(|x-y|)$, which we interpret as the kernel of an integral operator. Consider the integral equation $$\iint_D c(x,y) f(y) \ \mathrm{d}y = \lambda f(x).$$ By Mercer's theorem, we can find a sequence of non-negative eigenvalues $\lambda_k$ and eigenfunctions $f_k(x)$ such that:

• The sequence $f_k$ is an orthonormal basis for $L^2(D)$,

• The eigenfunctions $f_k$ corresponding to non-negative eigenvalues are continuous on $D$, and

• We have the representation $c(x,y) = \sum \lambda_k f_k(x) f_k(y)$.

Question: How can we numerically approximate the first $K$ eigenfunctions $f_1, \cdots, f_K$ and keep the total error within a desired error? (Polynomial or Fourier approximations would suffice)

Edit:A number of commentators have described the natural technique: discretize space, turn the covariance into a matrix, easily solve for the eigenvectors, then reconstruct the functions $f_k(t)$ by interpolation.

Our motivation is to generate a Gaussian random field $\xi$ using the Karhunen-Loève transformation: let $Z_1, \cdots, Z_K$ be independent Gaussian RVs, and define the random field $$\xi(t) = \sum_{k=1}^K Z_k f_k(t).$$ Since the covariance is so smooth, the eigenfunctions will be smooth, and the Gaussian field will be smooth. I'm sure that a competent numerical analyst could use the discretization technique to reconstruct the eigenfunctions up to arbitrary derivatives. Unfortunately, I am an incompetent numerical analyst, and I would surely introduce systematic errors along the way. Yet the problem is so common that I am still hoping somebody could point us to a simple reference which does this well.

To make the question more mathematical precise: the numerical approximations $f_k$ should be within an arbitrarily small $C^3$-error of the actual eigenfunctions.

Answer 1

This is my first answer here and I know it is not very precise. However, I hope my thoughts might help you find an answer to your questions.

As I understand you, your goal is to sample from a Gaussian process with covariance kernel $c(t,t')$. Assuming your process has mean zero, the standard way to sample from $\:\xi$ at base points $t_1,...,t_n$ is to compute the matrix $C_{ij} = c(t_i,t_j)$, compute its Cholesky decomposition and multiply that with samples from i.i.d. Gaussian random variables $Z_1,...,Z_n$. The problem is that once you add a new base point you have to repeat the whole procedure which is computationally costly (I presume this is the motivation for your question).

I don't think that it is possible to numerically compute the eigenfunctions with bounded error. The reason is the following: If you look into the literature about kernel principal component analysis, you find that for a finite discretization the eigenfunctions can be written as a linear combination $f_k(t) = \sum_i \alpha_i c(t_i,t)$. Now, in your example, the kernel is shift-invariant which means that its eigenfunctions are the Fourier basis. I do not see how those would be well approximated by functions of the form $f_k(t) = \sum_i \alpha_i c(t_i,t)$.

However, if you really want to use the kernel you specified above you know that the eigenfunctions are the Fourier basis, so you could just use that. The corresponding eigenvalues can be computed via a Fourier transform on $c(t,t')=c(t-t')$.