# Numerically finding a Mercer expansion for a given covariance kernel

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.

• Not quite sure I follow you, but maybe this will help: en.wikipedia.org/wiki/… – Steve Huntsman Aug 16 '11 at 22:39
• How about discretising the problem and finding the eigenvectors of the Gram matrix of $c$? – Simon Lyons Aug 17 '11 at 9:42
• I am interested to know if Tom has something in mind different than what Simon suggests in his comment -- isn't discretization the standard thing to do in this case? But the upvotes for the question suggest to me that I'm missing some aspect of the problem. – R Hahn Aug 17 '11 at 13:02
• Thanks for the responses, everybody. Let me make an edit explaining why I am not taking this approach. – Tom LaGatta Aug 18 '11 at 1:40
• I guess you're setting $Z_i \sim \mathcal{N}(0,\lambda_i)$, are you? – Simon Lyons Aug 18 '11 at 16:46

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')$.