接地气可视化讲解聚类算法（一）：K-Medoids

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import make_blobs
from sklearn.cluster import KMeans
from sklearn_extra.cluster import KMedoids
from sklearn.metrics import pairwise_distances_argmin_min

# 生成主簇 + 多个离群点
X, _ = make_blobs(n_samples=300, centers=3, cluster_std=0.9, random_state=42)
outliers = np.random.uniform(low=15, high=25, size=(30, 2))
X_with_outliers = np.vstack([X, outliers])

# KMeans
kmeans = KMeans(n_clusters=3, random_state=42).fit(X_with_outliers)
kmeans_labels = kmeans.labels_
kmeans_centers = kmeans.cluster_centers_

# KMedoids
kmedoids = KMedoids(n_clusters=3, method='alternate', random_state=42).fit(X_with_outliers)
kmedoids_labels = kmedoids.labels_
kmedoids_centers = kmedoids.cluster_centers_

# 标签颜色对齐
mapping = pairwise_distances_argmin_min(kmedoids_centers, kmeans_centers)[0]
aligned_kmedoids_labels = np.array([mapping[label] for label in kmedoids_labels])

# 可视化（美观版）
fig, axes = plt.subplots(1, 2, figsize=(14, 6))
titles = ['K-Means 聚类（受异常点影响）', 'K-Medoids 聚类（更稳健）']
colors = ['#1f77b4', '#ff7f0e', '#2ca02c']

for ax, labels, centers, title, marker, center_color in zip(
    axes,
    [kmeans_labels, aligned_kmedoids_labels],
    [kmeans_centers, kmedoids_centers],
    titles,
    ['X', 'D'],
    ['red', 'blue']
):
    for i in range(3):
        cluster_points = X_with_outliers[labels == i]
        ax.scatter(cluster_points[:, 0], cluster_points[:, 1], s=30, color=colors[i], alpha=0.7)
    ax.scatter(centers[:, 0], centers[:, 1], c=center_color, marker=marker, s=250, edgecolor='white', linewidth=2, label='中心点')
    ax.set_title(title, fontsize=14)
    ax.set_xlim(-10, 30)
    ax.set_ylim(-10, 30)
    ax.grid(True, linestyle='--', alpha=0.5)
    ax.legend()

plt.tight_layout()
plt.show()

import random
import numpy as np

def manhattan_distance(a, b):
    return np.sum(np.abs(a - b))  # 你也可以改成欧氏距离

def total_cost(data, medoids, clusters):
    cost = 0
    for i, cluster in enumerate(clusters):
        for point in cluster:
            cost += manhattan_distance(data[point], data[medoids[i]])
    return cost

def assign_points(data, medoids):
    clusters = [[] for _ in medoids]
    for idx, point in enumerate(data):
        distances = [manhattan_distance(point, data[medoid]) for medoid in medoids]
        min_index = np.argmin(distances)
        clusters[min_index].append(idx)
    return clusters

def k_medoids(data, k, max_iter=100):
    n = len(data)
    medoids = random.sample(range(n), k)
    for _ in range(max_iter):
        clusters = assign_points(data, medoids)
        best_cost = total_cost(data, medoids, clusters)
        changed = False

        for i in range(k):
            for j in range(n):
                if j in medoids:
                    continue
                new_medoids = medoids[:]
                new_medoids[i] = j
                new_clusters = assign_points(data, new_medoids)
                new_cost = total_cost(data, new_medoids, new_clusters)
                if new_cost < best_cost:
                    medoids = new_medoids
                    best_cost = new_cost
                    changed = True
        if not changed:
            break
    return medoids, assign_points(data, medoids)

#include <iostream>
#include <vector>
#include <cstdlib>
#include <ctime>
#include <cmath>
#include <limits>
using namespace std;

double distance(const vector<double>& a, const vector<double>& b) {
    double sum = 0;
    for (size_t i = 0; i < a.size(); i++)
        sum += abs(a[i] - b[i]);
    return sum;
}

double total_cost(const vector<vector<double>>& data, const vector<int>& medoids, const vector<vector<int>>& clusters) {
    double cost = 0;
    for (size_t i = 0; i < clusters.size(); i++)
        for (int idx : clusters[i])
            cost += distance(data[idx], data[medoids[i]]);
    return cost;
}

vector<vector<int>> assign(const vector<vector<double>>& data, const vector<int>& medoids) {
    vector<vector<int>> clusters(medoids.size());
    for (int i = 0; i < data.size(); i++) {
        double min_dist = numeric_limits<double>::max();
        int best = 0;
        for (int j = 0; j < medoids.size(); j++) {
            double d = distance(data[i], data[medoids[j]]);
            if (d < min_dist) {
                min_dist = d;
                best = j;
            }
        }
        clusters[best].push_back(i);
    }
    return clusters;
}

pair<vector<int>, vector<vector<int>>> k_medoids(const vector<vector<double>>& data, int k) {
    srand(time(0));
    int n = data.size();
    vector<int> medoids;
    while (medoids.size() < k) {
        int r = rand() % n;
        if (find(medoids.begin(), medoids.end(), r) == medoids.end())
            medoids.push_back(r);
    }

    for (int iter = 0; iter < 100; iter++) {
        auto clusters = assign(data, medoids);
        double best_cost = total_cost(data, medoids, clusters);
        bool changed = false;

        for (int i = 0; i < k; i++) {
            for (int j = 0; j < n; j++) {
                if (find(medoids.begin(), medoids.end(), j) != medoids.end())
                    continue;
                vector<int> new_medoids = medoids;
                new_medoids[i] = j;
                auto new_clusters = assign(data, new_medoids);
                double new_cost = total_cost(data, new_medoids, new_clusters);
                if (new_cost < best_cost) {
                    medoids = new_medoids;
                    best_cost = new_cost;
                    changed = true;
                }
            }
        }
        if (!changed) break;
    }

    return {medoids, assign(data, medoids)};
}

use rand::seq::SliceRandom;
use rand::thread_rng;

fn manhattan(a: &[f64], b: &[f64]) -> f64 {
    a.iter().zip(b).map(|(x, y)| (x - y).abs()).sum()
}

fn assign(data: &[Vec<f64>], medoids: &[usize]) -> Vec<Vec<usize>> {
    let mut clusters = vec![vec![]; medoids.len()];
    for (i, point) in data.iter().enumerate() {
        let (min_idx, _) = medoids.iter()
            .enumerate()
            .map(|(j, &m)| (j, manhattan(point, &data[m])))
            .min_by(|a, b| a.1.partial_cmp(&b.1).unwrap())
            .unwrap();
        clusters[min_idx].push(i);
    }
    clusters
}

fn total_cost(data: &[Vec<f64>], medoids: &[usize], clusters: &[Vec<usize>]) -> f64 {
    clusters.iter().enumerate().map(|(i, cluster)| {
        cluster.iter().map(|&j| manhattan(&data[j], &data[medoids[i]])).sum::<f64>()
    }).sum()
}

fn k_medoids(data: &[Vec<f64>], k: usize) -> (Vec<usize>, Vec<Vec<usize>>) {
    let mut rng = thread_rng();
    let mut medoids: Vec<usize> = (0..data.len()).collect();
    medoids.shuffle(&mut rng);
    medoids.truncate(k);

    for _ in 0..100 {
        let clusters = assign(data, &medoids);
        let mut best_cost = total_cost(data, &medoids, &clusters);
        let mut updated = false;

        for i in 0..k {
            for j in 0..data.len() {
                if medoids.contains(&j) { continue; }
                let mut new_medoids = medoids.clone();
                new_medoids[i] = j;
                let new_clusters = assign(data, &new_medoids);
                let new_cost = total_cost(data, &new_medoids, &new_clusters);
                if new_cost < best_cost {
                    medoids = new_medoids;
                    best_cost = new_cost;
                    updated = true;
                }
            }
        }

        if !updated { break; }
    }

    (medoids, assign(data, &medoids))
}

def manhattan_distance(a, b):
    return np.sum(np.abs(a - b))

def assign_points(data, medoids):
    clusters = [[] for _ in medoids]
    ...

for i in range(k):
    for j in range(n):
        if j in medoids: continue

double distance(const vector<double>& a, const vector<double>& b)

vector<vector<int>> assign(const vector<vector<double>>& data, const vector<int>& medoids)

for (int i = 0; i < k; i++) {
    for (int j = 0; j < n; j++) {
        ...

fn manhattan(a: &[f64], b: &[f64]) -> f64 {
    a.iter().zip(b).map(|(x, y)| (x - y).abs()).sum()
}

fn assign(data: &[Vec<f64>], medoids: &[usize]) -> Vec<Vec<usize>> { ... }

medoids.shuffle(&mut rng);
medoids.truncate(k);

distance(a, b) = √((a1 - b1)² + (a2 - b2)² + ...)

let old = medoids[i];
medoids[i] = j; // 替换
...
medoids[i] = old; // 还原，避免 clone

use arrayvec::ArrayVec; // 限定容量，避免 heap 分配
let mut medoids: ArrayVec<usize, 10> = (0..n).collect();