I am trying to implement the Kernel Ridge Regression algorithm but I am getting some quite strange results. I am afraid that I have made some silly mistakes, so I need your help to find out how to fix them. I quote my code bellow. All suggestions are welcome.
function [MSE_train_mean,MSE_train_std,MSE_test_mean,MSE_test_std,w_star] = perform_krr(X,Y)
% Set number of runs.
runs = 2;
% Set number of folds.
folds = 5;
% Set values and number of powers.
power = 0:9;
powers = size(power,2);
% Set values and number of gammas.
gamma = 2 .^ (-power);
gammas = size(gamma,2);
% Set values and number of sigmas.
sigma = 2 .^ (power);
sigmas = size(sigma,2);
% Initialize auxiliary matrices.
MSE_train = zeros(folds,sigmas,gammas);
MSE_train_optimum = zeros(runs,1);
MSE_test = zeros(folds,sigmas,gammas);
MSE_test_optimum = zeros(runs,1);
w_star_optimum = zeros(runs,round(size(X,1)*4/5));
% Perform runs.
for r=1:runs,
% Split original dataset into training and test set.
split_assignments = cross_val_split(folds,size(X,1));
% Perform folds.
for f=1:folds,
% Assign explanatory variables (x).
x_train = X((split_assignments(:,f)==0),:);
x_test = X((split_assignments(:,f)==1),:);
% Assign respond variable (y).
y_train = Y((split_assignments(:,f)==0),:);
y_test = Y((split_assignments(:,f)==1),:);
% Retrieve size of matrices.
[l_train_n,l_train_m] = size(x_train);
[l_test_n,l_test_m] = size(x_test);
% Perform sigmas.
for s=1:sigmas,
% Construct the Gaussian Kernel for given sigma.
[kernel_train] = calculate_krr_gaussiankernel(x_train,x_train,sigma(s));
% Perform gammas.
for g=1:gammas,
% Compute the Dual Weights for given gamma.
a_star = calculate_krr_dualversion(kernel_train,y_train,gamma(g));
% Apply the learned weights on training set and compute the corresponding
% MSE for the given values of sigma and gamma.
MSE_train(f,s,g) = calculate_krr_dualcost(kernel_train,y_train,a_star);
end
end
end
% Average over folds.
MSE_train_mean = mean(MSE_train,1);
MSE_train_std = std(MSE_train,0,1);
% Reshape both mean and std datasets.
MSE_train_mean = reshape(MSE_train_mean,[sigmas gammas]);
MSE_train_std = reshape(MSE_train_std,[sigmas gammas]);
% Plot mean and std of training set MSEs as a function of sigma and gamma.
figure, mesh(MSE_train_mean), title('Mean of Training Set MSEs vs Gamma and Sigma');
figure, mesh(MSE_train_std), title('Std of Training Set MSEs vs Gamma and Sigma');
% Choose the optimum values for the regularization parameters, given the
% fact that we want a minimal training set error; we use the training set
% error as a guidance to select the regularization parameters.
[V,I] = min(MSE_train_mean(:));
[S,G] = ind2sub(size(MSE_train_mean),I);
% Retrieve the optimal parameters.
optimum_sigma = sigma(S);
optimum_gamma = gamma(G);
% Construct the training and test sets: let the first 4/5 entries be the
% training set and the last 1/5 entries be the test set.
x_train = X(1:l_train_n,1:l_train_m);
x_test = X(l_train_n+1:end,1:l_train_m);
y_train = Y(1:l_train_n,1:1);
y_test = Y(l_train_n+1:end,1:1);
% Construct the optimum Gaussian Kernel of both training and test sets.
[kernel_train_optimum] = calculate_krr_gaussiankernel(x_train,x_train,optimum_sigma);
[kernel_test_optimum] = calculate_krr_gaussiankernel(x_test,x_test,optimum_sigma);
% Compute the optimum Dual Weights of both training and test sets.
a_star_train_optimum = calculate_krr_dualversion(kernel_train_optimum,y_train,optimum_gamma);
a_star_test_optimum = calculate_krr_dualversion(kernel_test_optimum,y_test,optimum_gamma);
w_star_optimum(r,:) = a_star_train_optimum';
% Apply the optimum learned weights on both training and test sets and
% compute the corresponding MSEs.
MSE_train_optimum(r,1) = calculate_krr_dualcost(kernel_train_optimum,y_train,a_star_train_optimum);
MSE_test_optimum(r,1) = calculate_krr_dualcost(kernel_test_optimum,y_test,a_star_test_optimum);
end
% % Calculate the training set's (in-sample) residuals.
% y_hat_train = x_train * a_star_optimum;
% e_train = y_train - y_hat_train;
% % Check that they sum up to zero.
% fprintf('Training set`s (in-sample) residuals sum up to %i.\n',sum(e_train));
% % Finally, plot them.
% figure, scatter(1:1:l_train_n,e_train);
% % Calculate the test set's (out-of-sample) residuals.
% y_hat_test = x_test * a_star_optimum;
% e_test = y_test - y_hat_test;
% % Check that they sum up to zero.
% fprintf('Test set`s (out-of-sample) residuals sum up to %i.\n',sum(e_test));
% % Finally, plot them.
% figure, scatter(1:1:l_test_n,e_test);
% Compute mean and standard deviation of both training and test set MSEs.
MSE_train_mean = mean(MSE_train_optimum,1);
MSE_train_std = std(MSE_train_optimum,0,1);
MSE_test_mean = mean(MSE_test_optimum,1);
MSE_test_std = std(MSE_test_optimum,0,1);
w_star = mean(w_star_optimum);
% Print summary statistics.
fprintf('Over %i runs, %i folds, %i gammas, %i sigmas, and regarding the training set, the mean and standard deviation is %.2f and %.2f respectively.\n',runs,folds,gammas,sigmas,MSE_train_mean,MSE_train_std);
fprintf('Over %i runs, %i folds, %i gammas, %i sigmas, and regarding the test set, the mean and standard deviation is %.2f and %.2f respectively.\n',runs,folds,gammas,sigmas,MSE_test_mean,MSE_test_std);
end
function [K] = calculate_krr_gaussiankernel(Xi,Xj,S)
K = zeros(size(Xi,1),size(Xj,1));
for Ixi = 1:size(Xi,1),
for Ixj = 1:size(Xj,1),
K(Ixi,Ixj) = exp((-norm(Xi(Ixi,:) - Xj(Ixj,:)) .^ 2) ./ (2 * (S .^ 2)));
end
end
end
function [A] = calculate_krr_dualversion(K,Y,G)
A = (K + G * size(Y,1) * eye(size(Y,1))) \ Y;
end
function [C] = calculate_krr_dualcost(K,Y,A)
C = (1 / size(Y,1)) * ((K * A - Y)') * (K * A - Y);
end