SVM 之 MATLAB 实现代码

MATLAB 中 SVM 实现

直接上代码

• main.m

 
%% Initialize data
clear, clc, close all;
load('data.mat');
 
y(y == 0) = -1;
% X_train = X(1:35, :);
% y_train = y(1:35);
% X_test = X(36:51, :);
% y_test = y(36:51);
 
%% Visualize data
% jhplotdata(X_train, y_train);
 
%% Training a SVM(Support Vector Machine) Classifier
C = 10;
svm = jhsvmtrain(X, y, C, 'Linear');
result = jhsvmtest(svm, X);
fprintf('Accuracy: %f\n', mean(double(result.y_pred == y)));
jhplotdata(X, y);
hold on;
x1_min = min(X(:, 1)) - 1;
x1_max = max(X(:, 1)) + 1;
x2_min = min(X(:, 2)) - 1;
x2_max = max(X(:, 2)) + 1;
 
[XX, YY] = meshgrid(x1_min:0.02:x1_max, x2_min:0.02:x2_max);
Z = jhsvmtest(svm, [XX(:) YY(:)]);
Z = reshape(Z.y_pred, size(XX));
contour(XX, YY, Z);
hold off;

• jhsvmtrain.m

 
function [model] = jhsvmtrain(X, y, C, kernel_type)
%% 函数的核心就是对拉格朗日对偶式的二次规划问题, 通过返回的alpha得到我们需要的支持向量
 
% convert the primal problem to a dual problem, the dual problem is written
% below.
 
% the number of training examples.
m = length(y);
 
% NOTE!! The following two statements, which represent the 
% target function, are fixed cause our target function is fixed.
H = y * y' * jhkernels(X', X', kernel_type);
f = -ones(m, 1);
A = [];
b = [];
Aeq = y';
beq = 0;
lb = zeros(m, 1);
% C is the regularization parameter which means that our model has already
% taken the error into the consideration.
ub = C * ones(m, 1);
alphas0 = zeros(m, 1);
 
epsilon = 1e-8;
options = optimset('LargeScale', 'off', 'Display', 'off');
alphas1 = quadprog(H, f, A, b, Aeq, beq, lb, ub, alphas0, options);
 
logic_vector = abs(alphas1) > epsilon;
model.vec_x = X(logic_vector, :);
model.vec_y = y(logic_vector);
model.alphas = alphas1(logic_vector);
 
end

• jhsvmtest.m

 
function result = jhsvmtest(model, X)
% 在svmTrain中我们主要计算的就是那几个支持向量, 对应的, 核心就是alpha
% 现在有了alpha, 我们通过公式可以轻而易举地计算出w, 我们还不知道b的值, 也即是超平面偏差的值
% 所有先将我们的支持向量代入到公式中, 计算出一个临时的w
% 对于一直的支持向量来说, 我们已经知道了它的label, 所有可以计算出b, 将超平面拉过来, 再将这个b运用到测试集中即可
 
% 带入公式w = \sum_{i=1}^{m}\alpha^{(i)}y^{(i)}x^{(i)}^Tx
% x是输入需要预测的值
tmp = (model.alphas' .* model.vec_y' * jhkernels(model.vec_x', model.vec_x', 'Linear'))';
% 计算出偏差, 也就是超平面的截距
total_bias = model.vec_y - tmp;
bias = mean(total_bias);
 
% 我们已经得到了apha, 因为w是由alpha表示的, 所以通过alpha可以计算出w
% w = sum(alpha .* y_sv)*kernel(x_sv, x_test)
% 其中y_sv是sv的标签, x_sv是sv的样本, x_test是需要预测的数据
w = (model.alphas' .* model.vec_y' * jhkernels(model.vec_x', X', 'Linear'))';
result.w = w;
result.y_pred = sign(w + bias);
result.b = bias;
end

• jhkernel.m

 
function K = jhkernels(X1, X2, kernel_type)
 
switch lower(kernel_type)
    
    case 'linear'
        K = X1' * X2;
    
    case 'rbf'
        K = X1' * X2;
        fprintf("I am sorry about that the rbg kernel is not implemented yet, here we still use the linear kernel to compute\n");    
end
 
end

• jhplotdata.m

 
function jhplotdata(X, y, caption, labelx, labely, color1, color2)
 
if ~exist('caption', 'var') || isempty(caption)
   caption = 'The relationship between X1 and X2';
end
 
if ~exist('labelx', 'var') || isempty(labelx)
   labelx = 'X1';
end
 
if ~exist('labely', 'var') || isempty(labely)
   labely = 'X2';
end
 
if ~exist('color1', 'var') || isempty(color1)
   color1 = 'r'; 
end
 
if ~exist('color2', 'var') || isempty(color2)
   color2 = 'r'; 
end
 
% JHPLOTDATA is going to plot two dimentional data
positive = find(y == 1);
negative = find(y == -1);
 
plot(X(positive, 1), X(positive, 2), 'ro', 'MarkerFace', color1);
hold on;
 
plot(X(negative, 1), X(negative, 2), 'bo', 'MarkerFace', color2);
title(caption);
xlabel(labelx);
ylabel(labely);
legend('Positive Data', 'Negative Data');
 
hold off;
end