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