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