test.m

%% 主程序 - SVM实验
clear all; close all; clc;

% 设置随机数种子确保可重复性
rng(42);

% 生成数据
[X_sep, y_sep] = generate_data(true);  % 线性可分数据
[X_nonsep, y_nonsep] = generate_data(false); % 线性不可分数据

% 实验1: 线性可分SVM
fprintf('%s\n', repmat('=', 1, 50));
fprintf('线性可分SVM实验\n');
fprintf('%s\n', repmat('=', 1, 50));
svm_sep = train_svm(X_sep, y_sep, 1.0, 1.0, 0.001, 5000);
plot_svm(X_sep, y_sep, svm_sep, '线性可分SVM');

% 实验2: 线性不可分SVM
fprintf('\n%s\n', repmat('=', 1, 50));
fprintf('线性不可分SVM实验\n');
fprintf('%s\n', repmat('=', 1, 50));
svm_nonsep = train_svm(X_nonsep, y_nonsep, 1.0, 1.0, 0.001, 5000);
plot_svm(X_nonsep, y_nonsep, svm_nonsep, '线性不可分SVM');

% 实验3: 参数C的影响分析
fprintf('\n%s\n', repmat('=', 1, 50));
fprintf('参数C对SVM的影响分析\n');
fprintf('%s\n', repmat('=', 1, 50));
parameter_experiment(X_nonsep, y_nonsep);

% 附加实验: 增广拉格朗日法参数分析
fprintf('\n%s\n', repmat('=', 1, 50));
fprintf('增广拉格朗日法参数分析\n');
fprintf('%s\n', repmat('=', 1, 50));
rho_values = [0.01, 0.1, 1, 10];
eta_values = [0.001, 0.01, 0.1];

for rho = rho_values
    for eta = eta_values
        fprintf('\n参数组合: rho=%.3f, eta=%.3f\n', rho, eta);
        try
            svm = train_svm(X_sep, y_sep, 1.0, rho, eta, 1000);
            fprintf('收敛迭代次数: %d\n', svm.iterations);
        catch ME
            fprintf('参数不收敛: %s\n', ME.message);
        end
    end
end

%% 数据生成函数
function [X, y] = generate_data(separable)
    if separable
        centers = [1, 1; 6, 6];  % 线性可分
    else
        centers = [1, 1; 3, 3];  % 线性不可分
    end
    
    % 每类生成100个样本
    data1 = mvnrnd(centers(1,:), [1, 0; 0, 1], 100);
    data2 = mvnrnd(centers(2,:), [1, 0; 0, 1], 100);
    
    % 合并数据并添加标签
    X = [data1; data2];
    y = [ones(100,1)*-1; ones(100,1)]; % 第一类为-1，第二类为1
end

%% SVM训练函数（增广拉格朗日法）
function model = train_svm(X, y, C, rho, eta, max_iter)
    % 正确获取数据维度
    n_samples = size(X, 1);  % 样本数量（行数）
    n_features = size(X, 2); % 特征数量（列数）
    
    % 初始化参数
    alphas = zeros(n_samples, 1);
    lambda_val = 0.0;
    
    % 预计算核矩阵（线性核）
    K = X * X';
    
    % 增广拉格朗日优化
    converged = false;
    for iter = 1:max_iter
        alphas_prev = alphas;
        
        % 计算梯度
        grad = ones(n_samples, 1) - (K * (alphas .* y)) .* y;
        grad = grad + lambda_val * y + rho * sum(alphas .* y) * y;
        
        % 更新alpha
        alphas = alphas - eta * grad;
        
        % 投影到[0, C]区间
        alphas = max(0, min(C, alphas));
        
        % 更新lambda
        lambda_val = lambda_val + rho * sum(alphas .* y);
        
        % 检查收敛
        diff = norm(alphas - alphas_prev);
        if diff < 1e-4
            fprintf('收敛于 %d 次迭代\n', iter);
            model.iterations = iter;
            converged = true;
            break;
        end
    end
    
    if ~converged
        model.iterations = max_iter;
        fprintf('达到最大迭代次数 %d\n', max_iter);
    end
    
    % 计算权重向量
    w = sum((alphas .* y) .* X, 1)';
    
    % 计算偏置项b (使用支持向量)
    sv_indices = alphas > 1e-5;
    sv_X = X(sv_indices, :);
    sv_y = y(sv_indices);
    sv_alphas = alphas(sv_indices);
    
    b_vals = zeros(sum(sv_indices), 1);
    for i = 1:length(sv_y)
        b_vals(i) = sv_y(i) - dot(w, sv_X(i,:));
    end
    b = mean(b_vals);
    
    % 保存模型参数
    model.w = w;
    model.b = b;
    model.alphas = alphas;
    model.C = C;
    model.rho = rho;
    model.eta = eta;
    model.X_train = X;
    model.y_train = y;
end

%% 决策函数
function scores = decision_function(model, X)
    if isfield(model, 'w')
        % 自定义模型
        scores = X * model.w + model.b;
    else
        % MATLAB内置模型
        [~, scores] = predict(model, X);
        scores = scores(:,2); % 获取第二类的分数
    end
end

%% 可视化函数
function plot_svm(X, y, model, title_str)
    figure('Position', [100, 100, 800, 600]);
    hold on;
    
    % 绘制数据点
    scatter(X(y == -1, 1), X(y == -1, 2), 50, 'o', 'MarkerEdgeColor', 'b', 'LineWidth', 1.5);
    scatter(X(y == 1, 1), X(y == 1, 2), 100, '*', 'MarkerEdgeColor', 'r', 'LineWidth', 1.5);
    
    % 获取坐标轴范围
    x_min = min(X(:,1)) - 1; x_max = max(X(:,1)) + 1;
    y_min = min(X(:,2)) - 1; y_max = max(X(:,2)) + 1;
    
    % 创建网格进行预测
    [xx, yy] = meshgrid(linspace(x_min, x_max, 100), linspace(y_min, y_max, 100));
    grid_points = [xx(:), yy(:)];
    
    % 计算决策边界和间隔
    Z = decision_function(model, grid_points);
    Z = reshape(Z, size(xx));
    
    % 绘制决策边界和间隔
    contour(xx, yy, Z, [-1, 0, 1], 'LineWidth', 2);
    
    % 标记支持向量
    if isfield(model, 'w')
        % 自定义模型
        sv_indices = model.alphas > 1e-5;
        sv_X = X(sv_indices, :);
        scatter(sv_X(:,1), sv_X(:,2), 150, 'o', 'MarkerEdgeColor', 'k', 'LineWidth', 1.5, 'MarkerFaceColor', 'none');
        
        fprintf('支持向量数量: %d\n', sum(sv_indices));
        fprintf('偏置项 b: %.4f\n', model.b);
    else
        % MATLAB内置模型
        sv = model.SupportVectors;
        scatter(sv(:,1), sv(:,2), 150, 'o', 'MarkerEdgeColor', 'k', 'LineWidth', 1.5, 'MarkerFaceColor', 'none');
    end
    
    % 设置图例和标题
    legend('Class -1', 'Class 1', 'Support Vectors', 'Location', 'best');
    title(title_str, 'FontSize', 14);
    xlabel('Feature 1', 'FontSize', 12);
    ylabel('Feature 2', 'FontSize', 12);
    grid on;
    axis equal;
    hold off;
end

%% 参数实验分析
function parameter_experiment(X, y)
    C_values = [0.01, 0.1, 1, 10, 100];
    
    figure('Position', [200, 200, 1200, 800]);
    sgtitle('参数 C 对支持向量机和决策边界的影响', 'FontSize', 16);
    
    for i = 1:length(C_values)
        C = C_values(i);
        
        % 使用MATLAB内置SVM训练
        model = fitcsvm(X, y, 'KernelFunction', 'linear', 'BoxConstraint', C);
        
        % 绘制结果
        subplot(2, 3, i);
        hold on;
        
        % 绘制数据点
        scatter(X(y == -1, 1), X(y == -1, 2), 50, 'o', 'MarkerEdgeColor', 'b', 'LineWidth', 1.5);
        scatter(X(y == 1, 1), X(y == 1, 2), 100, '*', 'MarkerEdgeColor', 'r', 'LineWidth', 1.5);
        
        % 获取坐标轴范围
        x_min = min(X(:,1)) - 1; x_max = max(X(:,1)) + 1;
        y_min = min(X(:,2)) - 1; y_max = max(X(:,2)) + 1;
        
        % 创建网格进行预测
        [xx, yy] = meshgrid(linspace(x_min, x_max, 100), linspace(y_min, y_max, 100));
        grid_points = [xx(:), yy(:)];
        
        % 计算决策边界
        Z = decision_function(model, grid_points);
        Z = reshape(Z, size(xx));
        
        % 绘制决策边界
        contour(xx, yy, Z, [-1, 0, 1], 'LineWidth', 2);
        
        % 标记支持向量
        sv = model.SupportVectors;
        scatter(sv(:,1), sv(:,2), 150, 'o', 'MarkerEdgeColor', 'k', 'LineWidth', 1.5, 'MarkerFaceColor', 'none');
        
        title(sprintf('C = %.2f\nSV: %d', C, size(sv,1)), 'FontSize', 12);
        grid on;
        axis equal;
        hold off;
    end
end