AP聚类算法实现三维数据点分类

基于亲和传播(Affinity Propagation, AP)聚类算法的三维数据点分类MATLAB实现。AP聚类是一种无监督学习算法，无需预先指定聚类数量，通过数据点之间的"消息传递"自动确定聚类中心。

function ap_clustering_3d()
    % 生成三维测试数据
    [data, true_labels] = generate_3d_data();
    
    % 计算相似度矩阵
    S = compute_similarity(data);
    
    % 执行AP聚类
    [idx, centers, exemplars] = apcluster(S);
    
    % 可视化聚类结果
    visualize_clusters(data, idx, centers, true_labels);
    
    % 评估聚类性能
    evaluate_clustering(data, idx, true_labels, exemplars);
end

function [data, true_labels] = generate_3d_data()
    % 生成三维测试数据
    rng(42); % 设置随机种子保证可重复性
    
    % 定义4个高斯分布的参数
    num_points = 300; % 总点数
    num_clusters = 4; % 真实聚类数量
    
    % 聚类中心
    centers = [1 1 1; 
               -1 -1 1; 
               1 -1 -1; 
               -1 1 -1];
    
    % 生成数据点
    data = [];
    true_labels = [];
    for i = 1:num_clusters
        cluster_points = centers(i,:) + 0.3 * randn(num_points/num_clusters, 3);
        data = [data; cluster_points];
        true_labels = [true_labels; i * ones(num_points/num_clusters, 1)];
    end
    
    % 添加一些噪声点
    noise_points = 4 * (rand(10, 3) - 2);
    data = [data; noise_points];
    true_labels = [true_labels; zeros(size(noise_points, 1), 1)];
    
    % 绘制原始数据
    figure('Name', '原始三维数据', 'Position', [100, 100, 800, 600]);
    scatter3(data(:,1), data(:,2), data(:,3), 30, true_labels, 'filled');
    title('原始三维数据（带真实标签）');
    xlabel('X'); ylabel('Y'); zlabel('Z');
    colormap(jet(num_clusters+1));
    colorbar('Ticks', 0:num_clusters, 'TickLabels', ['噪声'; cellstr(num2str((1:num_clusters)'))]);
    grid on; axis equal;
    view(30, 30); % 设置视角
end

function S = compute_similarity(data)
    % 计算相似度矩阵（负的欧氏距离平方）
    n = size(data, 1);
    S = zeros(n, n);
    
    for i = 1:n
        for j = 1:n
            if i ~= j
                % 计算欧氏距离平方
                dist_sq = sum((data(i,:) - data(j,:)).^2);
                % 相似度为负的欧氏距离平方
                S(i,j) = -dist_sq;
            else
                % 对角线设置为偏好值（preference）
                % 使用中位数作为偏好值
                S(i,j) = median(S(:));
            end
        end
    end
    
    % 使用中位数作为偏好值
    median_val = median(S(S < 0)); % 只考虑负值部分
    S(1:n+1:end) = median_val;
end

function [idx, centers, exemplars] = apcluster(S)
    % AP聚类算法实现
    % 输入: S - 相似度矩阵
    % 输出: 
    %   idx - 每个点的聚类标签
    %   centers - 聚类中心索引
    %   exemplars - 每个点的代表点索引
    
    % 算法参数
    max_iter = 1000; % 最大迭代次数
    conv_iter = 50;  % 收敛检查迭代次数
    lambda = 0.9;    % 阻尼系数
    
    n = size(S, 1);
    
    % 初始化信息矩阵
    R = zeros(n, n); % 责任矩阵
    A = zeros(n, n); % 可用性矩阵
    
    % 存储收敛状态
    conv_history = zeros(1, max_iter);
    
    % 迭代更新
    for iter = 1:max_iter
        % 保存旧的R和A用于收敛检查
        R_old = R;
        A_old = A;
        
        % 更新责任矩阵 (Responsibility)
        AS = A + S;
        [Y, I] = max(AS, [], 2);
        for i = 1:n
            AS(i, I(i)) = -Inf;
        end
        Y2 = max(AS, [], 2);
        
        R_new = S - repmat(Y, 1, n);
        for i = 1:n
            R_new(i, I(i)) = S(i, I(i)) - Y2(i);
        end
        
        % 应用阻尼
        R = (1 - lambda) * R_new + lambda * R_old;
        
        % 更新可用性矩阵 (Availability)
        Rp = max(R, 0);
        Rp(1:n+1:end) = R(1:n+1:end); % 对角线特殊处理
        
        A_new = repmat(sum(Rp, 1), n, 1) - Rp;
        dA = diag(A_new);
        A_new = min(A_new, 0);
        A_new(1:n+1:end) = dA;
        
        % 应用阻尼
        A = (1 - lambda) * A_new + lambda * A_old;
        
        % 计算收敛指标
        conv_history(iter) = mean(abs(R(:) - R_old(:))) + mean(abs(A(:) - A_old(:)));
        
        % 检查收敛（连续conv_iter次变化小于阈值）
        if iter > conv_iter
            recent_changes = conv_history(iter-conv_iter+1:iter);
            if all(recent_changes < 1e-6)
                fprintf('AP聚类在 %d 次迭代后收敛\n', iter);
                break;
            end
        end
    end
    
    if iter == max_iter
        fprintf('达到最大迭代次数 %d\n', max_iter);
    end
    
    % 确定聚类中心和分配点
    E = A + R; % 判断矩阵
    exemplars = diag(E) > 0; % 代表点（聚类中心）
    centers = find(exemplars); % 聚类中心索引
    
    % 分配每个点到其聚类中心
    [~, idx] = max(S(:, centers), [], 2);
    
    % 处理未分配点（理论上AP会分配所有点）
    unassigned = find(isnan(idx));
    if ~isempty(unassigned)
        fprintf('警告: %d 个点未被分配\n', length(unassigned));
        idx(unassigned) = length(centers) + 1; % 分配到新聚类
    end
    
    fprintf('发现 %d 个聚类中心\n', length(centers));
end

function visualize_clusters(data, idx, centers, true_labels)
    % 可视化聚类结果
    figure('Name', 'AP聚类结果', 'Position', [100, 100, 1200, 500]);
    
    % 子图1：聚类结果
    subplot(1, 2, 1);
    cluster_centers = data(centers, :);
    unique_labels = unique(idx);
    num_clusters = length(unique_labels);
    
    % 为每个聚类分配颜色
    colors = lines(num_clusters);
    
    hold on;
    for i = 1:num_clusters
        cluster_points = data(idx == unique_labels(i), :);
        scatter3(cluster_points(:,1), cluster_points(:,2), cluster_points(:,3), ...
                30, colors(i,:), 'filled');
    end
    
    % 绘制聚类中心
    scatter3(cluster_centers(:,1), cluster_centers(:,2), cluster_centers(:,3), ...
            150, 'k', 'p', 'LineWidth', 2);
    
    title(sprintf('AP聚类结果 (%d个聚类)', num_clusters));
    xlabel('X'); ylabel('Y'); zlabel('Z');
    grid on; axis equal;
    view(30, 30);
    
    % 添加图例
    legend_str = arrayfun(@(x) sprintf('聚类 %d', x), 1:num_clusters, 'UniformOutput', false);
    legend([legend_str, {'聚类中心'}], 'Location', 'best');
    
    % 子图2：与真实标签对比
    subplot(1, 2, 2);
    scatter3(data(:,1), data(:,2), data(:,3), 30, idx, 'filled');
    title('聚类分配');
    xlabel('X'); ylabel('Y'); zlabel('Z');
    colormap(lines(num_clusters));
    colorbar('Ticks', 1:num_clusters, 'TickLabels', arrayfun(@num2str, 1:num_clusters, 'UniformOutput', false));
    grid on; axis equal;
    view(30, 30);
    
    % 创建单独的聚类中心图
    figure('Name', '聚类中心与数据分布', 'Position', [100, 100, 800, 600]);
    hold on;
    
    % 绘制所有数据点（灰色）
    scatter3(data(:,1), data(:,2), data(:,3), 20, [0.7, 0.7, 0.7], 'filled');
    
    % 高亮聚类中心
    scatter3(cluster_centers(:,1), cluster_centers(:,2), cluster_centers(:,3), ...
             200, 'r', 'p', 'LineWidth', 2);
    
    % 连接聚类中心
    for i = 1:size(cluster_centers, 1)
        for j = i+1:size(cluster_centers, 1)
            line([cluster_centers(i,1), cluster_centers(j,1)], ...
                 [cluster_centers(i,2), cluster_centers(j,2)], ...
                 [cluster_centers(i,3), cluster_centers(j,3)], ...
                 'Color', 'k', 'LineStyle', '--', 'LineWidth', 1);
        end
    end
    
    title('聚类中心与数据分布');
    xlabel('X'); ylabel('Y'); zlabel('Z');
    legend('数据点', '聚类中心', 'Location', 'best');
    grid on; axis equal;
    view(30, 30);
end

function evaluate_clustering(data, idx, true_labels, exemplars)
    % 评估聚类性能
    
    % 计算轮廓系数
    s = silhouette(data, idx);
    avg_silhouette = mean(s);
    
    % 计算Calinski-Harabasz指数
    eval = evalclusters(data, idx, 'CalinskiHarabasz');
    
    % 计算戴维斯-布尔丁指数
    db = evalclusters(data, idx, 'DaviesBouldin');
    
    % 显示评估指标
    fprintf('\n聚类性能评估:\n');
    fprintf('平均轮廓系数: %.4f\n', avg_silhouette);
    fprintf('Calinski-Harabasz指数: %.4f\n', eval.CriterionValues);
    fprintf('戴维斯-布尔丁指数: %.4f\n', db.CriterionValues);
    
    % 显示聚类中心信息
    centers = find(exemplars);
    cluster_centers = data(centers, :);
    fprintf('\n聚类中心坐标:\n');
    disp(cluster_centers);
    
    % 绘制轮廓系数图
    figure('Name', '聚类质量评估', 'Position', [100, 100, 800, 400]);
    subplot(1, 2, 1);
    silhouette(data, idx);
    title('轮廓系数');
    
    % 绘制聚类大小分布
    subplot(1, 2, 2);
    cluster_counts = histcounts(idx, 'BinMethod', 'integers');
    bar(cluster_counts);
    title('聚类大小分布');
    xlabel('聚类标签');
    ylabel('点数');
    grid on;
    
    % 计算并显示每个聚类的统计信息
    unique_labels = unique(idx);
    num_clusters = length(unique_labels);
    
    fprintf('\n聚类统计信息:\n');
    fprintf('%-10s %-10s %-10s %-10s\n', '聚类ID', '点数', '平均轮廓', '中心坐标');
    
    for i = 1:num_clusters
        cluster_points = data(idx == unique_labels(i), :);
        cluster_size = size(cluster_points, 1);
        cluster_silhouette = mean(s(idx == unique_labels(i)));
        
        if i <= size(cluster_centers, 1)
            center_coord = sprintf('(%.2f, %.2f, %.2f)', ...
                                  cluster_centers(i,1), ...
                                  cluster_centers(i,2), ...
                                  cluster_centers(i,3));
        else
            center_coord = 'N/A';
        end
        
        fprintf('%-10d %-10d %-10.4f %-10s\n', ...
                unique_labels(i), cluster_size, cluster_silhouette, center_coord);
    end
end

% 运行主函数
ap_clustering_3d();

AP聚类算法详解

1. 算法核心思想

亲和传播(AP)聚类是一种基于"消息传递"的聚类算法，其主要特点包括：

• 无需预先指定聚类数量：算法自动确定最佳聚类数
• 以实际数据点作为聚类中心：聚类中心称为"exemplars"
• 处理非凸形状聚类：适用于各种形状的数据分布
• 对初始点不敏感：结果稳定，不受随机初始化的影响

2. 核心数学原理

AP聚类通过两种消息传递进行迭代：

1. 责任(Responsibility)消息 R(i,k)：

   • 表示点i对点k作为其代表的累积证据
   • 计算公式：\(R(i,k) = S(i,k) - \max_{k' \neq k} \{A(i,k') + S(i,k')\}\)

2. 可用性(Availability)消息 A(i,k)：

   • 表示点k适合作为点i代表的累积证据
   • 计算公式：\(A(i,k) = \min \left(0, R(k,k) + \sum_{i' \notin \{i,k\}} \max(0, R(i',k)) \right)\)

其中S(i,k)是相似度矩阵，通常定义为负的欧氏距离平方：\(S(i,k) = -\|x_i - x_k\|^2\)

参考代码 ap聚类算法实现三维数据点的分类 www.youwenfan.com/contentcng/78315.html

3. 算法实现关键步骤

1. 数据生成：

   • 创建4个三维高斯分布作为基础聚类
   • 添加随机噪声点模拟真实场景

2. 相似度矩阵计算：

   • 使用负的欧氏距离平方作为相似度度量
   • 对角线设置为偏好值(preference)，控制聚类数量

3. 消息传递迭代：

   • 交替更新责任矩阵R和可用性矩阵A
   • 应用阻尼因子(λ=0.9)确保收敛
   • 连续50次迭代变化小于阈值则判定收敛

4. 聚类确定：

   • 结合A和R矩阵确定聚类中心(exemplars)
   • 将每个点分配到相似度最高的聚类中心

5. 结果可视化：

   • 三维散点图展示聚类结果
   • 突出显示聚类中心
   • 对比真实标签与聚类结果

6. 性能评估：

   • 轮廓系数：衡量聚类内聚度和分离度
   • Calinski-Harabasz指数：类间方差与类内方差之比
   • Davies-Bouldin指数：类内距离与类间距离之比

4. 参数调整建议

1. 偏好值(preference)：

   • 控制聚类数量，值越大聚类越多
   • 默认使用相似度矩阵的中位数
   • 可尝试：median(S)、min(S)或max(S)

2. 阻尼因子(lambda)：

   • 范围0.5-1.0，防止振荡
   • 默认0.9，收敛慢时可增大到0.95

3. 收敛条件：

   • conv_iter：连续迭代次数
   • 变化阈值：1e-6（可根据数据规模调整）

5. 实际应用扩展

1. 大规模数据优化：

   • 使用稀疏相似度矩阵
   • 采用近似最近邻(ANN)算法

2. 不同相似度度量：

   • 余弦相似度：S(i,j) = dot(x_i,x_j)/(norm(x_i)*norm(x_j))
   • 曼哈顿距离：S(i,j) = -sum(abs(x_i - x_j))

3. 处理分类变量：

   • 使用Jaccard相似度或汉明距离

4. 动态聚类：

   • 增量式AP聚类处理流数据
   • 结合时间衰减因子处理时序数据

graph TD A[生成三维测试数据] --> B[计算相似度矩阵] B --> C[AP聚类迭代] C --> D[确定聚类中心] D --> E[分配点到聚类] E --> F[可视化结果] F --> G[评估聚类质量]

该实现完整展示了AP聚类在三维数据分类中的应用，适用于各种三维点云数据处理场景，如：

• 三维物体识别与分割
• 地理空间数据分析
• 生物医学图像处理
• 分子结构分析
• 运动轨迹聚类