function [logmax,dL]=nolin1d_grad_optag(x,M,y)

% NOLIN1D_GRAD_OPTAG            This optimisation script uses only GRADIENT
% information to perform Gaussian Process prior model optimisation for
% data [M y] with one-dimensional explanatory variable, M, 
% and y as the target values. Hessian is approximated 
% by finite-differencing using MATLAB operation.
% The hyperparamters, g represents the length-scale, and V, the noise
% variance.
%
% This optimisation script uses only Gradient information to perform
% optimisation for Gaussian process prior models, by adapting
% hyperparameters to minimise negative the log-likelihood function,
%
%                   logmax = log(det(Q)) + N0*log(y'*inv(Q)*y)
%
% The prior covariance function used is of the form,
%
%                   C(zi,zj) = exp(g*[zi-zj]^2) + v
%
% [LOGMAX,DL] = NOLIN1D_GRAD_OPTAG(x,M,y)
%
% Inputs:   M   : one-dimensional input column vector
%           y   : one-dimensional output column vector
%           x   : hyperparameter for GP, x = [g v] for training
%
% Outputs:  LOGMAX  = log likelihood function value
%           DL      = gradient vector value
%
%
% (C) Gaussian Process Toolbox 1.0
% (C) Copyright 2005-2007, Keith Neo <keith.neo@nuim.ie>
%                          http://www.hamilton.ie

N0=size(M,1);
g=exp(x(1));
V=exp(x(2));
expC=g*(repmat(permute(M,[1,3,2]),[1,N0,1])-repmat(permute(M',[3,2 1]),[N0,1,1])).^2;
expA=exp(-0.5*sum(expC,3));
Q=diag_add(expA,V,N0);
invQ=inv(Q);
invQy=invQ*y;
yQy=y'*invQy;
[cholC,pp]=chol(Q);
if pp==0, logdet=2*sum(log(diag(cholC)));
else logdet=sum(log(svd(Q)));
end
logmax=logdet+N0*log(yQy);
clear cholC;
Q=-0.5*expC.*expA;
pp=trace(invQ);
yQQy=invQy'*invQy;
Ny=N0/yQy;

dL=zeros(2,1);
dL(1,1)=traceAxxB(invQ,Q)-Ny*(invQy'*Q*invQy);
dL(2,1)=V*(pp-Ny*yQQy);

function Q=diag_add(Q,v,N0)
for pp=1:N0
    Q(pp,pp)=Q(pp,pp)+v;
end

function trace_AB=traceAxxB(A,B)
trace_AB=0;
dimAB=size(A);dimAB=dimAB(1,1);
for i=1:dimAB,
    trace_AB=trace_AB+A(i,:)*B(:,i);
end;