rand('state',2);
n=100;
d = 10;
noiselevel = 0.1;
x = rand(n,d);
ftrue = 4 * x(:,1) /3 + pi * sin(pi * x(:,3)) + 8 * x(:,6).^5 + 2/(2.71828 - 1) * exp(x(:,8)) - 5;
y = ftrue + normrnd(0,noiselevel,[n 1]);

K3darray = zeros(n,n,d);
for i = 1:d
K3darray(:,:,i) = newsplinematrix(x(:,i),n);
end

theta = ones(d,1);
Kth = zeros(n,n);
for i = 1:d
Kth = Kth + theta(i) * K3darray(:,:, i);
end

m = 30;
sqrtmse = zeros(1,m);
for j = 1:m
lambda0 = 2^( - j);
bigKth = [Kth + lambda0 * eye(n), ones(n,1); ones(1,n), 0];
cb = bigKth \ [y;0];
fest = [Kth, ones(n,1)] * cb;
sqrtmse(j) = norm(fest - ftrue);
end

[C, bestj] = min(sqrtmse);
lambda0 = 2^( - bestj);


m = 20;
sqrtmse = zeros(1,m);
for iter = 1:m
iter
iter
iter
iter
iter
iter
iter
iter
lambda = 2^(15 - iter);

theta = ones(d,1);
normdiffth = 1;
while normdiffth > 0.01

Kth = zeros(n,n);
for inneriter = 1:d
Kth = Kth + theta(inneriter) * K3darray(:,:, inneriter);
end

bigKth = [Kth + lambda0 * eye(n), ones(n,1); ones(1,n), 0];
cb = bigKth \ [y;0];
c = cb(1:n);
b = cb(n+1);

yb = y - b;

G = zeros(n,d);
for inneriter = 1:d
G(:, inneriter) = K3darray(:,:,inneriter) * c;
end

J = G' * c;
Jlambda = lambda0 * J + lambda * ones(d,1);
H = G' * G;
a = Jlambda / 2 - G' * yb;

oldtheta = theta;
theta = quadprog(H, a, zeros(1,d), Inf, zeros(1,d), 0, zeros(d,1));
normdiffth = norm(theta - oldtheta);
end

fest = [Kth, ones(n,1)] * cb;
sqrtmse(iter) = norm(fest - ftrue);
end

[thesqrtmse, bestiter] = min(sqrtmse);
lambda = 2^(15 - bestiter);
theta = ones(d,1);
normdiffth = 1;
while normdiffth > 0.01

Kth = zeros(n,n);
for inneriter = 1:d
Kth = Kth + theta(inneriter) * K3darray(:,:, inneriter);
end

bigKth = [Kth + lambda0 * eye(n), ones(n,1); ones(1,n), 0];
cb = bigKth \ [y;0];
c = cb(1:n);
b = cb(n+1);

yb = y - b;

G = zeros(n,d);
for inneriter = 1:d
G(:, inneriter) = K3darray(:,:,inneriter) * c;
end

J = G' * c;
Jlambda = lambda0 * J + lambda * ones(d,1);
H = G' * G;
a = Jlambda / 2 - G' * yb;

oldtheta = theta;
theta = quadprog(H, a, zeros(1,d), Inf, zeros(1,d), 0, zeros(d,1));
normdiffth = norm(theta - oldtheta);
end

fest = [Kth, ones(n,1)] * cb;
theta










