Continuing the discussion from [Approximate GPs with Spectral Stuff]

(Approximate GPs with Spectral Stuff)

Checking the code of this wonderful post: Its correctly implemented. I’m struggling

with the **inverting** the approximation matrix produced by `approx_L.R`

.

If I calculated the eigenvalues with `eigen(L%*%t(L))`

. They are off

the ones in the paper. The approximation however is *accurate*.

I need to use the inverse to calculate the `Analytical Form of Joint Predictive Inference`

.

The `Predictive Inference with a Gaussian Process`

doesn’t lead to good results

in my case.

Where is my mistake?

The paper is referenced here:

Here’s my verification so far.

```
library(mpoly)
cov_exp_quad = function(x, sigma, l) {
outer(x, x, function(xi, xj) { sigma^2 * exp(-(xi - xj)^2 / (2 * l^2)) })
}
approx_L.R = function(M, scale, x, sigma, l) {
epsilon = sqrt(1 / (2 * l^2))
alpha = 1 / scale
beta = (1 + (2 * epsilon / alpha)^2)^0.25
delta = sqrt(alpha^2 * (beta^2 - 1) / 2)
N = length(x)
Ht = matrix(0, nrow = N, ncol = M)
xp = alpha * beta * x
f = sqrt(epsilon^2 / (alpha^2 + delta^2 + epsilon^2))
Ht[, 1] = sqrt(sqrt(alpha^2 / (alpha^2 + delta^2 + epsilon^2))) * sqrt(beta) * exp(-delta^2 * x * x)
Ht[, 2] = f * sqrt(2) * xp * Ht[, 1]
for(n in 3:M) {
Ht[, n] = f * sqrt(2 / (n - 1)) * xp * Ht[, n - 1] - f^2 * sqrt((n - 2) / (n - 1)) * Ht[, n - 2]
}
sigma * Ht
}
l <- 1.0 /sqrt(2)
scale <- 1.0
epsilon = 1 / (sqrt(2) * l);
alpha = 1 / scale;
beta = (1.0 + (2.0 * epsilon / alpha)^2)^0.25;
beta <- (1.0 + (2.0*epsilon/alpha)^2)^0.25
delta = alpha* sqrt((beta^2 - 1.0) / 2.0);
M <- 10
la<-sqrt(alpha^2/(alpha^2+delta^2+epsilon^2))*(epsilon^2/(alpha^2+delta^2+epsilon^2))^(0:(M-1))
X <- matrix(NA, length(x), M)
for(n in 1:M) {
X[,n] <- sqrt(beta/2^(n-1)/gamma(n))*exp(-delta^2*x^2)*as.function(hermite(n-1,kind="h") )(alpha*beta*x)
}
L2<-X %*% diag(la) %*% t(X)
L2.s<- X %*% diag(1/la) %*% t(X)
L<-approx_L.R(M,scale,x,1.0,l)
cov_exp_quad(x,1.0,l) - L %*% t(L)
L2 - cov_exp_quad(x,1.0,l)
L2.inv <- (cbind(X,matrix(0,100,90))) %*% diag(c(1/la,rep(0,90))) %*% t( cbind(X,matrix(0,100,90)))
L2 %*% L2.inv
```