Commit a32df447 authored by Emile Contal's avatar Emile Contal

initial commit

parents
This diff is collapsed.
% This file is part of GpOptimization.
%
% GpOptimization is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, version 3 of the License.
%
% GpOptimization is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with GpOptimization. If not, see <http://www.gnu.org/licenses/>.
%
% Copyright (c) by Emile Contal, 2015
function [B] = basis_cst(X)
%%
% Constant mean function
%% Syntax
% B = basis_cst(X)
%% Arguments
% * _X_ matrix _(n, d)_ where _n_ is the number of data points and _d_ is the dimension
%% Outputs
% * _B_ vector _(n, 1)_ of ones
%% See also
% <basis_none.html basis_none>
B = [ones(size(X,1),1)];
end
% This file is part of GpOptimization.
%
% GpOptimization is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, version 3 of the License.
%
% GpOptimization is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with GpOptimization. If not, see <http://www.gnu.org/licenses/>.
%
% Copyright (c) by Emile Contal, 2015
function [B] = basis_none(X)
%%
% Zero mean GP inferance
%% Syntax
% B = basis_none(X)
%% Arguments
% * _X_ matrix _(n, d)_ where _n_ is the number of data points and _d_ is the dimension
%% See also
% <basis_cst.html basis_cst>
B = [];
end
% This file is part of GpOptimization.
%
% GpOptimization is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, version 3 of the License.
%
% GpOptimization is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with GpOptimization. If not, see <http://www.gnu.org/licenses/>.
%
% Copyright (c) by Emile Contal, 2015
function [HP, kernel, noise, fval] = bfgs_search_prior(Xt, Yt, HPini, kfun, basis)
%%
% Optimize prior hyper-parameter with respect to the pseudo-likelihood, using the BFGS algorithm
%% Syntax
% HP = bfgs_search_prior(Xt, Yt, HPini, kfun, basis)
%% Arguments
% * _Xt_ matrix _(n, d)_ where _n_ is the number of data points and _d_ is the dimension
% * _Yt_ vector _(n, 1)_ of noisy observations
% * _HPini_ vector _(h, 1)_ of initial hyper-parameters formatted as : _[log(sf2) log(sn) log(w1) log(w2)]_
% * _kfun_ kernel function such as _<kernel_se.html kernel_se>_
% * _basis_ basis function such as _<basis_none.html basis_none>_
%% Outputs
% * _HP_ vector _(h, 1)_ of locally optimal log hyper-parameters
% * _kernel_ found kernel function
% * _noise_ found noise standard deviation
%% See also
% <gp_loolik.html gp_loolik>
NelderMeadIters = 50;
Ht = basis(Xt);
f = @(hp) nllcost(Xt, Yt, hp(1), hp(2), hp(3:end), kfun, Ht);
% Starts with few Nelder-Mead iterations
[hpopt, fvalNM] = fminsearch(f, HPini, optimset('maxFunEvals',NelderMeadIters,'display','off'));
% BGFS search
[hpopt, fval] = fminunc(f, hpopt,optimset('display','off','LargeScale','off'));
HP = hpopt;
kernel = @(x,y) kfun(x,y, exp(HP(1)), exp(HP(3:end)));
noise = exp(HP(1)+HP(2));
end
function [nll] = nllcost(Xt, Yt, sf2, rsn, W, kfun, Ht)
Ktt = kfun(Xt, Xt, exp(sf2), exp(W));
noise = exp(sf2+rsn);
BayesInv = gp_inf(Ktt, Yt, noise, Ht);
nll = gp_loolik(Ktt, Yt, noise, BayesInv, Ht);
end
\ No newline at end of file
% This file is part of GpOptimization.
%
% GpOptimization is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, version 3 of the License.
%
% GpOptimization is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with GpOptimization. If not, see <http://www.gnu.org/licenses/>.
%
% Copyright (c) by Emile Contal, 2015
function [eNet, PiX, DeltaT, U] = chaining_tree(D2, dmax, step, iMin, iMax, varargin)
%%
% Compute the chaining tree given the canonical distance matrix between point of _X_
%% Syntax
% [eNet, PiX, DeltaT, U] = chaining_tree(D2, dmax, step, iMin, iMax)
% [eNet, PiX, DeltaT, U] = chaining_tree(..., 'Name',Value)
%% Arguments
% * _D2_ matrix _(n, n)_ of canonical squared distance
% * _dmax_ scalar of initial diameter of _X_
% * _step_ scalar > 1 for the geometric decay of $\epsilon_i$
% * _iMin_ integer for the first level to consider
% * _iMax_ integer for the last level to consider
%% Name-Value Pair Arguments
% * _a_ scalar > 1 power to use for the geometric decay in the union bound, e.g. 2
% * _lza_ scalar of logarithm of the Riemann zeta of a, e.g. _log(pi^2/6)_
%% Outputs
% * _eNet_ vector _(1, N)_ of indices of points in the final $\epsilon$-net
% * _PiX_ matrix _(h, n)_ of indices of the closest element of the net for all _h=iMax-iMin_ levels
% * _DeltaT_ matrix _(h, N)_ of diameters of the cells of the net for all _h_ levels
% * _U_ vector _(h, 1)_ of negative log probabilities w.r.t the union bounds for all _h_ levels
%% See also
% <gp_dist.html gp_dist> | <enet_greedy.html enet_greedy>
ip = inputParser;
ip.addOptional('a', 2.2);
ip.addOptional('lza', 0.399141);
ip.parse(varargin{:});
opt = ip.Results;
n = size(D2,1);
Ti = [];
nSteps = iMax-iMin;
PiX = zeros(nSteps,n,'uint32');
DeltaT = inf(nSteps,n);
U = inf(nSteps,1);
for ind=1:nSteps
i = ind+iMin-1;
ei = dmax*step^-i;
[ni,Ti] = enet_greedy(D2,ei,Ti);
U(ind) = log(ni+1) + opt.a*log(i) + opt.lza;
[TiDist2, TiInd] = min(D2(Ti,:), [], 1);
PiX(ind,:) = TiInd;
for j=1:length(Ti)
DeltaTi_j = sqrt(max(TiDist2(TiInd==j)));
if ~isempty(DeltaTi_j)
DeltaT(ind,j) = DeltaTi_j;
end
end
end
eNet = Ti;
end
\ No newline at end of file
% This file is part of GpOptimization.
%
% GpOptimization is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, version 3 of the License.
%
% GpOptimization is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with GpOptimization. If not, see <http://www.gnu.org/licenses/>.
%
% Copyright (c) by Emile Contal, 2015
function [ucb] = chaining_ucb(D2, s2, u, varargin)
%%
% Compute the chaining UCB given the canonical distance matrix between point of _X_
%% Syntax
% ucb = chaining_ucb(D2, s2, u)
% ucb = chaining_ucb(..., 'Name',Value)
%% Arguments
% * _D2_ matrix _(n, n)_ of canonical squared distance
% * _s2_ matrix _(n, 1)_ of posterior variance
% * _u_ scalar for negative log probability
%% Name-Value Pair Arguments
% * _step_ scalar > 1 power to use for the geometric decay of $\epsilon$
% * _cdfinv_ function to upper bound the tail probabilities e.g. _@(u,d) d.*sqrt(2*u)_
%% Ouputs
% * _ucb_ vector _(n, 1)_ such that $P[\sup_{x^\star}f(x^\star)-f(x)-\mu(x^\star)+\mu(x^\star) > ucb(x)+cst] < \exp(-u)$
%% See also
% <gpucb.html gpucb> | <chaining_tree.html chaining_tree>
ip = inputParser;
ip.addOptional('step', 2);
ip.addOptional('cdfinv', @(u,d) sqrt(2)*d.^2*erfcinv(2*exp(-u)));
ip.parse(varargin{:});
opt = ip.Results;
n = size(D2,1);
% chaining tree
dmax = max(sqrt(2*s2));
deltaInd = ceil(-log(sqrt(s2)./dmax)/log(opt.step)); % indices such that an upper bound on Delta(X) is < s(X)
[Tree,PiX,DeltaT,U] = chaining_tree(D2,dmax,opt.step,min(deltaInd),max(deltaInd));
% compute ucb(X)
ucb = zeros(n,1);
L1 = false(n,1); % i > min{ i: Delta_i(X) < s(X) }
DeltaXi1 = inf(n,1);
for i=1:size(PiX,1)
ui = U(i)+u;
DeltaXi = DeltaT(i,PiX(i,:))';
LT = DeltaXi < sqrt(s2);
L0 = LT & (~L1);
ucb(L0) = opt.cdfinv(ui, sqrt(s2(Tree(PiX(i,L0))))); % i = min{ i: Delta_i(X) < s(X) }
ucb(L1) = ucb(L1) + opt.cdfinv(ui, DeltaXi1(L1)); % i > min{ }
DeltaXi1 = DeltaXi;
L1(LT) = true;
end
end
\ No newline at end of file
% This file is part of GpOptimization.
%
% GpOptimization is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, version 3 of the License.
%
% GpOptimization is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with GpOptimization. If not, see <http://www.gnu.org/licenses/>.
%
% Copyright (c) by Emile Contal, 2015
function [R] = cholpsd(X)
%%
% Upper Cholesky decomposition of psd matrix
% Tries to approach numerically $\lim_{\epsilon \to 0} chol(X+\epsilon I)$
%% Syntax
% [R] = cholpsd(X)
%% Arguments
% * _X_ psd matrix _(n, n)_
%% Outputs
% * _R_ upper triangular matrix _(n, n)_ such that _X=R'R_
%% See also
% <solve_chol.html solve_chol>
m = min(min(X));
if m<0; fprintf('Error, X is negative in cholpsd\n'); error(''); end
m = max(eps, m*1e-14);
e = m;
I = eye(size(X));
ok = false;
while ~ok
try
R = chol(X);
ok = true;
catch
% if the Cholesky decomposition failed, try to add a small epsilon on the diagonal
X = X+e*I;
if e > 1e6 * m
fprintf('Warning, adding %f for cholpsd\n', e);
end
e = 10*e;
end
end
end
\ No newline at end of file
% This file is part of GpOptimization.
%
% GpOptimization is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, version 3 of the License.
%
% GpOptimization is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with GpOptimization. If not, see <http://www.gnu.org/licenses/>.
%
% Copyright (c) by Emile Contal, 2015
function M = cummax(R)
%%
% Cumulative maximum as used in simple regrets
%% Syntax
% M = cummax(R)
%% Arguments
% * _R_ matrix _(n,t)_ for _n_ run of length _t_
%% Outputs
% * _M_ matrix _(n,t)_ of maximum along the column _(:,1:i)_ for each row
[d,n] = size(R);
M = zeros(d,n);
m = R(:,1);
M(:,1) = m;
for i=2:n
m = max(m, R(:,i));
M(:,i) = m;
end
end
\ No newline at end of file
% This file is part of GpOptimization.
%
% GpOptimization is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, version 3 of the License.
%
% GpOptimization is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with GpOptimization. If not, see <http://www.gnu.org/licenses/>.
%
% Copyright (c) by Emile Contal, 2015
function [N, eNet] = enet_greedy (D2, e, eNetInit)
%%
% Compute an $\epsilon$-net of _Xt_ given the distance matrix
%% Syntax
% [N, eNet] = enet_greedy(D2, e, eNetInit)
%% Arguments
% * _D2_ matrix _(n, n)_ of canonical squared distance
% * _e_ scalar for $\epsilon$
% * _eNetInit_ vector _(1, N)_ of indices of points in the initial $\epsilon$-net
%% Outputs
% * _N_ scalar for the size of the final $\epsilon$-net
% * _eNet_ vector _(1, N)_ of indices of points in the final $\epsilon$-net
%% See also
% <gp_dist.html gp_dist> | <chaining_tree.html chaining_tree>
e2 = e^2;
[n,~] = size(D2);
L = sparse(D2<=e2);
eNet = eNetInit;
if eNet
nc = ~any(L(:,eNet),2); % not covered
else
nc = true(n,1);
end
while any(nc)
nc_ind = find(nc)';
P = sum(L(nc,nc), 2); % number of points in the ball of center X_i and radius e
[~,maxP] = max(P); % we choose the maximum
isolated = find(P==1)'; % and all the isolated points
eNet = [eNet nc_ind([maxP isolated])];
nc(nc) = ~any(L(nc,eNet),2);
end
N = length(eNet);
end
% This file is part of GpOptimization.
%
% GpOptimization is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, version 3 of the License.
%
% GpOptimization is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with GpOptimization. If not, see <http://www.gnu.org/licenses/>.
%
% Copyright (c) by Emile Contal, 2015
function [D2, Kuv_post] = gp_dist(Kuv, Ktu, Ktv, dKuu, dKvv, BayesInv, Ht, Hu, Hv)
%%
% Canonical GP distance $d^2(u,v) = Var[f(u)-f(v) \mid Xt] = \sigma_t^2(u)+\sigma_t^2(v)-2k(u,v)$
%% Syntax
% D2 = gp_dist(Kuv, Ktu, Ktv, dKuu, dKvv, BayesInv)
% D2 = gp_dist(Kuv, Ktu, Ktv, dKuu, dKvv, BayesInv, Ht, Hu, Hv)
%% Arguments
% * _Kuv_ kernel matrix _(nu, nv)_ between the points of _U_ and _V_
% * _Ktu_ kernel matrix _(nt, nu)_ between the points of _Xt_ and _U_
% * _Ktv_ kernel matrix _(nt, nv)_ between the points of _Xt_ and _V_
% * _dKuu_ vector _(nu, 1)_ of the diagonal kernel between the points of _U_
% * _dKvv_ vector _(nv, 1)_ of the diagonal kernel between the points of _V_
% * _BayesInv_ struct array as returned by _<gp_inf.html gp_inf>(Ht, Ktt, Yt, noise)_
% * _Ht_ matrix _(nt, b)_ basis for the points of _Xt_
% * _Hu_ matrix _(nu, b)_ basis for the points of _U_
% * _Hv_ matrix _(nv, b)_ basis for the points of _V_
%% Outputs
% * _D2_ matrix _(nu, nv)_ of squared distance between _U_ and _V_
%% See also
% <gp_pred.html gp_pred>
if nargin<7; Ht = []; end
if nargin<8; Hu = []; end
if nargin<9; Hv = []; end
[~,s2u] = gp_pred(Ktu, dKuu, BayesInv, Ht, Hu);
[~,s2v] = gp_pred(Ktv, dKvv, BayesInv, Ht, Hv);
% TODO update with basis
if Ht
error('gp_dist is not implemented with basis functions.')
end
Kuv_post = Kuv - Ktu'*solve_chol(BayesInv.RC, Ktv);
D2 = bsxfun(@plus,s2u,bsxfun(@minus,s2v',2*Kuv_post));
D2(D2<eps) = 0;
end
\ No newline at end of file
% This file is part of GpOptimization.
%
% GpOptimization is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, version 3 of the License.
%
% GpOptimization is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with GpOptimization. If not, see <http://www.gnu.org/licenses/>.
%
% Copyright (c) by Emile Contal, 2015
function [mu, sigma2] = gp_downdate(Ktt, Yt, i, BayesInv, Ht)
%%
% Posterior $\mu(x_i)$ and $\sigma^2(x_i)$ given $X_t \setminus \{x_i\}$
%% Syntax
% [mu, sigma2] = gp_downdate(Ktt, Yt, i, BayesInv)
% [mu, sigma2] = gp_downdate(Ktt, Yt, i, BayesInv, Ht)
%% Arguments
% * _Ktt_ kernel matrix _(nt, nt)_ between the points of _Xt_
% * _Yt_ vector _(nt, 1)_ of observations
% * _i_ indice of removed observation
% * _BayesInv_ struct array returned by _<gp_inf.html gp_inf>(Ht, Ktt, Yt, noise)_
% * _Ht_ matrix _(nt, b)_ of basis data as returned by _<basis_cst.html basis_cst>(Xt)_
%% Outputs
% * _mu_ scalar _(1,1)_ posterior mean $E[f(x_i) \mid X_t\setminus x_i, Y_t \setminus y_i]$
% * _sigma2_ scalar _(1,1)_ posterior variance $V[f(x_i) \mid X_t\setminus x_i, Y_t \setminus y_i]$
%% See also
% <gp_pred.html gp_pred> | <gp_loolik.html gp_loolik>
n = size(Ktt,1);
T = [1:i-1 i+1:n];
Yt1 = Yt(T);
% Covariance
Kti = Ktt(T,i);
Kii = Ktt(i,i);
% Cholsky downdates (cf Osborne2010 p216)
RC = BayesInv.RC;
RC11 = RC(1:i-1,1:i-1);
RC13 = RC(1:i-1,i+1:end);
S23 = RC(i,i+1:end);
S33 = RC(i+1:end,i+1:end);
RC33 = cholupdate(S33, S23');
RC1 = [RC11 RC13; zeros(size(RC13')) RC33];
if Ht
Ht1 = Ht(T,:);
Hi = Ht(i,:);
RHCH1 = cholpsd(Ht1'*solve_chol(RC1,Ht1));
% System resolution (cf RasmussenWilliams2006 Ch2 p28 Eq2.42)
Ri = Hi - Ht1'*solve_chol(RC1,Kti);
bet = solve_chol(RHCH1, (Ht1'*solve_chol(RC1, Yt1)));
invCY = solve_chol(RC1, (Yt1 - Ht1*bet));
mu = Hi'*bet + Kti'*invCY;
else
invCY = solve_chol(RC1, Yt1);
mu = Kti'*invCY;
bet = [];
end
% sigma2
if nargout > 1
Vf = RC1'\Kti;
covf = Kii - sum(Vf.*Vf, 1)';
if Ht
Vb = RHCH1'\Ri;
covb = sum(Vb.*Vb, 1)';
sigma2 = covb + covf;
else
sigma2 = covf;
end
end
end
% This file is part of GpOptimization.
%
% GpOptimization is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, version 3 of the License.
%
% GpOptimization is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with GpOptimization. If not, see <http://www.gnu.org/licenses/>.
%
% Copyright (c) by Emile Contal, 2015
function [BayesInv] = gp_inf(Ktt, Yt, noise, Ht)
%%
% Bayesian system resolution for computing posterior of GP given the observations _Yt_ at _Xt_
%% Syntax
% BayesInv = gp_inf(Ktt, Yt, noise)
% BayesInv = gp_inf(Ktt, Yt, noise, Ht)
%% Arguments
% * _Ktt_ matrix _(nt, nt)_ of kernel between the points of _Xt_
% * _Yt_ vector _(nt, 1)_ of observations
% * _noise_ noise standard deviation $\eta^2$
% * _Ht_ matrix _(nt, b)_ of basis data as returned by _<basis_cst.html basis_cst>(Xt)_
%% Outputs
% struct array containing:
%
% * _RC_ upper triangular matrix _(nt,nt)_ of Cholesky decomposition of _Ktt+noise*I_
% * _invCY_ vector _(nt,1)_ solution of $(K + \eta^2 I)^{-1} Y$
% * _beta_ vector _(b,1)_ solution of the basis system
%% See also
% <gp_pred.html gp_pred> | <gp_inf_update.html gp_inf_update>
if nargin<4; Ht = []; end
C = Ktt + noise*eye(size(Ktt)); % (nt x nt)
% Cholesky decomposition
RC = cholpsd(C); % (nt x nt)
if Ht
RHCH = cholpsd(Ht' * (solve_chol(RC, Ht))); % (b x b)
% system resolution
bet = solve_chol(RHCH, (Ht' * solve_chol(RC, Yt))); % (b x 1)
invCY = solve_chol(RC, (Yt - Ht*bet));
else
invCY = solve_chol(RC, Yt);
bet = [];
RHCH = [];
end
BayesInv = struct('RC',RC, 'invCY',invCY, 'bet',bet, 'RHCH',RHCH);
end
% This file is part of GpOptimization.
%
% GpOptimization is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, version 3 of the License.
%
% GpOptimization is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with GpOptimization. If not, see <http://www.gnu.org/licenses/>.
%
% Copyright (c) by Emile Contal, 2015
function BayesInv = gp_inf_update(Ktt12, Ktt22, Yt, noise, BayesInv1, Ht)
%%
% Bayesian system update with new observations at _Xt2_
%% Syntax
% BayesInv = gp_inf_update(Ktt12, Ktt22, Yt, noise, BayInv1)
% BayesInv = gp_inf_update(Ktt12, Ktt22, Yt, noise, BayInv1, Ht)
%% Arguments
% * _Ktt12_ matrix _(nt1, nt2)_ of kernel between the points of _Xt1_ and _Xt2_
% * _Ktt22_ matrix _(nt2, nt2)_ of kernel between the points of _Xt2_
% * _Yt_ vector _(nt1+nt2, 1)_ of observations
% * _noise_ noise standard deviation $\eta^2$
% * _BayInv1_ old system resolution for _Xt1_
% * _Ht_ matrix _(nt1+nt2, b)_ of basis data
%% Outputs
% struct array containing:
%
% * _RC_ upper triangular matrix _(nt1+nt2,nt1+nt2)_ of Cholesky decomposition of _Ktt+noise*I_
% * _invCY_ vector _(nt1+nt2,1)_ solution of $(K + \eta^2 I)^{-1} Y$
% * _beta_ vector _(b,1)_ solution of the basis system
%% See also
% <gp_inf.html gp_inf>
[nt1,nt2] = size(Ktt12);
if size(Ktt22,1)~=nt2 | size(Ktt22,2)~=nt2 | size(Yt,1)~=nt1+nt2
error('Wrong arguments size: Ktt12:(%d,%d), Ktt22:(%d,%d), Yt:(%d,%d)',size(Ktt12),size(Ktt22),size(Yt));
end
if size(BayesInv1.RC,1)~=nt1
error('Arguments size incompatible with previous BayesInv: Ktt12:(%d,%d), RC1:(%d,%d)',size(Ktt12),size(BayesInv1.RC));
end
C12 = Ktt12;
if length(noise)==1
C22 = Ktt22 + noise*eye(size(Ktt22)); % (nt x nt)
else
C22 = Ktt22 + diag(noise); % (nt x nt)
end
% Cholsky updates (cf Osborne2010 p214)
RC1 = BayesInv1.RC;
RC12 = RC1'\C12;
RC22 = chol(C22 - RC12'*RC12);
RC = [RC1 RC12; zeros(size(RC12')) RC22]; % (nt x nt)
if Ht
RHCH = cholpsd(Ht'*(solve_chol(RC, Ht))); % (b x b)
% system resolution (cf RasmussenWilliams2006 Ch2 p28 Eq2.42)
bet = solve_chol(RHCH, (Ht' * solve_chol(RC, Yt))); % (b x 1)
invCY = solve_chol(RC, (Yt-Ht*bet));
else
invCY = solve_chol(RC, Yt);
bet = [];
RHCH = [];
end
BayesInv = struct('RC',RC, 'invCY',invCY, 'bet',bet, 'RHCH',RHCH);
end
% This file is part of GpOptimization.
%
% GpOptimization is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, version 3 of the License.
%
% GpOptimization is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with GpOptimization. If not, see <http://www.gnu.org/licenses/>.
%
% Copyright (c) by Emile Contal, 2015
function [nll] = gp_lik(Ktt, Yt, BayesInv, Ht)
%%
% Negative log likelihood given oberservations _Yt_ at _Xt_
%% Syntax
% nll = gp_lik(Ktt, Yt, BayesInv)
% nll = gp_lik(Ktt, Yt, BayesInv, Ht)
%% Arguments
% * _Ktt_ matrix _(nt, nt)_ of kernel between the points of _Xt_
% * _Yt_ vector _(nt, 1)_ of observations
% * _BayesInv_ structure array returned by _<gp_inf.html gp_inf>(Ht, Ktt, Yt, noise)_
% * _Ht_ matrix _(nt, b)_ of basis data as returned by _<basis_cst.html basis_cst>(Xt)_
%% Outputs
% * _nll_ float, negative of the logarithm of the likelihood
%% See also
% <gp_loolik.html gp_loolik> | <bfgs_search_prior.html bfgs_search_prior>
if nargin<4; Ht = []; end
n = size(Yt,1);
m = rank(Ht);
if m
RK = cholpsd(Ktt);
A = Ht'*(solve_chol(RK, Ht));
RA = cholpsd(A);
KHAHK = solve_chol(RK, Ht*solve_chol(RA, Ht'*solve_chol(RK, eye(n))));
nll = .5*Yt'*BayesInv.invCY - .5*Yt'*KHAHK*Yt ...
+ sum(log(diag(RK))) + sum(log(diag(A))) + .5*(n-m)*log(2*pi);
else
nll = .5*Yt'*(BayesInv.invCY) + .5*sum(log(diag(BayesInv.RC))) + .5*n*log(2*pi);
end
end
% This file is part of GpOptimization.
%
% GpOptimization is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, version 3 of the License.
%
% GpOptimization is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with GpOptimization. If not, see <http://www.gnu.org/licenses/>.
%
% Copyright (c) by Emile Contal, 2015
function [nll] = gp_loolik(Ktt, Yt, noise, BayesInv, Ht)
%%
% Negative log leave-one-out likelihood also called Pseudo-likelihood given oberservations _Yt_ at _Xt_
%% Syntax
% nll = gp_loolik(Ktt, Yt, noise, BayesInv)
% nll = gp_loolik(Ktt, Yt, noise, BayesInv, Ht)
%% Arguments