MATLAB implementaion of kernel ridge regression [closed]

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 sets (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 sets (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

closed as not a real question by gung♦, whuber♦Sep 29 '12 at 20:45

It's difficult to tell what is being asked here. This question is ambiguous, vague, incomplete, overly broad, or rhetorical and cannot be reasonably answered in its current form. For help clarifying this question so that it can be reopened, visit the help center. If this question can be reworded to fit the rules in the help center, please edit the question.

• I would suggest to provide additional details, especially regarding the "quite strange results" you got. – chl Aug 5 '12 at 20:50
• When I said "quite strange" I meant quite higher than the least squares regression as well as the ridge regression algorithms. But, what do you think about the algorithm itself? – user2295350 Aug 6 '12 at 14:15
• Why do I get a very large table of betas and not one as big as the number of features? – user2295350 Aug 6 '12 at 15:34