Skip to content

Commit 7210ea5

Browse files
amarquandamarquand
committed
updates to BLR
- Add non-ARD case (single lengthscale) - Fixed cosmetic bug with the interpretation of the input parameters (precision matrix was referred to as a covariance matrix)
1 parent ec3dff6 commit 7210ea5

15 files changed

+974
-12
lines changed

‎blr.m

Lines changed: 97 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,97 @@
1+
function [varargout] = blr(hyp, X, t, xs)
2+
3+
% Bayesian linear regression
4+
%
5+
% Fits a bayesian linear regression model, where the inputs are:
6+
% hyp : vector of hyperparmaters. hyp = [log(beta); log(alpha)]
7+
% X : N x D data matrix
8+
% t : N x 1 vector of targets
9+
% xs : Nte x D matrix of test cases
10+
%
11+
% The hyperparameter beta is the noise precision and alpha is the precision
12+
% over lengthscale parameters. This can be either a scalar variable (a
13+
% common lengthscale for all input variables), or a vector of length D (a
14+
% different lengthscale for each input variable, derived using an automatic
15+
% relevance determination formulation).
16+
%
17+
% Two modes are supported:
18+
% [nlZ, dnlZ, post] = blr(hyp, x, y); % report evidence and derivatives
19+
% [mu, s2, post] = blr(hyp, x, y, xs); % predictive mean and variance
20+
%
21+
% Written by A. Marquand
22+
23+
if nargin<3 || nargin>4
24+
disp('Usage: [nlZ dnlZ] = blr(hyp, x, y);')
25+
disp(' or: [mu s2 ] = blr(hyp, x, y, xs);')
26+
return
27+
end
28+
29+
[N,D] = size(X);
30+
beta = exp(hyp(1)); % noise precision
31+
alpha = exp(hyp(2:end)); % weight precisions
32+
Nalpha = length(alpha);
33+
if Nalpha ~= 1 && Nalpha ~= D
34+
error('hyperparameter vector has invalid length');
35+
end
36+
37+
if Nalpha == D
38+
%Sigma = diag(alpha); % weight prior precision
39+
%iSigma = diag(1./alpha); % weight prior covariance
40+
Sigma = diag(1./alpha); % weight prior covariance
41+
iSigma = diag(alpha); % weight prior precision
42+
else
43+
Sigma = 1./alpha*eye(D); % weight prior covariance
44+
iSigma = alpha*eye(D); % weight prior precision
45+
end
46+
47+
XX = X'*X;
48+
A = beta*XX + iSigma; % posterior precision
49+
Q = A\X';
50+
m = beta*Q*t; % posterior mean
51+
52+
if nargin == 3
53+
nlZ = -0.5*( N*log(beta) - N*log(2*pi) - log(det(Sigma)) ...
54+
- beta*(t-X*m)'*(t-X*m) - m'*iSigma*m - log(det(A)) );
55+
56+
if nargout > 1 % derivatives?
57+
dnlZ = zeros(size(hyp));
58+
b = (eye(D) - beta*Q*X)*Q*t;
59+
60+
% noise precision
61+
dnlZ(1) = -( N/(2*beta) - 0.5*(t'*t) + t'*X*m + beta*t'*X*b - 0.5*m'*XX*m ...
62+
- beta*b'*XX*m - b'*iSigma*m -0.5*trace(Q*X) )*beta;
63+
64+
% variance parameters
65+
for i = 1:Nalpha
66+
if Nalpha == D % use ARD?
67+
dSigma = zeros(D);
68+
%dSigma(i,i) = 1 % if alpha is the variance
69+
dSigma(i,i) = -alpha(i)^-2; % if alpha is the precision
70+
else
71+
dSigma = -alpha(i)^-2*eye(D);
72+
end
73+
74+
F = -iSigma*dSigma*iSigma;
75+
c = -beta*F*X'*t;
76+
77+
dnlZ(i+1) = -( -0.5*trace(iSigma*dSigma) + beta*t'*X*c - beta*c'*XX*m ...
78+
- c'*iSigma*m - 0.5*m'*F*m - 0.5*trace(A\F) )*alpha(i);
79+
end
80+
post.m = m;
81+
post.A = A;
82+
end
83+
if nargout > 1
84+
varargout = {nlZ, dnlZ, post};
85+
else
86+
varargout = {nlZ};
87+
end
88+
89+
else % prediction mode
90+
ys = xs*m;
91+
s2 = 1/beta + diag(xs*(A\xs'));
92+
post.m = m;
93+
post.A = A;
94+
varargout = {ys, s2, post};
95+
end
96+
97+
end

‎blr.m~

Lines changed: 93 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,93 @@
1+
function [varargout] = blr(hyp, X, t, xs)
2+
3+
% Bayesian linear regression
4+
%
5+
% Fits a bayesian linear regression model, where the inputs are:
6+
% hyp : vector of hyperparmaters. hyp = [log(sn2); log(alpha)]
7+
% X : N x D data matrix
8+
% t : N x 1 vector of targets
9+
% xs : Nte x D matrix of test cases
10+
%
11+
% The post
12+
%
13+
% Two modes are supported:
14+
% [nlZ, dnlZ, post] = blr(hyp, x, y); % report evidence and derivatives
15+
% [mu, s2, post] = blr(hyp, x, y, xs); % predictive mean and variance
16+
%
17+
% Written by A. Marquand
18+
19+
if nargin<3 || nargin>4
20+
disp('Usage: [nlZ dnlZ] = blr(hyp, x, y);')
21+
disp(' or: [mu s2 ] = blr(hyp, x, y, xs);')
22+
return
23+
end
24+
25+
[N,D] = size(X);
26+
beta = exp(hyp(1)); % noise precision
27+
alpha = exp(hyp(2:end)); % weight precisions
28+
Nalpha = length(alpha);
29+
if Nalpha ~= 1 && Nalpha ~= D
30+
error('hyperparameter vector has invalid length');
31+
end
32+
33+
if Nalpha == D
34+
%Sigma = diag(alpha); % weight prior precision
35+
%iSigma = diag(1./alpha); % weight prior covariance
36+
Sigma = diag(1./alpha); % weight prior covariance
37+
iSigma = diag(alpha); % weight prior precision
38+
else
39+
Sigma = 1./alpha*eye(D); % weight prior covariance
40+
iSigma = alpha*eye(D); % weight prior precision
41+
end
42+
43+
XX = X'*X;
44+
A = beta*XX + iSigma; % posterior precision
45+
Q = A\X';
46+
m = beta*Q*t; % posterior mean
47+
48+
if nargin == 3
49+
nlZ = -0.5*( N*log(beta) - N*log(2*pi) - log(det(Sigma)) ...
50+
- beta*(t-X*m)'*(t-X*m) - m'*iSigma*m - log(det(A)) );
51+
52+
if nargout > 1 % derivatives?
53+
dnlZ = zeros(size(hyp));
54+
b = (eye(D) - beta*Q*X)*Q*t;
55+
56+
% noise precision
57+
dnlZ(1) = -( N/(2*beta) - 0.5*(t'*t) + t'*X*m + beta*t'*X*b - 0.5*m'*XX*m ...
58+
- beta*b'*XX*m - b'*iSigma*m -0.5*trace(Q*X) )*beta;
59+
60+
% variance parameters
61+
for i = 1:Nalpha
62+
if Nalpha == D % use ARD?
63+
dSigma = zeros(D);
64+
%dSigma(i,i) = 1 % if alpha is the variance
65+
dSigma(i,i) = -alpha(i)^-2; % if alpha is the precision
66+
else
67+
dSigma = -alpha(i)^-2*eye(D);
68+
end
69+
70+
F = -iSigma*dSigma*iSigma;
71+
c = -beta*F*X'*t;
72+
73+
dnlZ(i+1) = -( -0.5*trace(iSigma*dSigma) + beta*t'*X*c - beta*c'*XX*m ...
74+
- c'*iSigma*m - 0.5*m'*F*m - 0.5*trace(A\F) )*alpha(i);
75+
end
76+
post.m = m;
77+
post.A = A;
78+
end
79+
if nargout > 1
80+
varargout = {nlZ, dnlZ, post};
81+
else
82+
varargout = {nlZ};
83+
end
84+
85+
else % prediction mode
86+
ys = xs*m;
87+
s2 = 1/beta + diag(xs*(A\xs'));
88+
post.m = m;
89+
post.A = A;
90+
varargout = {ys, s2, post};
91+
end
92+
93+
end

‎blr.old.m

Lines changed: 66 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,66 @@
1+
function [varargout] = blr(hyp, X, t, xs)
2+
3+
% Bayesian linear regression
4+
%
5+
% Fits a bayesian linear regression model, where the inputs are:
6+
% X is an N x D data matrix
7+
% t is an N x 1 vector of targets
8+
% xs is an Nte x D matrix of test cases
9+
10+
if nargin<3 || nargin>4
11+
disp('Usage: [nlZ dnlZ] = blr(hyp, x, y);')
12+
disp(' or: [mu s2 ] = blr(hyp, x, y, xs);')
13+
return
14+
end
15+
16+
[N,D] = size(X);
17+
beta = exp(hyp(1)); % noise precision
18+
alpha = exp(hyp(2:end)); % weight precisions
19+
Sigma = diag(alpha); % weight prior covariance
20+
iSigma = diag(1./alpha); % weight prior precision
21+
22+
if size(X,2) ~= D
23+
24+
end
25+
26+
XX = X'*X;
27+
A = beta*XX + iSigma; % posterior precision
28+
Q = A\X';
29+
m = beta*Q*t; % posterior mean
30+
31+
if nargin == 3
32+
nlZ = -0.5*( N*log(beta) - N*log(2*pi) - log(det(Sigma)) ...
33+
- beta*(t-X*m)'*(t-X*m) - m'*iSigma*m - log(det(A)) );
34+
35+
if nargout > 1 % derivatives?
36+
dnlZ = zeros(size(hyp));
37+
b = (eye(D) - beta*Q*X)*Q*t;
38+
39+
% noise precision
40+
dnlZ(1) = -( N/(2*beta) - 0.5*(t'*t) + t'*X*m + beta*t'*X*b - 0.5*m'*XX*m ...
41+
- beta*b'*XX*m - b'*iSigma*m -0.5*trace(Q*X) )*beta;
42+
43+
% variance parameters
44+
for i = 1:D
45+
dSigma = zeros(D); dSigma(i,i) = 1;
46+
47+
F = -iSigma*dSigma*iSigma;
48+
c = -beta*F*X'*t;
49+
50+
dnlZ(i+1) = -( -0.5*trace(iSigma*dSigma) + beta*t'*X*c - beta*c'*XX*m ...
51+
- c'*iSigma*m - 0.5*m'*F*m - 0.5*trace(A\F) )*alpha(i);
52+
end
53+
post.m = m;
54+
post.A = A;
55+
end
56+
varargout = {nlZ, dnlZ, post};
57+
58+
else % prediction mode
59+
ys = xs*m;
60+
s2 = 1/beta + diag(xs*(A\xs'));
61+
post.m = m;
62+
post.A = A;
63+
varargout = {ys, s2, post};
64+
end
65+
66+
end

‎sp_blr_cluster_job.m

Lines changed: 67 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,67 @@
1+
function [NLML, DNLML, Hyp, Yhat, S2, Yhattr, S2tr] = sp_blr_cluster_job(hyp0,X,Y,opt,Xs)
2+
3+
ones(10)*ones(10); % stupid hack to get matlab to work properly
4+
5+
T = size(Y,2); % number of tasks
6+
7+
% -----------------------------
8+
% defaults
9+
% -----------------------------
10+
try opt.type2ml; catch, opt.type2ml = true; end
11+
try opt.maxEval; catch, opt.maxEval = 100; end
12+
try opt.debug; catch, opt.debug = false;end
13+
14+
D = size(X,2);
15+
16+
Hyp = zeros(T,length(unwrap(hyp0)));
17+
NLML = zeros(T,1);
18+
DNLML = zeros(length(unwrap(hyp0)),T);
19+
20+
if nargin > 4 && nargout > 2
21+
N = size(X,1);
22+
Ns = size(Xs,1);
23+
Yhat = zeros(Ns,T);
24+
S2 = zeros(Ns,T);
25+
Yhattr = zeros(N,T);
26+
S2tr = zeros(N,T);
27+
end
28+
for t = 1:T
29+
if opt.debug; fprintf('processing case %d of %d ...\n',t,T); end
30+
y = Y(:,t);
31+
hyp = zeros(D+1,1);
32+
nlml = NaN;
33+
34+
if opt.type2ml
35+
try
36+
%[hyp,nlml] = minimize(hyp, @gp, opt.maxEval, opt.inf, opt.mean, opt.cov, opt.lik, X, y);
37+
[hyp,nlml] = minimize(zeros(D+1,1), @blr, opt.maxEval, X, y);
38+
39+
% check gradients
40+
fun = @(lh)blr(lh,X,y);
41+
[~,g] = blr(zeros(D+1,1),X,y);
42+
gnum = computeNumericalGradient(fun,zeros(D+1,1));
43+
catch
44+
warning('Optimisation failed. Using default values');
45+
end
46+
end
47+
if nargin > 4
48+
%[yhat, s2] = gp(hyp,opt.inf,opt.mean,opt.cov,opt.lik, X, y, Xs, zeros(Ns,1));
49+
[yhat, s2] = blr(hyp, X, y, Xs);
50+
51+
Yhat(:,t) = yhat;
52+
S2(:,t) = s2;
53+
if nargout > 5
54+
%[yhattr, s2tr] = gp(hyp,opt.inf,opt.mean,cov,opt.lik, X, y, X, zeros(N,1));
55+
[yhattr, s2tr] = blr(hyp, X, y, X);
56+
Yhattr(:,t) = yhattr;
57+
S2tr(:,t) = s2tr;
58+
end
59+
else % just report marginal likelihood and derivatives
60+
%[nlml, dnlml] = gp(hyp,opt.inf,opt.mean,opt.cov,opt.lik, X, y);
61+
%DNLML(:,t) = unwrap(dnlml);
62+
[nlml,DNLML(:,t)] = blr(hyp, X, y);
63+
end
64+
65+
NLML(t) = min(nlml);
66+
Hyp(t,:) = hyp';%unwrap(hyp)';
67+
end

‎sp_blr_cluster_job.sh

Lines changed: 32 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,32 @@
1+
#!/bin/sh
2+
# script for execution of deployed applications
3+
#
4+
# Sets up the MCR environment for the current $ARCH and executes
5+
# the specified command.
6+
#
7+
exe_name=$0
8+
exe_dir=`dirname "$0"`
9+
echo "------------------------------------------"
10+
if [ "x$1" = "x" ]; then
11+
echo Usage:
12+
echo $0 \<deployedMCRroot\> args
13+
else
14+
echo Setting up environment variables
15+
MCRROOT="$1"
16+
echo ---
17+
LD_LIBRARY_PATH=.:${MCRROOT}/runtime/glnxa64 ;
18+
LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:${MCRROOT}/bin/glnxa64 ;
19+
LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:${MCRROOT}/sys/os/glnxa64;
20+
export LD_LIBRARY_PATH;
21+
echo LD_LIBRARY_PATH is ${LD_LIBRARY_PATH};
22+
shift 1
23+
args=
24+
while [ $# -gt 0 ]; do
25+
token=$1
26+
args="${args} \"${token}\""
27+
shift
28+
done
29+
eval "\"${exe_dir}/sp_blr_cluster_job_comp\"" $args
30+
fi
31+
exit
32+

‎sp_blr_cluster_job_cfun.mat

353 Bytes
Binary file not shown.

‎sp_blr_cluster_job_comp

6.28 MB
Binary file not shown.

0 commit comments

Comments
 (0)