盧朓老師《數(shù)值方法：原理、算法及應(yīng)用》課堂 學生大作業(yè)~

供大家參考學習~

相關(guān)代碼和結(jié)果如下：

%The accuracy of the SVM model is 81.11%

修改之后的代碼可以在北太天元上運行, 所有的代碼如下:

load_plugin("optimization");
% 導(dǎo)入數(shù)據(jù)
data = readmatrix('heart disease.csv');
% 檢查并處理缺失值
if any(any(ismissing(data)))
data = rmmissing(data);
end
% 分離特征和標簽
X = data(:,1:end-1);
y = data(:,end);
% 主成分分析降維
[coeff,score,latent] = pca(X);
X_pca = score(:,1:2);
% 劃分訓(xùn)練集和測試集
cv = my_cvpartition(size(X_pca,1),'HoldOut',0.3);
idx = cv.test;
train_idx = setdiff(1: size(X_pca,1), idx);
X_train = X_pca(train_idx,:);
Y_train = y(train_idx,:);
X_test = X_pca(idx,:);
Y_test = y(idx,:);
% 訓(xùn)練SVM模型
C = 1e2; % 設(shè)置較大的C值以確保硬間隔分類（對于線性可分數(shù)據(jù)）
SVMModel = my_fitcecoc(X_train,Y_train,C);
% 預(yù)測測試集
[label,score] = predict(SVMModel,X_test);
% 計算準確率
accuracy = sum(label-1 == Y_test) / length(Y_test);
% 可視化
gscatter(X_pca(:,1),X_pca(:,2),y);
hold on;
% 繪制決策邊界c
d = 0.15;
[x1Grid,x2Grid] = meshgrid(min(X_pca(:,1)):d:max(X_pca(:,1)),...
min(X_pca(:,2)):d:max(X_pca(:,2)));
xGrid = [x1Grid(:),x2Grid(:)];
[~,scores] = predict(SVMModel,xGrid);
contour(x1Grid,x2Grid,reshape(scores(:,2),size(x1Grid)),[0 0],'k');
% 輸出準確率
fprintf('The accuracy of the SVM model is %.2f%%\n', accuracy * 100);
hold off;
function [coeff, score, latent] = pca(X)
% X: 數(shù)據(jù)矩陣，每一列是一個特征，每一行是一個樣本
% coeff: 主成分系數(shù)
% score: 表示主成分得分
% latent: 主成分對應(yīng)的特征值
% 標準化數(shù)據(jù)（均值為0，方差為1）
X = (X - mean(X)) ./ std(X);
% 計算協(xié)方差矩陣
CovMat = cov(X);
% 對協(xié)方差矩陣進行特征值分解
[V, D] = eig(CovMat);
% 將特征值按降序排序，并獲取對應(yīng)的特征向量
[latent, order] = sort(diag(D), 'descend');
coeff = V(:, order);
% 計算主成分得分
score = X * coeff;
end
function cv = my_cvpartition(n, method, param)
% n: 總樣本數(shù)
% method: 分區(qū)方法，目前僅支持 'HoldOut'
% param: 分區(qū)方法的參數(shù)，對于 'HoldOut'，該參數(shù)為留出的比例
% 初始化輸出結(jié)構(gòu)體
cv = struct('train', [], 'test', []);
if strcmp(method, 'HoldOut')
% 生成一個從 1 到 n 的整數(shù)數(shù)組
idx = 1:n;
% 隨機打亂 idx 數(shù)組
idx = idx(randperm(n));
% 計算留出樣本的數(shù)量
numHoldOut = floor(n * param);
% 分配訓(xùn)練和測試索引
cv.test = idx(1:numHoldOut);
cv.train = idx(numHoldOut+1:end);
else
error('Unsupported partition method.');
end
end
function SVMModel = my_fitcecoc(X_train, Y_train, C)
% X_train: 特征矩陣 (n x d)，其中 n 是樣本數(shù)，d 是特征維度
% Y_train: 標簽向量 (n x 1)
% C: 正則化參數(shù)，控制對錯分樣本的懲罰程度
% 獲取類別數(shù)量
uniqueClasses = unique(Y_train);
numClasses = length(uniqueClasses);
% 初始化 SVM 模型結(jié)構(gòu)體數(shù)組
SVMModel = struct('Classifiers', cell(numClasses*(numClasses-1)/2, 1), ...
'ClassPairs', []);
% 訓(xùn)練一對一 SVM 分類器
classifierIndex = 1;
for i = 1:numClasses
for j = (i+1):numClasses
% 提取當前類別對的訓(xùn)練數(shù)據(jù)
classI = Y_train == uniqueClasses(i);
classJ = Y_train == uniqueClasses(j);
X_train_ij = [X_train(classI, :); X_train(classJ, :)];
Y_train_ij = [ones(sum(classI), 1); -1*ones(sum(classJ), 1)];
% 訓(xùn)練 SVM 分類器
SVMModel.Classifiers{classifierIndex} = my_fitcsvm_soft(X_train_ij, Y_train_ij, C);
% 記錄當前分類器對應(yīng)的類別對
SVMModel.ClassPairs(classifierIndex, :) = [uniqueClasses(i), uniqueClasses(j)];
classifierIndex = classifierIndex + 1;
end
end
end
function [wb] = my_fitcsvm(X, Y)
% X: 特征矩陣 (n x d)，其中 n 是樣本數(shù)，d 是特征維度
% Y: 標簽向量 (n x 1)，取值為 +1 或 -1
load_plugin("optimization");
% 假設(shè)X和Y已經(jīng)定義，如之前的示例
N = size(X, 1); % 數(shù)據(jù)點的數(shù)量
D = size(X,2); % 數(shù)據(jù)的維度
% 將w和b組合成一個向量，以便使用quadprog
% 注意：這里我們將w放在前面，b放在最后
p = rand(D + 1, 1); % 初始猜測解（通常設(shè)為0）
Aeq = []; % 沒有等式約束
beq = [];
% 構(gòu)造二次規(guī)劃的目標函數(shù)和線性不等式約束
H = eye(D+1); % Hessian矩陣（目標函數(shù)的二次項系數(shù)）
f = zeros(D+1,1); % 目標函數(shù)的一次項系數(shù)（對于SVM原始問題，通常為負）
% 線性不等式約束 Ax <= b
% 對于每個數(shù)據(jù)點 (x_i, y_i)，我們有 y_i * (w' * x_i + b) >= 1
A = [-Y.*X, -Y]; % 不等式約束的系數(shù)矩陣
b = -ones(N, 1); % 不等式約束的右側(cè)向量
% 使用quadprog求解二次規(guī)劃問題
options = optimoptions('quadprog','Algorithm','interior-point');
[w_b, fval, exitflag, output] = quadprog(H, f, A, b, Aeq, beq, [], [], p, options);
% 分離出w和b
w = w_b(1:D);
b = w_b(D+1);
wb = struct('w',w,'b',b);
end
function [wb] = my_fitcsvm_soft(X, Y, C)
% X: 特征矩陣 (n x d)，其中 n 是樣本數(shù)，d 是特征維度
% Y: 標簽向量 (n x 1)，取值為 +1 或 -1
% C: 正則化參數(shù)，控制對錯分樣本的懲罰程度
% 假設(shè)X和Y已經(jīng)定義，如之前的示例
N = size(X, 1); % 數(shù)據(jù)點的數(shù)量
D = size(X, 2); % 數(shù)據(jù)的維度
% 初始化松弛變量（slack variables）
xi = zeros(N, 1);
% 將w和b以及松弛變量組合成一個向量，以便使用quadprog
p = zeros(D + 1 + N, 1); % 初始猜測解
% 構(gòu)造二次規(guī)劃的目標函數(shù)和線性不等式約束
H = [diag( [ones(1,D),0] ), zeros(D+1, N); zeros(N, D+1), zeros(N,N)]; % Hessian矩陣
f = [zeros(D+1, 1); C*ones(N, 1)]; % 目標函數(shù)的一次項系數(shù)
% 線性不等式約束 Ax <= b
% 對于每個數(shù)據(jù)點 (x_i, y_i)，我們有 y_i * (w' * x_i + b) >= 1 - xi_i
A = [-Y.*X, -Y, diag( -ones(1,N) )]; % 不等式約束的系數(shù)矩陣
b = -ones(N, 1); % 不等式約束的右側(cè)向量
Aeq = []; % 沒有等式約束
beq = [];
lb = [ -inf*ones(D+1,1); zeros(N,1)];
ub = inf*ones(D+1+N,1);
% 使用quadprog求解二次規(guī)劃問題
options = optimoptions('quadprog','Algorithm','interior-point');
[w_b_xi, fval, exitflag, output] = quadprog(H, f, A, b, Aeq, beq, lb, ub, p, options)
% 分離出w, b和xi
w = w_b_xi(1:D);
b = w_b_xi(D+1);
xi = w_b_xi(D+2:end);
% 構(gòu)造并返回結(jié)構(gòu)體wb，包含w和b
wb = struct('w', w, 'b', b);
end
function [label, score] = predict(SVMModel, X_test)
% SVMModel: 訓(xùn)練好的多類SVM模型，由my_fitcecoc函數(shù)返回
% X_test: 測試集特征數(shù)據(jù)
% label: 預(yù)測的類別標簽
% score: 預(yù)測的得分（可選，這里簡化為投票得分）
% 初始化預(yù)測標簽和得分
[numTest, ~] = size(X_test);
numClasses = length(SVMModel.Classifiers);
label = zeros(numTest, 1);
score = zeros(numTest, numClasses);
% 對每個測試樣本進行預(yù)測
for i = 1:numTest
sample = X_test(i, :);
votes = zeros(1, numClasses+1);
% 使用每個SVM分類器進行預(yù)測，并進行投票
for j = 1:numClasses
svmModel = SVMModel.Classifiers{j};
[~, ~, ~, output] = predictOneVsAll(svmModel, sample);
if output == 1
classIdx1 = SVMModel.ClassPairs(j, 1)+1;
votes(classIdx1) = votes(classIdx1) + 1;
else
classIdx2 = SVMModel.ClassPairs(j, 2)+1;
votes(classIdx2) = votes(classIdx2) + 1;
end
end
% 找到得票最多的類別作為預(yù)測標簽
[~, maxVoteIdx] = max(votes);
label(i) = maxVoteIdx;
score(i, maxVoteIdx) = votes(maxVoteIdx); % 將最高票數(shù)作為得分
end
end
function [predictedLabel, predictedScore, decisionValues, output] = predictOneVsAll(svmModel, sample)
% 使用單個SVM模型進行預(yù)測
% 這里簡化為直接計算決策函數(shù)的值，實際應(yīng)用中可能需要更復(fù)雜的處理
w = svmModel.w;
b = svmModel.b;
decisionValues = dot(w, sample) + b;
if decisionValues >= 0
predictedLabel = 1; % 屬于正類
else
predictedLabel = -1; % 屬于負類
end
predictedScore = abs(decisionValues); % 決策函數(shù)的絕對值作為得分
output = predictedLabel; % 用于投票的輸出
end

---

軟間隔SVM通過引入松弛變量（slack variables）$\xi_i$來允許一些樣本被錯誤分類，從而處理非線性可分的數(shù)據(jù)。這些松弛變量衡量了每個樣本違反約束條件的程度。軟間隔SVM的優(yōu)化問題可以寫成以下形式：

$$\min_{\mathbf{w}, b, \xi} \frac{1}{2} ||\mathbf{w}||^2 + C \sum_{i=1}^{m} \xi_i$$

約束條件為：

$$y_i(\mathbf{w}^T \mathbf{x}_i + b) \geq 1 - \xi_i, \quad i = 1, \ldots, m$$
$$\xi_i \geq 0, \quad i = 1, \ldots, m$$

其中，$\mathbf{w}$是權(quán)重向量，$b$是偏置項，$C$是一個正則化參數(shù)，用于控制對錯分樣本的懲罰程度。$\xi_i$是第$i$個樣本的松弛變量。

這個優(yōu)化問題仍然是一個二次規(guī)劃問題，但目標函數(shù)和約束條件都包含了松弛變量。當數(shù)據(jù)不是完全線性可分時，通過允許一些樣本被錯誤分類（即$\xi_i > 0$），優(yōu)化問題可以找到一個更好的分類超平面。

現(xiàn)在，我們來解釋如何克服可行域為空集的問題。在硬間隔SVM（即不允許任何錯誤分類的SVM）中，如果數(shù)據(jù)不是線性可分的，那么約束條件$y_i(\mathbf{w}^T \mathbf{x}_i + b) \geq 1$對于所有樣本可能都無法同時滿足，導(dǎo)致可行域為空集。

然而，在軟間隔SVM中，由于引入了松弛變量$\xi_i$，約束條件變?yōu)榱?y_i(\mathbf{w}^T \mathbf{x}_i + b) \geq 1 - \xi_i$。這意味著即使某些樣本不能滿足原始的硬間隔約束（即$y_i(\mathbf{w}^T \mathbf{x}_i + b) \geq 1$），但只要它們滿足新的軟間隔約束（即$y_i(\mathbf{w}^T \mathbf{x}_i + b) \geq 1 - \xi_i$），并且松弛變量$\xi_i$足夠小，那么這些樣本就不會對優(yōu)化問題的解產(chǎn)生太大的影響。

因此，通過引入松弛變量和軟間隔約束，軟間隔SVM能夠處理非線性可分的數(shù)據(jù)，并找到一個盡可能將數(shù)據(jù)正確分類的分類超平面。在實際應(yīng)用中，可以使用二次規(guī)劃求解器（如SMO算法）來求解這個優(yōu)化問題。