Noiseless Hilbert Space Gaussian Process?

Modeling
1

So the vanilla GP:

data{
	int n_x ;
	array[n_x] real x ;
	vector[n_x] y ;
}
parameters{
	real<lower=0> noise ;
	real<lower=0> amplitude ;
	real<lower=0> lengthscale ;
	real intercept ;
	vector[n_x] f ;
}
model{
	... // priors for noise, amplitude, lengthscale, & intercept here
	f ~ multi_normal( // n.b. "centered" parameterization, but often non-centered samples better
		rep_vector(intercept,n_x)
		, add_diag( gp_exp_quad_cov( x , amplitude , lengthscale ) , 1e-5 )
	) ;
	y ~ normal( f , noise ) ;
}

can be altered to remove measurement noise (something usually not a good idea, but useful in niche circumstances) to yield:

... // same data
parameters{
	// no `noise` parameter
	real<lower=0> amplitude ;
	real<lower=0> lengthscale ;
	real intercept ;
	// no `f` parameter
}
model{
	... // priors for amplitude, lengthscale, & intercept here
	y ~ multi_normal( 
		rep_vector(intercept,n_x)
		, add_diag( gp_exp_quad_cov( x , amplitude , lengthscale ) , 1e-5 )
	) ;
}

Now if one opts for the cool new Hilbert space GPs, the noise-included model would be:

functions{
	matrix PHI(int N, int n_basis_fns, real L, vector x) {
		return sin(diag_post_multiply(rep_matrix(pi()/(2*L) * (x+L), n_basis_fns), linspaced_vector(n_basis_fns, 1, n_basis_fns)))/sqrt(L);
	}
	vector diagSPD_EQ(real amplitude, real lengthscale, real L, int n_basis_fns) {
		return amplitude * sqrt(sqrt(2*pi()) * lengthscale) * exp(-0.25*(lengthscale*pi()/2/L)^2 * linspaced_vector(n_basis_fns, 1, n_basis_fns)^2);
	}
}

data{
	int n_x ;
	vector[n_x] x ;
	vector[n_x] y ;
	real<lower=0> boundary_factor ;
	int n_basis_fns ;
}

transformed data{
	vector[n_x] x_std = (x-mean(x))/sd(x) ;
	real L_f = boundary_factor*max(abs(x_std));
	matrix[n_x,n_basis_fns] PHI_f = PHI(n_x, n_basis_fns, L_f, x_std);
}

parameters{
	real<lower=0> noise ;
	real<lower=0> lengthscale ;
	real<lower=0> amplitude ;
	real intercept ;
	vector[n_basis_fns] beta ;
}

model{
	... // priors for noise, amplitude, lengthscale, & intercept here
	beta ~ normal( // n.b. "centered" parameterization
		rep_vector( 0 , n_basis_fns )
		, diagSPD_EQ( amplitude , lengthscale , L_f , n_basis_fns ) 
	) ;
	vector[n_x] f = intercept + PHI_f * beta ;
	y ~ normal( f , noise ) ;
}

However, I’m unclear what the “noiseless” version of this HSGP would be.

Any ideas?

2 Likes
2

You still have the noise 1e-5 there, so not yet completely removed.

In theory, just set noise to be sqrt(1e-5) (assuming that is close enough for noiseless as in the covariance matrix based code), but with explicit basis function coefficients beta that will produce a very challenging posterior geometry.

In case of normal observation model, you could integrate beta out analytically and present the linear model as GP with dot product covariance matrix. While the original gp_exp_quad_cov covariance matrix approach has a computation cost of O(n_x^3), the basis function dot product covariance matrix approach has a computation cost of O(n_basis_fns^3).