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);


% 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);


% % 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);


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)));

function [A] = calculate_krr_dualversion(K,Y,G)
    A = (K + G * size(Y,1) * eye(size(Y,1))) \ Y;

function [C] = calculate_krr_dualcost(K,Y,A)
    C = (1 / size(Y,1)) * ((K * A - Y)') * (K * A - Y);

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.

  • 2
    $\begingroup$ I would suggest to provide additional details, especially regarding the "quite strange results" you got. $\endgroup$ – chl Aug 5 '12 at 20:50
  • $\begingroup$ 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? $\endgroup$ – user2295350 Aug 6 '12 at 14:15
  • $\begingroup$ Why do I get a very large table of betas and not one as big as the number of features? $\endgroup$ – user2295350 Aug 6 '12 at 15:34

the number of "features" is the dimensionality in the kernel-space ("feature-space"), not the number of input variables.

You should get as many betas as you have features, which is likely far more than you have variables.

  • 1
    $\begingroup$ Couldn't you also technically have far less features, if you had a highly co-linear dataset? More features would surely run the risk of over fitting the dataset, although the original question may relate to an under fitting problem. $\endgroup$ – analystic Sep 29 '12 at 7:09

Not the answer you're looking for? Browse other questions tagged or ask your own question.