1. 项目概述：从数据到洞察，层次聚类的C++实践

在数据分析和机器学习的工具箱里，聚类算法扮演着将无序数据点分门别类的角色，而层次聚类（Hierarchical Clustering）因其直观的树状结构（通常称为树状图或谱系图）和无需预先指定聚类数量的特性，一直备受青睐。它不像K-Means那样需要你拍脑袋猜一个K值，而是通过计算数据点之间的“距离”，自底向上（凝聚法）或自顶向下（分裂法）地构建一个层次关系，最终结果是一棵完整的聚类树。你可以像切蛋糕一样，在树的任意高度“切一刀”，就能得到对应粒度的聚类结果，这种灵活性在探索性数据分析中尤其有用。

今天，我们不依赖任何第三方机器学习库，就用纯粹的C++来实现一个完整的凝聚型层次聚类算法。选择C++，一方面是因为它离硬件更近，执行效率高，处理大规模数据时优势明显；另一方面，对于学习算法本质而言，亲手实现一遍远比调用sklearn.cluster.AgglomerativeClustering来得深刻。你会清晰地看到距离矩阵如何计算、最近簇如何寻找与合并、树状结构如何一步步生成。这个过程会涉及指针、动态内存管理、标准模板库（STL）的灵活运用，以及对算法时间、空间复杂度的切身感受。无论你是正在啃《C++ Primer》的学生，还是希望夯实算法基础的开发者，亦或是好奇聚类算法内部机制的爱好者，这篇详解都能带你走通从理论到实现的全过程。

2. 算法核心原理与设计思路拆解

2.1 层次聚类的基本思想与流程

层次聚类的核心思想很简单：开始时，每个数据点自成一簇；然后，寻找距离最近的两个簇，将它们合并成一个新簇；重复这个过程，直到所有数据点合并成一个大簇。这个“合并史”就被记录成了一棵树，也就是树状图。

我们主要实现最常用的凝聚层次聚类（Agglomerative Hierarchical Clustering）。其标准流程可以概括为以下几步：

1. 初始化：将每个样本点视为一个独立的簇。同时，计算所有簇两两之间的距离，形成一个初始的距离矩阵。
2. 循环合并： a. 在当前距离矩阵中，找到距离最小的两个簇（记为簇A和簇B）。 b. 将簇A和簇B合并，形成一个新的簇AB。 c. 更新距离矩阵：删除代表簇A和簇B的行和列，新增一行一列代表新簇AB，并计算新簇AB与所有其他簇之间的距离。
3. 记录与终止：记录本次合并（谁和谁合并，合并时的距离是多少）。重复步骤2，直到所有样本合并为一个簇。

这个流程中，最关键的变量是“距离”的定义。它有两个层面：一是点与点之间的距离（如欧氏距离），二是簇与簇之间的距离。后者定义了不同的聚类“链接准则”，会显著影响最终的聚类形状。

2.2 关键设计决策：链接准则与数据结构

链接准则（Linkage Criterion）决定了如何计算两个簇之间的距离。常见的几种有：

• 单链接（Single Linkage）：两个簇中所有点对之间距离的最小值。容易产生“链式效应”，擅长发现非球形的、拉长的簇，但对噪声点敏感。
• 全链接（Complete Linkage）：两个簇中所有点对之间距离的最大值。倾向于产生紧凑的、大小相近的球状簇，对噪声相对稳健。
• 平均链接（Average Linkage）：两个簇中所有点对之间距离的平均值。是单链接和全链接的折中，也是最常用的一种。
• 沃德法（Ward‘s Method）：合并后导致的簇内方差平方和增量最小的两个簇。倾向于生成大小相近的簇。

提示：对于初次实现，我建议从平均链接开始。它在大多数情况下能取得不错的效果，且计算逻辑清晰，易于实现和调试。

数据结构的设计直接决定了程序的效率和代码的清晰度。我们需要表示：

1. 数据点（Data Point）：一个多维向量。可以用std::vector<double>来表示。
2. 簇（Cluster）：一个簇包含多个数据点的索引。随着合并，簇会动态变化。我们需要一种结构来高效地管理簇的成员，并支持快速的合并操作。
3. 距离矩阵（Distance Matrix）：存储所有簇两两之间的距离。由于簇的数量会从N减少到1，矩阵大小会变化。我们需要一种能高效查找最小值、删除行/列、插入新行/列的数据结构。

一个高效且清晰的策略是：

• 簇的表示：使用std::vector<int>存储该簇包含的所有原始数据点的索引。合并簇时，直接将一个簇的成员列表追加到另一个簇的末尾。同时，我们维护一个cluster_map，将簇的ID映射到其成员列表。
• 距离矩阵的表示：使用二维std::vector<double>。虽然每次合并需要删除两行两列并新增一行一列（O(n²)操作），但对于数据量不是特别大（例如几千个点）的教学和一般应用场景，其实现简单直观的优势更大。我们也可以用一个std::priority_queue（优先队列/最小堆）来辅助快速查找最小距离，但更新操作会变复杂。为了优先保证代码的可读性和教学性，我们首选用二维向量实现完整的矩阵更新。

3. 核心数据结构与工具函数实现

3.1 数据点与簇的定义

首先，我们定义一些类型别名和数据结构，让代码更易读。

#include <vector> #include <cmath> #include <limits> #include <algorithm> #include <iostream> #include <map> // 类型别名，提高代码可读性 using Point = std::vector<double>; // 一个数据点，即一个多维向量 using Dataset = std::vector<Point>; // 整个数据集 using Cluster = std::vector<int>; // 一个簇，存储该簇内数据点的索引 // 合并记录：用于最终生成树状图 struct MergeRecord { int cluster1_id; // 合并时，簇1的ID（在合并时刻的ID） int cluster2_id; // 合并时，簇2的ID double distance; // 合并时的距离 int new_cluster_id; // 合并后新簇的ID（通常可以用cluster1_id或新生成的ID） // 注意：在实际树状图绘制中，还需要记录合并后新簇的大小（包含的样本数） };

3.2 距离计算函数

距离计算是聚类的基础。我们先实现点与点之间的欧氏距离，这是最常用的距离度量。

/** * 计算两个点之间的欧几里得距离（L2范数） * @param p1 点1 * @param p2 点2 * @return 欧氏距离 */ double euclideanDistance(const Point& p1, const Point& p2) { // 确保两个点维度相同 if (p1.size() != p2.size()) { throw std::invalid_argument("Points must have the same dimensionality."); } double sum = 0.0; for (size_t i = 0; i < p1.size(); ++i) { double diff = p1[i] - p2[i]; sum += diff * diff; } return std::sqrt(sum); }

接下来是实现链接准则，即簇间距离的计算。我们以实现平均链接为例：

/** * 计算两个簇之间的平均链接距离 * @param cluster1 簇1，包含数据点索引 * @param cluster2 簇2，包含数据点索引 * @param dataset 原始数据集 * @return 两个簇间所有点对距离的平均值 */ double averageLinkageDistance(const Cluster& cluster1, const Cluster& cluster2, const Dataset& dataset) { if (cluster1.empty() || cluster2.empty()) { return std::numeric_limits<double>::max(); // 或者抛出异常 } double total_distance = 0.0; int pair_count = 0; // 遍历簇1和簇2中的所有点对 for (int idx1 : cluster1) { const Point& p1 = dataset[idx1]; for (int idx2 : cluster2) { const Point& p2 = dataset[idx2]; total_distance += euclideanDistance(p1, p2); pair_count++; } } // 防止除零，虽然理论上不会发生，因为簇非空 return pair_count > 0 ? total_distance / pair_count : std::numeric_limits<double>::max(); }

注意：计算簇间距离是一个双重循环，时间复杂度为O(|C1| * |C2|)。当簇变大时，计算成本会上升。在实际的高性能库中，会使用更高效的方法（如维护距离的递推公式）来更新距离矩阵，而不是每次都重新计算。但为了代码清晰，我们这里采用直接计算的方式。

3.3 距离矩阵的初始化与表示

我们需要一个类或结构来管理距离矩阵以及簇的信息。这里我们设计一个ClusterManager类。

class ClusterManager { private: Dataset dataset_; // 原始数据 std::map<int, Cluster> clusters_; // 当前活跃的簇，键为簇ID，值为簇成员索引 std::vector<std::vector<double>> distance_matrix_; // 当前簇间的距离矩阵 std::vector<int> active_cluster_ids_; // 当前活跃簇的ID列表，与distance_matrix_行列顺序对应 std::vector<MergeRecord> merge_history_; // 合并历史，用于构建树状图 // 内部函数：根据活跃簇ID列表，重新计算或更新距离矩阵 void computeDistanceMatrix() { size_t n = active_cluster_ids_.size(); distance_matrix_.assign(n, std::vector<double>(n, 0.0)); for (size_t i = 0; i < n; ++i) { for (size_t j = i + 1; j < n; ++j) { // 只计算上三角，矩阵是对称的 int id_i = active_cluster_ids_[i]; int id_j = active_cluster_ids_[j]; double dist = averageLinkageDistance(clusters_[id_i], clusters_[id_j], dataset_); distance_matrix_[i][j] = dist; distance_matrix_[j][i] = dist; // 对称赋值 } // 对角线为0，自己到自己的距离为0，初始化时已经是0 } } public: // 构造函数：用数据集初始化，每个点作为一个簇 ClusterManager(const Dataset& dataset) : dataset_(dataset) { int n = dataset.size(); for (int i = 0; i < n; ++i) { clusters_[i] = {i}; // 簇ID为i，初始只包含索引i active_cluster_ids_.push_back(i); } computeDistanceMatrix(); // 初始化距离矩阵 } // ... 其他成员函数将在后续实现 };

这个ClusterManager在构造时完成了算法的第一步：初始化。它创建了N个单点簇，并计算了初始的N×N距离矩阵。

4. 层次聚类核心算法的C++实现

4.1 查找最小距离与簇合并

算法的核心循环是查找当前距离矩阵中的最小距离（对角线除外），然后合并对应的两个簇。我们需要在ClusterManager中添加相应的方法。

首先，实现一个查找最小距离项的函数：

// 在ClusterManager类内部添加 /** * 在当前距离矩阵中查找非零最小距离及其对应的簇索引（在active_cluster_ids_中的位置） * @param min_i [输出] 最小距离对应的第一个簇的矩阵行索引 * @param min_j [输出] 最小距离对应的第二个簇的矩阵列索引 * @param min_dist [输出] 最小距离值 * @return 是否成功找到（当活跃簇少于2个时返回false） */ bool findMinDistance(size_t& min_i, size_t& min_j, double& min_dist) const { min_dist = std::numeric_limits<double>::max(); min_i = 0; min_j = 0; size_t n = active_cluster_ids_.size(); if (n < 2) return false; // 不足两个簇，无法合并 for (size_t i = 0; i < n; ++i) { for (size_t j = i + 1; j < n; ++j) { // 只遍历上三角 if (distance_matrix_[i][j] < min_dist) { min_dist = distance_matrix_[i][j]; min_i = i; min_j = j; } } } return min_dist < std::numeric_limits<double>::max(); }

接下来是实现合并操作。这是算法中最关键也最需要小心的一步，涉及数据结构的更新。

// 在ClusterManager类内部添加 /** * 执行一次合并操作：找到距离最近的两个簇，合并它们，并更新所有相关数据结构。 * @return 如果成功合并返回true，如果只剩一个簇则返回false。 */ bool performMerge() { size_t i, j; double dist; if (!findMinDistance(i, j, dist)) { return false; // 没有可合并的簇了 } // 1. 获取要合并的两个簇的ID int id_i = active_cluster_ids_[i]; int id_j = active_cluster_ids_[j]; // 2. 创建合并记录（重要：在修改数据结构前记录） MergeRecord record; record.cluster1_id = id_i; record.cluster2_id = id_j; record.distance = dist; // 新簇的ID，我们可以约定使用较小的那个ID，或者生成一个新ID。这里使用id_i作为新ID。 record.new_cluster_id = id_i; merge_history_.push_back(record); // 3. 合并簇成员：将簇j的成员全部加入到簇i中 Cluster& cluster_i = clusters_[id_i]; const Cluster& cluster_j = clusters_[id_j]; cluster_i.insert(cluster_i.end(), cluster_j.begin(), cluster_j.end()); // 删除旧的簇j（或者标记为无效，这里直接从map中删除） clusters_.erase(id_j); // 4. 更新活跃簇ID列表：删除j，保留i // 注意：因为我们要从vector中删除元素，且j > i，先删除j再删除i可以防止索引错乱 // 更安全的方法是先保存要删除的索引值，然后按索引从大到小删除 // 这里j > i，所以先删j，再删i active_cluster_ids_.erase(active_cluster_ids_.begin() + j); active_cluster_ids_.erase(active_cluster_ids_.begin() + i); // 将新合并的簇（id_i）的ID加入到列表末尾（或任意位置，但为了简化，放末尾） active_cluster_ids_.push_back(id_i); // 5. 重新计算距离矩阵（这是最耗时的部分，但逻辑清晰） computeDistanceMatrix(); return true; }

4.2 主循环与聚类过程驱动

有了单次合并的函数，主循环就非常简单了：持续合并，直到只剩下一个簇。

// 在ClusterManager类内部添加 /** * 执行完整的层次聚类过程。 */ void fit() { merge_history_.clear(); // 清空历史记录 // 持续合并，直到performMerge返回false（只剩一个活跃簇） while (performMerge()) { // 可以在这里添加一些调试信息，比如打印当前合并步骤 // std::cout << "Merged, remaining clusters: " << active_cluster_ids_.size() << std::endl; } std::cout << "Hierarchical clustering completed. Total merge steps: " << merge_history_.size() << std::endl; }

4.3 获取聚类结果与树状图信息

算法运行完毕后，我们通常需要两种输出：

1. 指定聚类数量K时的扁平化聚类结果：在树状图的某个高度切割，得到K个簇。
2. 完整的合并历史（树状图数据）：用于可视化或进一步分析。

首先，实现根据merge_history_来获取指定簇数量的结果。我们需要从合并历史中反向推导出在形成K个簇时的状态。

// 在ClusterManager类内部添加 /** * 根据合并历史，获取当聚类数量为n_clusters时的聚类标签。 * 标签是一个vector，长度等于样本数，每个元素表示该样本所属的簇ID（从0到n_clusters-1）。 * @param n_clusters 期望的聚类数量 * @return 每个样本的簇标签 */ std::vector<int> getClusterLabels(int n_clusters) const { int n_samples = dataset_.size(); if (n_clusters <= 0 || n_clusters > n_samples) { throw std::invalid_argument("n_clusters must be between 1 and number of samples."); } // 初始化：每个点属于自己的簇，簇ID为其原始索引 std::vector<int> labels(n_samples); for (int i = 0; i < n_samples; ++i) { labels[i] = i; } // 我们需要知道在合并历史的哪一步之后，簇的数量变成了n_clusters。 // 初始有n_samples个簇，每合并一次减少一个簇。 // 因此，当我们执行了 (n_samples - n_clusters) 次合并后，就得到了n_clusters个簇。 int merges_needed = n_samples - n_clusters; // 确保合并历史足够长 merges_needed = std::min(merges_needed, static_cast<int>(merge_history_.size())); // 应用前merges_needed次合并，来更新标签 for (int step = 0; step < merges_needed; ++step) { const MergeRecord& record = merge_history_[step]; int c1 = record.cluster1_id; int c2 = record.cluster2_id; int new_id = record.new_cluster_id; // 这里我们约定新ID为cluster1_id // 将所有属于簇c2的样本的标签，改为new_id（即c1） // 注意：我们需要一个从旧簇ID到样本索引的映射。但我们只有标签数组。 // 更高效的方法是维护一个并查集(Union-Find)结构。这里为了清晰，使用遍历。 // 这是一个O(n_samples * steps)的操作，对于教学和小数据可以接受。 // 在实际应用中，强烈建议使用并查集。 for (int& label : labels) { if (label == c2) { label = new_id; } } // 注意：如果c2在后续合并中又作为new_id出现，我们的标签仍然指向它，这没问题。 // 但为了正确性，我们还需要处理一种情况：后续合并中，c1可能又被合并到其他簇。 // 上面的简单遍历方法在多次合并后可能会出错，因为标签没有“追溯”到最终的根簇。 // 因此，我们需要一个更鲁棒的方法：使用并查集在合并历史中维护簇的归属关系。 } // 此时，labels中的ID是经过merges_needed次合并后的簇ID，但这些ID可能不是连续的0到K-1。 // 我们需要将它们重新映射到连续的整数标签。 std::map<int, int> id_to_label; int new_label = 0; for (int& label : labels) { if (id_to_label.find(label) == id_to_label.end()) { id_to_label[label] = new_label++; } label = id_to_label[label]; } // 理论上，此时new_label应该等于n_clusters return labels; }

实操心得：上面getClusterLabels函数中的标签更新部分是一个简化实现，它在合并步骤较多时可能无法得到正确结果，因为它没有处理簇ID的传递性（即A合并到B，B又合并到C，那么A的标签最终应该是C）。一个工业级的实现必须使用并查集（Union-Find）数据结构。并查集可以在近乎常数时间内完成“查找根节点”和“合并集合”的操作，完美适配层次聚类的合并过程。在fit()过程中，每合并一次，就在并查集中做一次union操作。最后，要得到K个簇的标签，只需对每个样本查找其在并查集中的根，然后将根映射到连续的标签即可。这是层次聚类实现中关于效率与正确性的一个关键点。

其次，获取合并历史以供可视化。

// 在ClusterManager类内部添加 /** * 获取完整的合并历史，用于绘制树状图。 * 通常，树状图绘制库需要：合并顺序、合并时的距离、以及每次合并涉及的簇大小。 * @return 合并历史记录 */ const std::vector<MergeRecord>& getMergeHistory() const { return merge_history_; } // 同时，我们需要在MergeRecord中添加簇大小的信息，这对绘制有意义的树状图很重要。 // 修改MergeRecord结构体： struct MergeRecord { int cluster1_id; int cluster2_id; double distance; int new_cluster_id; int cluster1_size; // 新增：合并时簇1的大小 int cluster2_size; // 新增：合并时簇2的大小 }; // 并在performMerge()中记录这些大小： // record.cluster1_size = clusters_[id_i].size(); // 在合并前记录 // record.cluster2_size = clusters_[id_j].size();

5. 完整代码整合与示例测试

现在，我们将所有部分整合到一个完整的程序中，并用一个简单的二维数据集进行测试。

#include <vector> #include <cmath> #include <limits> #include <algorithm> #include <iostream> #include <map> #include <fstream> #include <sstream> #include <iomanip> // ... 此处插入之前定义的 Point, Dataset, Cluster, MergeRecord, euclideanDistance, averageLinkageDistance ... class ClusterManager { private: Dataset dataset_; std::map<int, Cluster> clusters_; std::vector<std::vector<double>> distance_matrix_; std::vector<int> active_cluster_ids_; std::vector<MergeRecord> merge_history_; void computeDistanceMatrix() { size_t n = active_cluster_ids_.size(); distance_matrix_.assign(n, std::vector<double>(n, 0.0)); for (size_t i = 0; i < n; ++i) { for (size_t j = i + 1; j < n; ++j) { int id_i = active_cluster_ids_[i]; int id_j = active_cluster_ids_[j]; double dist = averageLinkageDistance(clusters_[id_i], clusters_[id_j], dataset_); distance_matrix_[i][j] = dist; distance_matrix_[j][i] = dist; } } } bool findMinDistance(size_t& min_i, size_t& min_j, double& min_dist) const { min_dist = std::numeric_limits<double>::max(); min_i = 0; min_j = 0; size_t n = active_cluster_ids_.size(); if (n < 2) return false; for (size_t i = 0; i < n; ++i) { for (size_t j = i + 1; j < n; ++j) { if (distance_matrix_[i][j] < min_dist) { min_dist = distance_matrix_[i][j]; min_i = i; min_j = j; } } } return true; } public: ClusterManager(const Dataset& dataset) : dataset_(dataset) { int n = dataset.size(); for (int i = 0; i < n; ++i) { clusters_[i] = {i}; active_cluster_ids_.push_back(i); } computeDistanceMatrix(); } bool performMerge() { size_t i, j; double dist; if (!findMinDistance(i, j, dist)) { return false; } int id_i = active_cluster_ids_[i]; int id_j = active_cluster_ids_[j]; MergeRecord record; record.cluster1_id = id_i; record.cluster2_id = id_j; record.distance = dist; record.new_cluster_id = id_i; record.cluster1_size = clusters_[id_i].size(); record.cluster2_size = clusters_[id_j].size(); merge_history_.push_back(record); // 合并簇 Cluster& cluster_i = clusters_[id_i]; const Cluster& cluster_j = clusters_[id_j]; cluster_i.insert(cluster_i.end(), cluster_j.begin(), cluster_j.end()); clusters_.erase(id_j); // 更新活跃ID列表：先删索引大的j，再删i active_cluster_ids_.erase(active_cluster_ids_.begin() + j); active_cluster_ids_.erase(active_cluster_ids_.begin() + i); active_cluster_ids_.push_back(id_i); computeDistanceMatrix(); return true; } void fit() { merge_history_.clear(); int step = 0; while (performMerge()) { step++; if (step % 10 == 0) { // 每10步打印一次进度 std::cout << "Merge step " << step << ", clusters left: " << active_cluster_ids_.size() << std::endl; } } std::cout << "Clustering finished after " << step << " merges." << std::endl; } // 使用并查集改进的getClusterLabels (简化版并查集实现) std::vector<int> getClusterLabels(int n_clusters) const { int n_samples = dataset_.size(); if (n_clusters <= 0 || n_clusters > n_samples) { throw std::invalid_argument("Invalid n_clusters"); } // 初始化并查集，每个样本的父节点是自己 std::vector<int> parent(n_samples); for (int i = 0; i < n_samples; ++i) { parent[i] = i; } // 查找根节点的函数（带路径压缩） std::function<int(int)> find = [&](int x) { if (parent[x] != x) { parent[x] = find(parent[x]); // 路径压缩 } return parent[x]; }; // 应用前 (n_samples - n_clusters) 次合并 int merges_to_apply = n_samples - n_clusters; merges_to_apply = std::min(merges_to_apply, static_cast<int>(merge_history_.size())); for (int step = 0; step < merges_to_apply; ++step) { const MergeRecord& record = merge_history_[step]; int root1 = find(record.cluster1_id); int root2 = find(record.cluster2_id); if (root1 != root2) { // 合并：将root2的根指向root1 (按照记录中的new_cluster_id，这里简化都指向cluster1_id的根) parent[root2] = root1; } } // 现在，每个样本的最终簇标签是其并查集的根 std::vector<int> raw_labels(n_samples); for (int i = 0; i < n_samples; ++i) { raw_labels[i] = find(i); } // 将根ID映射到连续的0到K-1的标签 std::map<int, int> root_to_label; int label_counter = 0; std::vector<int> final_labels(n_samples); for (int i = 0; i < n_samples; ++i) { int root = raw_labels[i]; if (root_to_label.find(root) == root_to_label.end()) { root_to_label[root] = label_counter++; } final_labels[i] = root_to_label[root]; } // 理论上label_counter应等于n_clusters return final_labels; } const std::vector<MergeRecord>& getMergeHistory() const { return merge_history_; } // 一个辅助函数，打印距离矩阵（用于调试） void printDistanceMatrix() const { std::cout << "Current distance matrix:" << std::endl; for (const auto& row : distance_matrix_) { for (double val : row) { std::cout << std::setw(10) << std::setprecision(4) << val << " "; } std::cout << std::endl; } } }; // 示例：生成一个简单的二维数据集（三个高斯分布） Dataset generateSampleData() { Dataset data; // 簇1 data.push_back({1.0, 1.0}); data.push_back({1.2, 0.9}); data.push_back({0.8, 1.1}); // 簇2 data.push_back({5.0, 5.0}); data.push_back({5.5, 4.5}); data.push_back({4.5, 5.5}); // 簇3 data.push_back({9.0, 1.0}); data.push_back({9.5, 1.5}); data.push_back({8.5, 0.5}); return data; } int main() { // 1. 准备数据 Dataset data = generateSampleData(); std::cout << "Dataset size: " << data.size() << " points." << std::endl; // 2. 创建聚类管理器并执行聚类 ClusterManager manager(data); manager.fit(); // 3. 获取3个簇的结果 int n_clusters = 3; std::vector<int> labels = manager.getClusterLabels(n_clusters); // 4. 输出结果 std::cout << "\nCluster labels for n_clusters=" << n_clusters << ":" << std::endl; for (size_t i = 0; i < data.size(); ++i) { std::cout << "Point " << i << " (" << data[i][0] << ", " << data[i][1] << ") -> Cluster " << labels[i] << std::endl; } // 5. 输出合并历史（树状图信息） std::cout << "\nMerge history (for dendrogram):" << std::endl; const auto& history = manager.getMergeHistory(); for (size_t i = 0; i < history.size(); ++i) { const auto& record = history[i]; std::cout << "Step " << i+1 << ": Merge cluster " << record.cluster1_id << "(size " << record.cluster1_size << ") and cluster " << record.cluster2_id << "(size " << record.cluster2_size << ") at distance " << record.distance << ", new cluster id: " << record.new_cluster_id << std::endl; } return 0; }

将上述代码保存为hierarchical_clustering.cpp，使用支持C++11及以上标准的编译器编译运行（例如g++ -std=c++11 -o hc hierarchical_clustering.cpp）。你会看到程序输出每个点被分配到的簇标签，以及详细的合并历史。对于我们的示例数据，三个分离的高斯点团应该被完美地分到三个不同的簇中。

6. 性能优化与高级话题探讨

我们上面实现的是一个清晰但效率不高的原型。当数据量超过几百个点时，O(n³)的时间复杂度（因为每次合并都要重新计算整个距离矩阵，而合并要进行n-1次）会使其变得非常慢。以下是几个关键的优化方向：

6.1 优化距离矩阵更新

每次合并后重新计算整个距离矩阵是最大的性能瓶颈。实际上，当簇A和簇B合并成新簇C后，只需要计算C与其他所有簇D的新距离，而其他簇对之间的距离保持不变。我们可以利用链接准则的递推公式来高效计算新距离。

以平均链接为例，假设簇C由A和B合并而成，那么C与另一个簇D的距离可以通过A、B、D的大小和原有距离计算得到：dist(C, D) = (|A| * dist(A, D) + |B| * dist(B, D)) / (|A| + |B|)

这样，更新距离矩阵的时间复杂度就从O(n²)降到了O(n)。我们需要在ClusterManager中维护每个簇的大小，并在合并时应用这个公式来更新distance_matrix_中与新簇相关的行和列，而不是全部重算。

6.2 使用优先队列查找最小距离

我们当前使用双重循环遍历上三角矩阵来查找最小距离，时间复杂度是O(n²)。可以使用一个最小堆（优先队列）来存储所有可能的簇对及其距离。每次合并后，只需要移除与被合并簇相关的所有对，并插入与新簇相关的所有新对。查找最小值的操作可以降到O(1)，但更新堆的操作需要仔细设计以保持正确性。一种常见的策略是使用“惰性删除”，即从堆顶弹出元素时检查该簇对是否仍然有效（两个簇都还存在）。

6.3 实现其他链接准则

我们只实现了平均链接。要支持单链接、全链接和沃德法，需要修改averageLinkageDistance函数，或者设计一个策略模式，将链接准则抽象成一个可插拔的类。例如：

class Linkage { public: virtual double compute(const Cluster& c1, const Cluster& c2, const Dataset& data) = 0; virtual ~Linkage() = default; }; class AverageLinkage : public Linkage { public: double compute(const Cluster& c1, const Cluster& c2, const Dataset& data) override { // ... 实现平均链接计算 } }; // 类似地实现 SingleLinkage, CompleteLinkage, WardsLinkage

然后在ClusterManager中持有一个Linkage的指针或引用，在计算距离时调用它。这样，用户可以在运行时指定不同的链接准则。

6.4 处理大规模数据与近似算法

对于海量数据（例如数十万以上），即使优化了矩阵更新，标准的层次聚类仍然不可行。此时需要考虑近似算法，例如：

• BIRCH算法：特别适用于大规模数值数据集，它通过构建一个聚类特征树（CF Tree）来对数据进行初步压缩和摘要，然后在摘要数据上执行层次聚类。
• 使用采样：先对数据进行随机采样，在样本上运行层次聚类，然后将剩余点分配到最近的簇中。
• 分布式计算：将数据分片，在不同节点上计算局部聚类，再合并结果。

7. 常见问题、调试技巧与实战建议

7.1 距离矩阵的对角线与初始化

注意：距离矩阵的对角线元素（自己到自己的距离）应该始终为0。在初始化时，我们显式或隐式地将其设为0。在查找最小距离时，必须跳过对角线，否则算法会在第一步就尝试将自己与自己合并。

7.2 浮点数精度问题

在比较距离大小时，直接使用==或<可能会因浮点数精度问题导致意外结果。建议使用一个极小的epsilon值进行比较：

if (std::abs(a - b) < 1e-12) { // 认为a等于b } if (a < b - 1e-12) { // 认为a小于b }

在查找最小距离时，如果出现多个相等的距离，需要定义打破平局的规则（例如，选择索引最小的对），以确保结果确定性。

7.3 簇ID的管理与并查集的使用

这是实现中最容易出错的地方之一。我们最初使用简单的遍历来更新标签，这在多次合并后会失效。务必使用并查集来管理簇的归属关系，如我们在改进版getClusterLabels中所做的那样。并查集的find操作（带路径压缩）和union操作非常高效，是处理这种动态连通性问题的标准工具。

7.4 内存使用与数据规模

我们的实现将整个距离矩阵存储在内存中，空间复杂度为O(n²)。对于10000个点，距离矩阵就需要存储大约1亿个双精度浮点数，占用约800MB内存。对于更大的数据，内存会成为瓶颈。此时，可以考虑：

• 使用float代替double以减少内存占用（会损失一些精度）。
• 使用稀疏矩阵格式存储距离（如果数据本身具有稀疏性，或者只存储最近邻距离）。
• 使用磁盘存储部分距离矩阵（但会极大增加I/O开销）。
• 直接转向更适合大数据的聚类算法（如BIRCH、Mini-Batch K-Means等）。

7.5 可视化树状图

获取merge_history_后，你可以使用Python的scipy.cluster.hierarchy和matplotlib来绘制树状图，这是分析层次聚类结果最直观的方式。你需要将合并历史转换为scipy需要的linkage matrix格式（一个(n-1)x4的矩阵）。每一行代表一次合并，包含：簇i索引、簇j索引、合并距离、新簇的大小（即cluster1_size + cluster2_size）。注意索引在scipy中通常是从0开始的。

7.6 确定最佳聚类数

层次聚类本身不给出“最佳”簇数。常用的确定方法有：

• 观察树状图：寻找树中合并距离突然增大的地方，在其下方切割。
• 肘部法则（Elbow Method）：绘制不同簇数下的簇内距离和（或类内方差），寻找拐点。
• 轮廓系数（Silhouette Score）：计算每个样本的轮廓系数并取平均，选择使平均轮廓系数最大的K值。

你可以用getClusterLabels获取不同K值下的聚类结果，然后计算这些指标来辅助决策。

从零开始用C++实现层次聚类，就像亲手搭建一座机械钟表，每一个齿轮（数据结构）的咬合，每一次指针（算法步骤）的跳动，都清晰可见。这个过程可能会让你对std::vector的erase操作保持警惕，对并查集的设计拍案叫绝，也对O(n³)的时间复杂度有了切肤之痛。但正是这些细节，构成了你对算法理解的深度。当你需要处理那些对性能有苛刻要求，或者环境限制无法引入庞大机器学习库的场景时，这份自己打造的、知其所以然的工具，会比任何黑盒调用都来得可靠。下一步，你可以尝试将链接准则抽象成接口，加入沃德法或更高效的距离矩阵更新，甚至用OpenMP为距离计算加上并行加速，让这个“轮子”跑得更快、更稳。