Машина опорных векторов, использующая трюк ядра с решателем quadprog

Question

Я публикую ниже свою собственную реализацию модели оптимизации SVM с мягким полем для классификации на наборе данных saheart.Для запуска кода необходимы Optimization Toolbox и набор данных saheart (набор данных доступен здесь https://web.stanford.edu/~hastie/ElemStatLearn/datasets/SAheart.data).. В большинстве уроков о kernel trick в качестве примера используются небольшие наборы фиктивных данных с двумя атрибутами. Для меня ониТрудно воспроизвести в реальном наборе данных, понять. Не могли бы вы внести изменения в мой код, чтобы он использовал, например, полиномиальное ядро? Я заметил, что есть, например, fitcsvm и т. д., но я бы не хотел его использовать.

% read data
filename = 'saheart.txt';
delimiterIn = ',';
headerlinesIn = 0;
dataSet = importdata(filename,delimiterIn,headerlinesIn);

% organise labels
Y = dataSet(:,10);
dataSet(:,10) = [];
for i = 1:462
    if Y(i,1) == 0
       Y(i,1) = -1; 
    end
end

% hessian, introduction of decision variables, f
decc = sym ('zeta', [462 1]);
syms omega1 omega2 omega3 omega4 omega5 omega6 omega7 omega8 omega9 b;

% tradeoff
C = 5;
f = 1/462*(sum(decc)) + C*power(omega1, 2) + C*power(omega2, 2) + ... 
    C*power(omega3, 2) + C*power(omega4, 2) + C*power(omega5, 2) + ...
    C*power(omega6, 2) + C*power(omega7, 2) + C*power(omega8, 2) + ...
    C*power(omega9, 2);
x = [omega1, omega2, omega3, omega4, omega5, omega6, omega7, omega8, omega9, b];
slackVariables = decc.';
x = [x slackVariables];
H = hessian(f,x);
H = double(H);
x = x.';
f = zeros(472, 1);
for i = 11:472
    f(i,1) = 1/462;
end

% subjects related to hyperplane parameters
AA = zeros(462, 10); 
for i = 1:462
    AA(i,1) = (-1)*Y(i,1)*dataSet(i,1);
    AA(i,2) = (-1)*Y(i,1)*dataSet(i,2);
    AA(i,3) = (-1)*Y(i,1)*dataSet(i,3);
    AA(i,4) = (-1)*Y(i,1)*dataSet(i,4);
    AA(i,5) = (-1)*Y(i,1)*dataSet(i,5);
    AA(i,6) = (-1)*Y(i,1)*dataSet(i,6);
    AA(i,7) = (-1)*Y(i,1)*dataSet(i,7);
    AA(i,8) = (-1)*Y(i,1)*dataSet(i,8);
    AA(i,9) = (-1)*Y(i,1)*dataSet(i,9);
    AA(i,10) = Y(i,1);
end

% subjects related to slack variables
AB = zeros(462, 462);
for i = 1:462
    AB(i,i) = -1;
end
A = [AA AB];

% contstants vector
b = (-1)*ones(462,1);

% lower bounds
lb = zeros(472,1);
lb(1,1) = -Inf;
lb(2,1) = -Inf;
lb(3,1) = -Inf;
lb(4,1) = -Inf;
lb(5,1) = -Inf;
lb(6,1) = -Inf;
lb(7,1) = -Inf;
lb(8,1) = -Inf;
lb(9,1) = -Inf;
lb(10,1) = -Inf;

% launching the optimization
options = optimoptions('quadprog', 'Algorithm','interior-point-convex','Display','off');
[x,fval,exitflag,output,lambda] = quadprog(H,f,A,b,[],[],lb,[],[],options);

% inserting results to the hyperplane equation
HYPERPLANE_EQUATION = zeros(462,1);
for i = 1:462
    HYPERPLANE_EQUATION(i,1) = x(1,1)*dataSet(i,1) + x(2,1)*dataSet(i,2) + ...    
    x(3,1)*dataSet(i,3) + x(4,1)*dataSet(i,4) + x(5,1)*dataSet(i,5) + ... 
    x(6,1)*dataSet(i,6) + x(7,1)*dataSet(i,7) + x(8,1)*dataSet(i,8) + ...
    x(9,1)*dataSet(i,9) - x(10,1);
end

% computing some kind of accuracy
comparsion = [HYPERPLANE_EQUATION Y];
accuracy = 0;
for i = 1:462
    if sign(comparsion(i,1)) == sign(comparsion(i,2))
        accuracy = accuracy + 1; 
    end
end
accuracy = accuracy / 462