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

11.1 アルゴリズム

k-平均法は、データセットをk個のクラスターに分類する代表的な反復型アルゴリズムである。各クラスターの「中心」を基準に、以下のプロセスを繰り返すことで、最適な分類を行う。

2つの主要ステップと並列化パターン

マップパターン
- 各データポイントに対し、最も距離が近い中心点を持つクラスターを特定し、そのクラスターに所属を更新する。
- 並列化の特性: 各データポイントの距離計算は他と完全に独立しているため、一斉に手分けして処理を行う「マップパターン」を適用できる。
リダクションパターン
- 各クラスターに属する全データポイントの座標の平均値を計算し、新しい中心点とする。
- 並列化の特性: 全データの情報を足し合わせて「集計・縮約」する必要があるため、複数のスレッドがバラバラに計算した結果を安全に合体させる「リダクションパターン」の手法が求められる。

フュージョンで処理を1つにまとめる

教科書通りの順番で作ると、「① 全データのグループ分け」がすべて終わった後に、もう一度最初からデータを読み直して「② 新しい中心の計算」をすることになる。これでは、重いデータを2回も読み込むことになり、オーバーヘッドが発生してしまう。

そこで、計算を限界まで速くするために、実際の並列化プログラムでは「グループ分け」と「リダクション」を1つのループ（do-while文）の中にフュージョンさせて実装する。こうすることで、データを読み込む回数を1回に減らし、圧倒的なスピードアップを実現できる。

// フュージョンによるループ構造のイメージ
do {
    // 1. 前回の結果から、各クラスターの仮の中心を決める
    // 2. 各データポイントに一番近い中心をみつける（マップ）
    // 3. みつけたと同時に、スレッドごとの手元の用紙に座標をメモする（リダクション）

} while ( データが別のクラスターに移動している間 ); // どのデータも動かなくなったら終了（収束）

11.2 k-平均法とCilk Plus

データ構造：`point` に関する最小限のルール（3つの仮定）

プログラムを並列で動かして、なおかつ計算を正しく・高速に行うためには、扱うデータ（point）が以下の3つの性質を必ず満たしている必要がある。

point() で「空」の状態を作れること
集計を始めるときに、最初に数値をリセットした「空の箱」を用意できる必要がある。
q += p で足し算加算が正しくできること
ポイント同士を安全に足し合わせて、合計座標を蓄積できる必要がある。
p / n で平均値が出せること
合計した座標をデータの個数nで割って、新しい「平均の中心座標」を計算できる必要がある。

リスト 11.1：Cilk Plus による k-平均法の実装

Cilk Plusの機能を使って、k-平均法の全体制御と並列化を行うメイン処理のコード。

void compute_k_means(point* points, size_t n, point* centroids, size_t k) {
    //初期クラスターの作成と割り当て
    cilk_for (size_t i = 0; i < n; ++i) {
        points[i].cluster_id = i % k;
    }

    //メインループ内の再割り当てステップ
    bool changed;
    do {
        changed = false;

        // ハイパーオブジェクト（リダクター配列）の宣言
        cilk::reducer_opadd<sum_and_count> sum[k];

        // 各データポイントの最近傍探索を並列実行
        cilk_for (size_t i = 0; i < n; ++i) {
            size_t id = reduce_min_ind(points[i], centroids, k);
            if (id != points[i].cluster_id) {
                points[i].cluster_id = id;
                changed = true; 
            }
            // スレッドごとの手元領域に座標と個数を安全に加算
            sum[id].ref().tally(points[i]);
        }

        // 新しい中心の計算
        for (size_t j = 0; j < k; ++j) {
            centroids[j] = sum[j].get_value().sum / sum[j].get_value().count;
        }

    } while (changed); 
}

このコードの役割

大量のデータポイントを、もっとも近い中心点（クラスター）にグループ分けし、全データの位置が変わらなくなるまで中心点を更新し続ける処理。

リスト 11.2：sum_and_count 構造体によるデータ定義

各クラスターの新しい中心（平均値）を計算するために、データを一時的に溜めておく「集計箱」となる構造体の定義コード。

struct sum_and_count {
    point sum;     // クラスターに属する全データポイントの「座標の合計」
    size_t count;  // クラスターに属するデータポイントの「個数の合計」

    // 最初に「ゼロ（初期状態）」を構築するためのコンストラクタ（最小限の仮定1を満たす）
    sum_and_count() : sum(point()), count(0) {}

    // データポイントを1つずつ集計領域に追加するためのメンバ関数
    void tally(const point& p) {
        sum += p;   // 座標を加算
        count++;    // 個数をカウントアップ
    }

    // スレッドごとの手元領域を自動統合（マージ）するための += 演算子
    // （最小限の仮定2を満たす）
    sum_and_count& operator+=(const sum_and_count& right) {
        sum += right.sum;
        count += right.count;
        return *this;
    }
};

このコードの役割

*　各クラスターに属するデータの「座標の合計」と「データの個数」を常にセットで記録・更新するための仕組み。

11.2.1 ハイパーオブジェクト

複数のスレッドが同時に同じクラスターの合計値を更新しようとすると、計算結果がぐちゃぐちゃに破壊されるデータ競合が発生する。Cilk Plusではこれを防ぐため、リダクターを使用する。

リダクターの宣言と並列ループ

// クラスター数 k 個分のリダクター配列を宣言
cilk::reducer_opadd<sum_and_count> sum[k]; 

cilk_for (size_t i = 0; i < n; ++i) {
    // 各データポイントに最も近いクラスターのID（0 ～ k-1）を計算
    size_t id = reduce_min_ind(points[i], centroids, k);

    // 計算したクラスターIDの集計領域に対して、安全にデータを加算
    sum[id].ref().tally(points[i]); 
}

11.2.1 のポイント

cilk::reducer_opadd<sum_and_count> sum[k];
- 何をしているか
  - クラスターの数kの分だけ、並列処理に対応した特別な「リダクター」の配列を作っている。
- なぜ必要なのか
  - 普通の配列だと、複数のスレッドが同じ場所に同時にデータを書き込んだときに、値が上書きされて消えてしまう「データ競合」が発生する。
  - このリダクターを使うことで、Cilk Plusがスレッドごとに独立したビューをメモリ上に自動で用意してくれる。その結果、スレッド同士がお互いに邪魔をすることなく、一斉に計算を進めることができる。

11.3 Intel TBB による実装と最適化

Intel TBB（Threading Building Blocks）を用いた実装では、コンパイラ拡張機能（cilk_for など）を使わずに、ライブラリのテンプレートクラスや構造体を利用して並列化を行う。そのため、メモリ管理や最終的なデータの集計プロセスを明示的に記述する必要がある。

スレッド間の干渉防止（TLSの定義）

リスト 11.3：Intel TBB による実装（TLSの定義）

データ競合を回避するために tbb::enumerable_thread_specific を用いて、各スレッド専用のワークスペースであるスレッド・ローカル・ストレージ（TLS）を構造体内部に定義する。

struct tbb_k_means {
    point* points;     // 全データポイントへのポインタ
    point* centroids;  // 各クラスターの現在の中心点の配列（points, centroidsは、全スレッドから参照される共有データだが、ループ内では読み取り専用として扱うため、データレースは発生しない）
    size_t k;          // 分類するクラスターの数
    bool changed;      // データポイントが移動したかを記録するフラグ

    // 各スレッドが独立してクラスターごとの座標の合計とデータ数をカウントする必要があるため、各スレッドに対して、クラスター数k分の要素を持つ std::vector をTLSとして自動分配
    tbb::enumerable_thread_specific<std::vector<sum_and_count>> local_sums;

    // コンストラクタ
    tbb_k_means(point* points_, point* centroids_, size_t k_)
        : points(points_), centroids(centroids_), k(k_), changed(false),
          local_sums(std::vector<sum_and_count>(k_)) {} // 最初にクラスター数分(k)の部屋を用意する
};

このコードの役割

何をしているか
- Intel TBBで並列ループ処理（parallel_for）を実行するために必要なデータや、スレッドごとの作業領域を一つにまとめた「タスク管理用の構造体」を定義している。
- 6行目で、TBB特有の安全な集計領域である local_sums（TLS）をメンバー変数として持たせている。

リスト11.4 Intel TBB におけるスレッド・ローカル・ビューの宣言

Intel TBBでk-平均法を動かす際、Cilk Plusのハイパーオブジェクトのような仕組みをライブラリの機能で構築するため、各スレッドが個別に持つ「ビュー」を明示的に定義する。

リスト 11.4：TBBにおけるスレッド・ローカル・ビューのtls_type型の宣言

データ競合を完全に回避するため、各スレッドが持つべき「集計用配列（sum_and_count）」と「移動検知用フラグ（change）」を1つにまとめた view クラス、およびそのTLS型の定義コード。

class view {
    view( const view& v );                    // 本来はクローンを作るための関数を、あえてprivateの場所に隠すことで、このviewがプログラムの途中で勝手にコピーされるのを絶対に禁止する。
    void operator=( const view& v );         // // 同様に、viewの中身が別のviewのデータで勝手に上書きされるのを絶対に禁止する。

public:
    sum_and_count* array;                    // クラスター数k個分の集計配列へのポインタ(記録メモ)
    size_t change;                           // スレッド内でのデータ移動回数をカウントするフラグ

    // コンストラクタ：新しくスレッド用メモが作られるときに、k個sum_and_countを用意しchangeをリセット
    view( size_t k ) : array(new sum_and_count[k]), change(0) {}

    // デストラクタ：処理が終わってスレッド用メモが破棄されるときに、メモリを綺麗に解放
    ~view() {delete[] array;}
};

// 定義した view クラスを、TLSとして扱うための型定義
typedef tbb::enumerable_thread_specific<view> tls_type;

このコードの役割

何をしているか
- 各スレッドが並列ループの中で他人に邪魔されずに計算結果を書き込めるように、スレッド専用の「記録メモ（array）」と「移動カウンター（change）」をセットにした view クラスを定義している。
なぜこの構造が必要か
- Cilk Plusでは、移動フラグの集計もリダクター（cilk::reducer_opadd<size_t> change）が勝手にやってくれましたが、TBBにはそれがない。
- そのため、「座標の合計（sum_and_count）」だけでなく、「そのスレッド内でデータが何回動いたか（change）」も、スレッドごとに完全に隔離した状態で一緒にメモして持ち回る必要があるため、このような構造が必要になる。

リスト 11.5：ローカルビューの変更の検出

各スレッドが個別にカウントしていた「データポイントの移動回数（change）」をすべて足し合わせて、全体の移動合計を算出する関数。

void reduce_local_counts_to_global_count( tls_type& tls, view& global )　//　tls:各スレッドに完全に独立して割り当てられるメモ
// global:並列ループがすべて終わった後に、最終的な答えを書き込むメモ
 {
    global.change = 0; // まず全体の移動カウンターをゼロに初期化

    // TLS領域を1つずつ順番に巡回するループ
    for ( auto i=tls.begin(); i!=tls.end(); ++i ) {
        view& v = *i;               // 現在スキャンしているスレッドのビューを取得
        global.change += v.change;  // スレッド内での移動回数を全体の合計に足し算
        v.change = 0;               // 次の反復に備えて、スレッド内のカウンターを0にクリア
    }
}

リスト11.5の役割

何をしているか
- 各スレッドが自分の担当タイル内で検知した「ポイントが動いた回数（v.change）」をループで全て回収し、全体のメモ（global.change）に足し算している。
- 回収し終わったスレッドの手元カウンターは、次の周回のために 0 にリセットしている。
なぜこの構造が必要か
- Cilk Plusのときは、ポイントが1つでも動いたかどうかの判定（change フラグの集計）をシステムが自動で処理してくれましたが、TBBではスレッドごとにバラバラのメモとして残ったままになる。
- そのため、まだ動いているポイントがあるかを正確に判定するために、このように明示的なループを回して手動で集計する必要がある。

リスト 11.6：ローカルビューからグローバル合計への累積（座標と個数の集計）

各スレッドの手元メモに溜まった「座標の合計（sum）」と「ポイントの個数（count）」を、グローバル配列へ安全に合体させる関数。

void reduce_local_sums_to_global_sum( size_t k, tls_type& tls, view& global ) {
    // TLS領域を1つずつ順番に巡回するループ
    for ( auto i=tls.begin(); i!=tls.end(); ++i ) {
        view& v = *i; // 現在スキャンしているスレッドのビューを取得

        // クラスターの数k個分だけループを回して、配列の中身を1つずつ足し合わせる
        for ( size_t j=0; j<k; ++j ) {
            global.array[j] += v.array[j]; // スレッドの部分合計を、全体の同じクラスター位置に加算
            v.array[j].clear();            // 次の周回に備えて、スレッドの手元配列をクリア
        }
    }
}

リスト11.6の役割

何をしているか
- 各スレッドが別々に計算したクラスターごとの「座標の足し算結果」と「カウントした個数」を、2重ループを使って全体の共通配列（global.array）にすべて集約している。
- 集計が終わったスレッドの配列は、次の周回のために clear() で綺麗にリセットしている。
なぜこの構造が必要か
- 並列処理中は、スレッド同士の衝突を防ぐために、わざとデータをバラバラのTLSに保存させていた。
- ループが終わった後は、それらを1つにまとめないと新しい中心座標が計算できないため、構造を作成する。

リスト 11.7：最も近い中心のインデックスを検出するルーチン

指定されたデータポイントから最も距離が近い中心点（クラスター）のIDを、愚直なループ（シリアル処理）で探索して返す補助関数。

int reduce_min_ind( const point centroid[], size_t k, point value )// centroid[]：すべてのクラスターの中心座標のリスト、k：クラスターの総数、value：現在のデータ1個の座標
 {
    int min = -1;                // 最も近いクラスターのIDを保存する変数をどこにも所属していない状態を表す -1にセット

    float mind = std::numeric_limits<float>::max();  // 最小距離を保持する変数

    // クラスター数kの分だけループを回し、一番近い中心点を探す
    for( int j = 0; j < k; ++j ) {
        float d = distance2(centroid[j], value);     // 中心点とデータポイントの間の距離を計算
        if( d < mind ) {                             // これまでの最小距離よりも短ければ更新
            mind = d;         // 最小距離を上書き
            min = j;         // そのときの一番近いクラスターIDを記録
        }
    }
    return min;                // 見つかった一番近いクラスターIDを返す
}

リスト11.7の役割

何をしているか
- 1つのデータポイントに注目し、すべてのクラスター中心（centroid）との距離を順番に計算して、「自分が今、どのグループに一番近いか」のインデックスを特定して返している。
なぜこの構造が必要か
- Cilk Plusのときは、コンパイラ専用の特殊な並列リダクション表現（__sec_reduce_min_ind）を使って1行で書けましたが、Intel TBBにはそのような専用のコンパイラ構文がない。
- そのため、通常のC++コードとして同じ動きをするシリアル（逐次）ループを明示的に自作して、TBBの並列ループ内から呼び出す必要がある。

リスト 11.8：TBB による k-平均法クラスタリング

これまで定義してきたスレッドローカル領域（tls）やマージ関数（リスト11.5, 11.6）をすべて組み合わせた、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;});                     // スレッドごとにk個の部屋を持つTLSを初期化
    view global(k);                                   // 最終集計用のグローバル領域（サイズ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]);      // スレッド手元の集計用紙に座標と個数を加算
            }
        }
    );

    // idが変更されなくなるまでメインループを回す（収束判定）
    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 ) {
                    // 最も近い新しい中心のIDを検索（リスト11.7の呼び出し）
                    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]);      // 新しいグループの手元用紙にデータを加算
                }
            }
        );

        // ローカルカウントからグローバルカウントへのレデュース（リスト11.5の呼び出し）
        reduce_local_counts_to_global_count( tls, global );

    } while( global.change != 0 );                     // どのスレッドでも移動が発生しなくなったら終了
}

リスト11.8の役割

何をしているか
- これまで個別に作ってきた「スレッド専用メモの準備（tls）」「最短距離の探索（reduce_min_ind）」「データの回収（reduce_local_...）」を1つの流れに繋ぎ合わせ、データが全く動かなくなるまで（収束するまで）ループを回してクラスタリングを完成させている。
なぜこの構造が必要か
- TBBではコンパイラが自動で裏側の同期をやってくれないため、処理のタイムラインを明示的に制御する必要がある。
- 「並列で手分けして計算する（parallel_for）」 $\rightarrow$ 「ループが終わったら一度みんなのメモを一箇所に集める（reduce）」 $\rightarrow$ 「集まったデータを元に中心点を更新する（centroid[j] = ...）」という、並列処理の正しいフェーズ管理を破綻なく実行させるために、この統括関数が必要になる

11.4 まとめ

k-平均法クラスタリングを並列化するにあたり、Cilk PlusとIntel TBBのどちらを使うかで、コードの書き方や裏側の仕組みが大きく変わる。

Cilk Plus と Intel TBB の決定的な違い

Cilk Plusは「すべて自動おまかせ」
- コンパイラの特別な機能を使っているため、コードが非常にシンプル。
- リダクターを宣言しておけば、並列ループが終わった瞬間に、各スレッドがバラバラにメモしていた結果をシステムが裏側で自動的に合体してくれる。
Intel TBBは「手動で明示的に書く」
- TBBは純粋なC++のライブラリなので、裏側で勝手にやってくれない。
- スレッドごとにTLSを自作し、並列ループが終わった後には、それらのメモを1つずつ手作業で回収して足し合わせるコードを自分で書く必要がない。

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

11.1 アルゴリズム

2つの主要ステップと並列化パターン

フュージョンで処理を1つにまとめる

11.2 k-平均法とCilk Plus

データ構造：point に関する最小限のルール（3つの仮定）

リスト 11.1：Cilk Plus による k-平均法の実装

このコードの役割

リスト 11.2：sum_and_count 構造体によるデータ定義

このコードの役割

11.2.1 ハイパーオブジェクト

リダクターの宣言と並列ループ

11.2.1 のポイント

11.3 Intel TBB による実装と最適化

スレッド間の干渉防止（TLSの定義）

リスト 11.3：Intel TBB による実装（TLSの定義）

このコードの役割

リスト11.4 Intel TBB におけるスレッド・ローカル・ビューの宣言

リスト 11.4：TBBにおけるスレッド・ローカル・ビューのtls_type型の宣言

このコードの役割

リスト 11.5：ローカルビューの変更の検出

リスト11.5の役割

リスト 11.6：ローカルビューからグローバル合計への累積（座標と個数の集計）

リスト11.6の役割

リスト 11.7：最も近い中心のインデックスを検出するルーチン

リスト11.7の役割

リスト 11.8：TBB による k-平均法クラスタリング

リスト11.8の役割

11.4 まとめ

Cilk Plus と Intel TBB の決定的な違い

データ構造：`point` に関する最小限のルール（3つの仮定）