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

% NOLIN1D_HESS_OPTAG_SCHUR_R    This optimisation script uses only GRADIENT
% information to perform Gaussian Process prior model optimisation for
% time-series data, i.e. one-dimensional data, [M,y] with M as the input
% and y as the output. Hessian information is supplied by user.
% The hyperparamters, g represents the length-scale, and V, the noise
% variance.
%
% This is a fast algorithm using Generalised Schur algorithm to perform
% computation task, impossible by standard MATLAB command. The Schur
% algorithm includes Hyperbolic rotation.
%
% [logmax,dL,he] = NOLIN1D_HESS_OPTAG_SCHUR_R(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
%           he      = hessian matrix
%
%
% (C) Gaussian Process Schur 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));
hyp=[g V];
G=get_eigen_block_trace_a(M,hyp);
[trdQ,bdQb,dQb,invb,logdet,trQ,yQy,yQQy]=schur_invbtrace3_r(G,y);
logmax=logdet+N0*log(yQy);
dL=zeros(2,1);
tmp1=bdQb/yQy;
tmp2=yQQy/yQy;
dL(1,1)=trdQ-N0*tmp1;
dL(2,1)=V*(trQ-N0*tmp2);
Ny=N0/yQy;

coef=0.01;
G=get_eigen_block_traceqq_a1(M,hyp,coef);
[trq2qqxx,dQQbxx]=schur_trace3qq_r(G,y);
trq2qq=(trq2qqxx-trQ)/coef;
dQQb=(dQQbxx-invb)/coef;
he=zeros(2,2);
he(2,2)=dL(2,1)-V*V*(trq2qq-2*Ny*(invb'*dQQb)+N0*tmp2*tmp2);
coef=0.1;
GBB=get_eigen_block_tracedqdq_a3(M,hyp,coef);   % coeff.norm=0.5
[trdqdq,dqdqb,trq3dqdq]=schur_trace3dqdq_ar(GBB,y,coef);
clear GBB;
trqdq=(trQ-trq3dqdq)/coef;
he(2,1)=-V*(trqdq+N0*tmp1*tmp2-2*Ny*(dQQb'*dQb));
he(1,2)=he(2,1);
G=get_eigen_block_trace_ahess(M,hyp);
[trdQ2,dQb2]=schur_invbtrq3_hessian_r(G,y);
he(1,1)=trdQ2-trdqdq+Ny*(2*(invb'*dqdqb)-(invb'*dQb2))-N0*tmp1*tmp1;