Matlab implementation for the following paper
- S. Remes, M. Heinonen, S. Kaski (2017). Non-stationary Spectral Kernels. Accepted for NIPS 2017. Preprint: https://arxiv.org/abs/1705.08736
To use the GSM kernel, we need to first define the kernels for the latent functions: frequencies mu(x), lengthscales ell(x) and mixture weights w(x). This is also provided in a helper function get_hyp_kernel.m.
ell = 1; sigma = 1; omega = 1e-4; % some standard values
hyp_kernel.ell = ell;
hyp_kernel.sigma = sigma;
hyp_kernel.omega = omega;
hyp_kernel.K_w = gausskernel(x, x, ell, sigma, omega);
hyp_kernel.Kw_inv = inv(hyp_kernel.K_w);
hyp_kernel.Lw = chol(hyp_kernel.K_w)';
hyp_kernel.K_mu = gausskernel(x, x, ell, sigma, omega);
hyp_kernel.Kmu_inv = inv(hyp_kernel.K_mu);
hyp_kernel.Lmu = chol(hyp_kernel.K_mu)';
hyp_kernel.K_sigma = gausskernel(x, x, ell, sigma, omega);
hyp_kernel.Ksigma_inv = inv(hyp_kernel.K_sigma);
hyp_kernel.Lsigma = chol(hyp_kernel.K_sigma)';Next we need to define the latent functions for A mixture components, let's use some simple constant functions for this example:
hyp.log_sigma = cell(A,1);
hyp.log_mu = cell(A,1);
hyp.log_noise = log(1e-1);
hyp.log_w = cell(A,1);
freq = [3.5,8];
w = [2 1];
for a = 1:A
hyp.log_mu{a} = logit(freq(a)*ones(size(x)), Fn);
hyp.log_sigma{a} = -1*ones(size(x));
hyp.log_w{a} = log(w(a)*ones(size(x)));
end
hyp.log_noise = log(1);A function init_inputdep.m is also provided for initializing the hyperparameters, which is based on fitting a GMM on the spectrogram of the data.
Now we can use the defined hyperparameters to compute the GSM kernel:
K = inputdep_gibbs(x, x, hyp);
Kxy = inputdep_gibbs(x, y, hyp, hyp_kernel);
% hyp_kernel needed to interpolate latent functions to new inputs "y"To learn the latent functions from data, we use minimize_v2.m. First whiten the latent functions, and unwhiten them after optimization.
hyp = inputdep_whiten_vars(hyp, hyp_kernel);
hyp = minimize_v2(hyp, @nlogp_gibbs, -100, u, x, hyp_kernel);% u is the data with inputs x
hyp = inputdep_unwhiten_vars(hyp, hyp_kernel);With Kronecker inference, we define our inputs x as a cell array:
x = {x1, x2, x3}; % for a 3d input grid
% u contains output values with length |x1|*|x2|*|x3|
hyp = init_kron(x, A, hyp_kernel, noise_var, max_freq); % random initialization with A components
hyp = minimize_v2(hyp, @nlogp_kronecker, -100, u, x, hyp_kernel);