MATLAB implementaion of kernel ridge regression [closed]

Question

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

Answer 1

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.