clear all, close all
dev = @(x,z) norm(x(:)-z(:)) / max([1,norm(x(:)),norm(z(:))]);
n = 350; x = 2*randn(n,1); sn = 0.02; x = sort(x); 
nu = 800; xu = linspace(min(x)-0.1,max(x)+0.1,nu)'; snu = sn/1e3;
dxu = xu(2)-xu(1);
cov = {@covSEiso}; ell = 0.5; sf = 1.5; hyp = [log(ell);log(sf)];
K = feval(cov{:},hyp,x);
subplot(221), imagesc(K), colorbar
subplot(222), imagesc(K), colorbar
ku = feval(cov{:},hyp,xu,xu(1));                            % Kuu = toeplitz(ku)

d = 3;
D = intermat(xu,x,1);                              % sparse interpolation matrix
Q = D*toeplitz(ku)*D';
g = ku(1)-diag(Q);
Kt = Q + diag(g);
subplot(223), imagesc(K-Kt), colorbar, title('interp')

Ku = feval(cov{:},hyp,xu,x);
ku(1) = ku(1) + snu^2;
Q = Ku'*(toeplitz(ku)\Ku);
g = ku(1)-diag(Q);
Kt = Q + diag(g);
subplot(224), imagesc(K-Kt), colorbar, title('FITC')

% 1d benchmark
ns = 32;   s = linspace(-4,4,ns);
nt = 2000; t = linspace(-3,3,nt);
f = exp(-s.*s/2);
E = intermat(s,t,d);
figure, plot(s,f,'-*'), hold on, plot(t,E*f','r')