MATLAB实现的SVDD算法

MATLAB实现的支持向量数据描述(Support Vector Data Description, SVDD)算法。SVDD是一种单分类方法，用于异常检测和新颖性检测。

代码

classdef SVDD
    %SVDD 支持向量数据描述算法实现
    %   用于单分类和异常检测
    
    properties
        % 模型参数
        C % 正则化参数
        kernelType % 核函数类型：'linear', 'gaussian', 'polynomial'
        kernelParam % 核函数参数
        alpha % 拉格朗日乘子
        SV % 支持向量
        SV_alpha % 支持向量对应的alpha
        radius % 球体半径
        center % 球体中心（特征空间中）
        bias % 决策函数偏置项
        K % 核矩阵
    end
    
    methods
        function obj = SVDD(C, kernelType, kernelParam)
            %SVDD 构造函数
            %   C: 正则化参数
            %   kernelType: 核函数类型
            %   kernelParam: 核函数参数
            
            if nargin < 1
                C = 1;
            end
            if nargin < 2
                kernelType = 'gaussian';
            end
            if nargin < 3
                if strcmp(kernelType, 'gaussian')
                    kernelParam = 1.0;
                else
                    kernelParam = 2;
                end
            end
            
            obj.C = C;
            obj.kernelType = kernelType;
            obj.kernelParam = kernelParam;
        end
        
        function obj = fit(obj, X)
            %FIT 训练SVDD模型
            %   X: 训练数据 (n×d矩阵，n个样本，d个特征)
            
            n = size(X, 1);
            
            % 计算核矩阵
            obj.K = obj.computeKernelMatrix(X, X);
            
            % 设置二次规划问题
            H = obj.K;
            f = -diag(obj.K);
            A = [];
            b = [];
            Aeq = ones(1, n);
            beq = 1;
            lb = zeros(n, 1);
            ub = obj.C * ones(n, 1);
            
            % 使用quadprog求解二次规划
            options = optimset('Display', 'off', 'LargeScale', 'off');
            alpha = quadprog(H, f, A, b, Aeq, beq, lb, ub, [], options);
            
            % 保存结果
            obj.alpha = alpha;
            
            % 找出支持向量
            sv_idx = find(alpha > 1e-6 & alpha < obj.C - 1e-6);
            obj.SV = X(sv_idx, :);
            obj.SV_alpha = alpha(sv_idx);
            
            % 计算半径
            if ~isempty(sv_idx)
                % 使用任意支持向量计算半径
                dist_sv = obj.computeDistanceToCenter(X(sv_idx(1), :), X, alpha);
                obj.radius = sqrt(dist_sv);
            else
                % 如果没有支持向量，使用最大alpha的样本
                [~, max_idx] = max(alpha);
                dist_max = obj.computeDistanceToCenter(X(max_idx, :), X, alpha);
                obj.radius = sqrt(dist_max);
            end
            
            % 计算中心（特征空间中）
            obj.center = (alpha' * X)';
        end
        
        function distances = computeDistanceToCenter(obj, x, X, alpha)
            %COMPUTEDISTANCETOCENTER 计算样本到中心的距离
            %   x: 待计算样本
            %   X: 训练数据
            %   alpha: 拉格朗日乘子
            
            k_x = obj.computeKernelVector(x, X);
            k_xx = obj.computeKernel(x, x);
            
            % 计算距离平方: K(x,x) - 2 * sum(alpha_i * K(x_i, x)) + sum(alpha_i alpha_j K(x_i, x_j))
            distances = k_xx - 2 * (alpha' * k_x') + alpha' * obj.K * alpha;
        end
        
        function [scores, labels] = predict(obj, X_test, threshold)
            %PREDICT 预测新样本
            %   X_test: 测试数据
            %   threshold: 决策阈值（可选）
            %   返回:
            %     scores: 异常得分（距离中心的距离）
            %     labels: 预测标签（1:正常, -1:异常）
            
            n_test = size(X_test, 1);
            scores = zeros(n_test, 1);
            
            for i = 1:n_test
                scores(i) = sqrt(obj.computeDistanceToCenter(X_test(i, :), obj.SV, obj.SV_alpha));
            end
            
            % 如果没有提供阈值，使用训练得到的半径
            if nargin < 3
                threshold = obj.radius;
            end
            
            % 生成标签
            labels = ones(n_test, 1);
            labels(scores > threshold) = -1;
        end
        
        function k = computeKernel(obj, x1, x2)
            %COMPUTEKERNEL 计算两个样本的核函数值
            
            switch obj.kernelType
                case 'linear'
                    k = x1 * x2';
                case 'gaussian'
                    sigma = obj.kernelParam;
                    k = exp(-norm(x1 - x2)^2 / (2 * sigma^2));
                case 'polynomial'
                    d = obj.kernelParam;
                    k = (x1 * x2' + 1)^d;
                otherwise
                    error('不支持的核函数类型');
            end
        end
        
        function K_vec = computeKernelVector(obj, x, X)
            %COMPUTEKERNELVECTOR 计算一个样本与一组样本的核函数值
            
            n = size(X, 1);
            K_vec = zeros(n, 1);
            
            for i = 1:n
                K_vec(i) = obj.computeKernel(x, X(i, :));
            end
        end
        
        function K_mat = computeKernelMatrix(obj, X1, X2)
            %COMPUTEKERNELMATRIX 计算两个样本集之间的核矩阵
            
            n1 = size(X1, 1);
            n2 = size(X2, 1);
            K_mat = zeros(n1, n2);
            
            for i = 1:n1
                for j = 1:n2
                    K_mat(i, j) = obj.computeKernel(X1(i, :), X2(j, :));
                end
            end
        end
        
        function plotDecisionBoundary(obj, X, X_test)
            %PLOTDECISIONBOUNDARY 可视化决策边界（仅适用于2维数据）
            
            if size(X, 2) ~= 2
                error('只能可视化2维数据');
            end
            
            % 创建网格
            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(:)];
            
            % 预测网格点的得分
            [scores, ~] = obj.predict(grid_points);
            
            % 绘制决策边界
            figure;
            contourf(xx, yy, reshape(scores, size(xx)), 20, 'LineStyle', 'none');
            hold on;
            colorbar;
            
            % 绘制训练数据
            scatter(X(:, 1), X(:, 2), 50, 'ko', 'filled');
            
            % 绘制支持向量
            if ~isempty(obj.SV)
                scatter(obj.SV(:, 1), obj.SV(:, 2), 100, 'ro', 'LineWidth', 2);
            end
            
            % 绘制测试数据（如果提供）
            if nargin > 2 && ~isempty(X_test)
                [~, labels] = obj.predict(X_test);
                normal_idx = labels == 1;
                anomaly_idx = labels == -1;
                
                if any(normal_idx)
                    scatter(X_test(normal_idx, 1), X_test(normal_idx, 2), 70, 'go', 'filled');
                end
                
                if any(anomaly_idx)
                    scatter(X_test(anomaly_idx, 1), X_test(anomaly_idx, 2), 70, 'mx', 'LineWidth', 2);
                end
            end
            
            % 绘制球体中心（在原始空间中）
            scatter(obj.center(1), obj.center(2), 150, 'bp', 'filled', 'LineWidth', 2);
            
            % 绘制球体边界（在原始空间中近似）
            theta = 0:0.01:2*pi;
            x_circle = obj.radius * cos(theta) + obj.center(1);
            y_circle = obj.radius * sin(theta) + obj.center(2);
            plot(x_circle, y_circle, 'b-', 'LineWidth', 2);
            
            title('SVDD决策边界');
            xlabel('特征1');
            ylabel('特征2');
            legend('异常得分', '训练数据', '支持向量', '正常测试点', '异常测试点', '中心', '决策边界');
            hold off;
        end
    end
end

使用

% 示例1: 简单二维数据
% 生成训练数据（正态分布）
rng(42); % 设置随机种子以确保可重复性
X_train = [mvnrnd([0, 0], eye(2), 100); 
           mvnrnd([3, 3], 0.5*eye(2), 50)];

% 生成测试数据（包含正常和异常样本）
X_test_normal = mvnrnd([0, 0], eye(2), 20);
X_test_anomaly = [mvnrnd([6, 6], eye(2), 10);
                  mvnrnd([-4, -4], eye(2), 10)];
X_test = [X_test_normal; X_test_anomaly];

% 创建并训练SVDD模型
svdd = SVDD(0.1, 'gaussian', 1.0);
svdd = svdd.fit(X_train);

% 预测
[scores, labels] = svdd.predict(X_test);

% 可视化
svdd.plotDecisionBoundary(X_train, X_test);

% 显示结果
fprintf('检测到的异常比例: %.2f%%\n', sum(labels == -1)/length(labels)*100);

% 示例2: 使用真实数据集（需要Statistics and Machine Learning Toolbox）
% 加载Fisher的鸢尾花数据集
load fisheriris;
X = meas(:, 1:2); % 只使用前两个特征以便可视化
y = strcmp(species, 'setosa'); % 将setosa视为正常类别

% 分割数据
rng(42);
idx = randperm(size(X, 1));
train_ratio = 0.7;
train_size = round(train_ratio * size(X, 1));

X_train = X(idx(1:train_size), :);
y_train = y(idx(1:train_size));

X_test = X(idx(train_size+1:end), :);
y_test = y(idx(train_size+1:end));

% 训练SVDD模型（只使用正常样本）
X_normal = X_train(y_train == 1, :);
svdd_iris = SVDD(0.1, 'gaussian', 0.5);
svdd_iris = svdd_iris.fit(X_normal);

% 预测
[scores_iris, labels_iris] = svdd_iris.predict(X_test);

% 评估性能
accuracy = sum(labels_iris == (y_test*2-1)) / length(y_test);
fprintf('准确率: %.2f%%\n', accuracy*100);

% 可视化
svdd_iris.plotDecisionBoundary(X_normal, X_test);

算法说明

SVDD基本原理

支持向量数据描述(SVDD)是一种单分类算法，其基本思想是通过在特征空间中找到一个最小超球体，使得该球体包含尽可能多的目标数据点。算法的核心优化问题可以表示为：

min R² + C∑ξ_i
s.t. ||φ(x_i) - a||² ≤ R² + ξ_i, ξ_i ≥ 0

其中：

• R是球体半径
• a是球心
• ξ_i是松弛变量，允许一些样本位于球体外
• C是正则化参数，控制松弛变量的惩罚程度

核技巧

通过使用核函数，SVDD可以处理非线性可分的数据。代码中实现了三种常见核函数：

1. 线性核：\(K(x, y) = x·y\)
2. 高斯核(RBF)：\(K(x, y) = exp(-||x-y||²/(2σ²))\)
3. 多项式核：\(K(x, y) = (x·y + 1)^d\)

参考代码 matlab实现的svdd算法 www.3dddown.com/cna/55138.html

关键步骤

1. 计算核矩阵：预先计算所有样本对之间的核函数值
2. 求解二次规划问题：使用MATLAB的quadprog函数求解对偶问题
3. 确定支持向量：选择拉格朗日乘子大于0的样本作为支持向量
4. 计算半径：使用任意支持向量计算球体半径
5. 预测：计算新样本到球心的距离，与半径比较判断是否为异常

参数选择

• C值：控制模型复杂度和泛化能力之间的平衡。较小的C值允许更多样本被分类为异常，较大的C值强制模型包含更多样本。
• 核参数：如高斯核的σ，影响模型的复杂度和灵活性。