第11章 K-平均法クラスタリング

11.1 K-平均法アルゴリズム

K-平均法とは？

K-平均法 (K-Means) アルゴリズム [Mac67, Llo82] は、データセット内のクラスターを検出するための手法です。
実装がシンプルであり、広く利用されているクラスタリング手法の一つです。

このアルゴリズムでは、データポイントをK個のクラスターに分類し、それぞれのクラスターの中心を反復的に更新することで、
データを適切にグループ化します。

アルゴリズムの手順

K-平均法アルゴリズムは、以下の2つのステップを収束するまで繰り返すことで動作します。

各クラスターの中心を計算する
各ポイントを最も近い中心のクラスターに再割り当てする

この流れを視覚的に表したのが、以下の図11.1です。

図11.1: K-平均法の1回の反復

クラスターが2つ、ポイントが6つのK-平均法アルゴリズムの1回の反復の様子
各反復で、クラスターの中心 (＋) を計算し、各ポイントを最も近い中心のクラスターに再割り当てします。

K-平均法の詳細な処理

初期クラスターの設定
初めにK個のクラスターを適当に選び、それぞれの中心点 (Centroid) を決定します。
この初期化方法によって、結果が異なることがあります。
クラスターの中心を計算
各クラスターの中心を、そこに属する全ポイントの平均値として再計算します。
データポイントの再割り当て
各データポイントを、最も近いクラスターの中心に移動させます。
つまり、ユークリッド距離が最も小さいクラスターを探し、そのクラスターに所属させます。
収束判定
ポイントの割り当てが前回と変化しなくなるまで、2と3の処理を繰り返します。

クラスターが空になる場合の対処法

場合によっては、クラスターが空になってしまうことがあります。
これを防ぐため、以下のような対策を取ることが一般的です。

現在のクラスター内で、最も遠いポイントを空のクラスターに割り当てる
クラスターの初期化を工夫し、極端な偏りを防ぐ

並列化と最適化

K-平均法アルゴリズムは並列化が可能ですが、以下のような並列化の課題があります。

クラスター数が少ないと並列スラック (並列処理の効果) が不足する
ポイントの割り当てが不均衡になる可能性がある
メモリのキャッシュ効率が低下する (キャッシュラインに無関係なデータが含まれる)

そのため、並列化ではポイントの次元で並列処理を行うことが重要になります。

アルゴリズムの擬似コード

K-平均法アルゴリズムを簡単な擬似コードで表すと、以下のようになります。

合計
do {
    分割 // クラスターの中心を再計算
    再割り当て(マップ) // 各ポイントを最も近いクラスターへ移動
    合計(リダクション) // 各クラスターの新しい中心を計算
} while(ポイントを別のクラスターに移動);

11.3 K-平均法とTBB

TBB（Threading Building Blocks）でのK-平均法アルゴリズムは、Cilk Plusと同様の基本的な構造を持ちますが、いくつかの重要な違いがあります。

主な違い

配列表記文: 配列表記はシリアルループとして記述され、ベクトル化の目的ではなく表記の簡潔さが優先されます。これにより、パフォーマンスへの影響はありません。
ハイパーオブジェクトとTLS: ハイパーオブジェクトはスレッド・ローカル・ストレージ（TLS）に置き換えられ、ローカルビューからグローバルビューへのマージは明示的なループを使って行われます。
タイル形式（tiling）: 並列ループにおいて、タイル（ブロック）ごとにスレッド・ローカル・ビューの処理が行われ、パフォーマンスを最適化します。

スレッド・ローカル・ビューの実装

スレッドローカルビューとは並列処理において各スレッドが独自に持つデータのコピーのこと。

TBBでは、tbb::enumerable_thread_specificテンプレートを使用してスレッド・ローカル・ビューを実装します。これにより、スレッドは各ビューにアクセスすることができます。

初期クラスターの設定
初めにK個のクラスターを適当に選び、それぞれの中心点 (Centroid) を決定します。
この初期化方法によって、結果が異なることがあります。

ビューの作成方法には3つのアプローチがあります：

デフォルト構築 (enumerable_thread_specific<T> a;)
コピー構築 (enumerable_thread_specific<T> b(x);)
ファンクターによる構築 (enumerable_thread_specific<T> c(f);)

以下のコード例は、TBBでのビューの宣言と使用方法を示しています。

class view {
    view(const view& v);  // コピー構築を拒否
    void operator=(const view& v);  // 割り当てを拒否
public:
    sum_and_count* array;  // 動的配列（各クラスターのデータを保持）
    size_t change;  // 変更回数を記録する変数
    view(size_t k) : array(new sum_and_count[k]), change(0) {}  // コンストラクタ
    ~view() {delete[] array;}  // デストラクタ
};

typedef tbb::enumerable_thread_specific<view> tls_type;  // スレッドごとに独立した `view` のインスタンスを持つための型定義

ローカルビューのリダクション

ローカルビューのリダクション操作は、可換かつ結合的でなければならず、シリアルループで行うことが多いため、スケーラビリティに課題が生じることがあります。この制限は、TBBのparallel_reduceテンプレートを使うことで解消できます。

ローカルビューの合計をグローバルビューに累積するための関数は以下のように記述できます。

void reduce_local_counts_to_global_count(tls_type& tls, view& global) {
    global.change = 0;  // グローバルビューのchangeを初期化（リセット）
    for (auto i = tls.begin(); i != tls.end(); ++i) {  // tlsのbeginからendまでイテレーションを開始（各スレッドのローカルビューを処理）
        view& v = *i;  // 現在のスレッドのローカルviewを参照（各スレッドの結果にアクセス）
        global.change += i->change;  // 現在のスレッドのchange値をグローバルのchangeに加算（リダクション）
        v.change = 0;  // ローカルのchangeをリセット（次回の反復処理のために初期化）
    }
}

K-平均法の並列化

tbb::parallel_forを使って、並列でK-平均法の計算を実行します。これにより、各スレッドがサブ範囲ごとにポイントを処理します。

tbb::parallel_for(  // 並列でforループを実行
    tbb::blocked_range<size_t>(0, n),  // 処理範囲を0からnまで分割して並列化
    [...](tbb::blocked_range<size_t> r) {  // ラムダ関数を定義（並列実行される処理）
        view& v = tls.local();  // スレッドローカルのビューを取得（各スレッド専用のデータ）
        for (size_t i = r.begin(); i != r.end(); ++i) {  // 範囲r内でデータを処理
            ... // ポイントiを処理（各スレッドが担当する処理）
        }
    }
);

最適化とパフォーマンス

parallel_forの明示的なタイル形式は、内ループの最適化に有利である場合に使用します。スレッドローカルビューを効率的にルックアップし、反復処理のパフォーマンスを向上させます。

リスト11.7：最も近いクラスタのインデックスを検出する関数

指定したポイントに最も近いインデックスを検出するルーチン。これはリスト11.1の25行目の_sec_reduce_min_indと等価なシリアル表現です。

int reduce_min_ind( const point centroid[], size_t k, point value ) { 
    int min = -1;  // 最も近いクラスターのインデックスを保持する変数minを初期化(どのクラスターも選ばれていない)
    float mind = std::numeric_limits<float>::max();  // 最小距離を保持する変数mindを初期化（無限大に設定）
    for( int j = 0; j < k; ++j ) {  // すべてのクラスター（centroid）について反復処理
        float d = distance2(centroid[j], value);  // クラスターの重心centroid[j]とデータポイントvalueとの距離を計算
        if( d < mind ) {  // 計算した距離dが現在の最小距離mindより小さい場合
            mind = d;  // mindを新しい最小距離に更新
            min = j;  // minを現在のクラスターインデックスjに更新
        }
    }
    return min;  // 最も近いクラスターのインデックスminを返す
}

例：K-平均法の実装

TBBでのK-平均法アルゴリズムでは、以下のコードを使って反復的にクラスタを計算します。各反復では、ローカルビューをグローバルビューに累積し、新しいクラスタを再計算します。

void compute_k_means(size_t n, const point points[], size_t k, cluster_id id[], point centroid[]) {
    tls_type tls([&]{return k;}); 
    view global(k);

    // クラスタの初期化と合計の計算
    tbb::parallel_for(
        tbb::blocked_range<size_t>(0, n),
        [=, &tls, &global](tbb::blocked_range<size_t> r) {
            view& v = tls.local();
            for (size_t i = r.begin(); i != r.end(); ++i) {
                id[i] = i % k;
                v.array[id[i]].tally(points[i]);
            }
        }
    );

    // 反復処理（クラスタの更新）
    size_t change;
    do {
        // ローカル合計をグローバル合計に累積
        reduce_local_sums_to_global_sum(k, tls, global);

        // 空のクラスタの修復
        repair_empty_clusters(n, points, id, k, centroid, global.array);

        // セントロイドの計算
        for (size_t j = 0; j < k; ++j) {
            centroid[j] = global.array[j].mean();
            global.array[j].clear();
        }

        // 新しいクラスタの計算
        tbb::parallel_for(
            tbb::blocked_range<size_t>(0, n),
            [=, &tls, &global](tbb::blocked_range<size_t> r) {
                view& v = tls.local();
                for (size_t i = r.begin(); i != r.end(); ++i) {
                    cluster_id j = reduce_min_ind(centroid, k, points[i]);
                    if (j != id[i]) {
                        id[i] = j;
                        ++v.change;
                    }
                    v.array[j].tally(points[i]);
                }
            }
        );

        // ローカル変更のグローバル合計への累積
        reduce_local_counts_to_global_count(tls, global);
    } while (global.change != 0);
}

TBBを使用したK-平均法の実装は、Cilk Plusと同様のアーキテクチャを採用しつつ、スレッド・ローカル・ストレージ（TLS）を活用して並列化されています。これにより、大規模データセットに対する高効率なクラスタリングが可能となります。

11.4　まとめ

K-平均法は簡単に実装できるクラスター分析手法
並列化においては、ポイント単位の並列処理が有効
TBB 実装ではスレッド・ローカル・ストレージを活用し、リダクション処理を最適化