第 14 章: K-means によるクラスタリング¶
14.1 はじめに¶
この章の クラスタリング は、正解ラベルの無いデータを、似たものどうしのかたまりに分ける教師なし学習です。使うのは K-means で、あらかじめクラスタの数 \(k\) を決め、「各点を最も近い中心に割り当てる」「クラスタごとの平均を新しい中心にする」を中心が動かなくなるまで繰り返します。
なでしこ3 には K-means のライブラリが無いので、Haskell 版と同じく 自作したものがそのまま最終実装になります（ADR 015）。突き合わせる相手は、Haskell 版など java.util.Random と同じ乱数で初期中心を選ぶ言語版の数値と、自作のデータで動かす Python の参照実装です。
なでしこ3 版では、次の 3 点に注目してください。
- gonako は 1 回の実行で 2 億命令を超えると止まる。Haskell 版と同じ書き方で K-means を書くと、報告を作り終える前に「実行した命令が多すぎます」で止まりました。この章の大半は、命令の数を数えて減らす話です
- 命令を減らしても、結果は 1 ビットも変えない。距離を測る点を減らす仕組み（Hamerly の方法）を入れたあとも、割り当ては「すべての中心と前から順に比べる」素朴な繰り返しと同じでなければなりません。これを、同点がわざと多くなる格子点のテストで確かめます
- Go で書かれた命令は 1 回で数えられる。
表列合計は 440 行を足しても命令 1 つ分です。足す順が自作の繰り返しと同じであることをソースで確かめてから使います
14.2 K-means の仕組み¶
- クラスタ数 \(k\) を決め、初期中心を \(k\) 個選ぶ
- 割り当て: 各点を、最も近い中心のクラスタに入れる
- 更新: クラスタごとに、属する点の平均を新しい中心にする
- 中心が動かなくなるまで 2 と 3 を繰り返す
「最も近い」は距離の 2 乗で測ります。K-means は初期中心に依存し、悪い置き方をすると局所解に落ちるので、初期中心を 10 通り試して SSE（各点と所属する中心の距離の 2 乗の合計）が最小のものを採ります。ここまでは Haskell 版と同じです。
14.3 題材とデータ¶
卸売業者の顧客ごとの年間支出額（Wholesale.csv、440 件）を使います。Channel と Region は区分を表す番号なので除き、残る 6 列（Fresh・Milk・Grocery・Frozen・Detergents_Paper・Delicassen）を、列ごとに平均 0・標準偏差 1（件数で割る標準偏差）にそろえてからクラスタリングします。
振る舞いはすべて Haskell 版にそろえます。
| 項目 | 決めごと |
|---|---|
| 初期中心 | シード付きの乱数（第 2 章の自作の線形合同法）で点の番号を並べ替え、先頭から \(k\) 個 |
| 試す回数 | シードを 0 から 1 ずつずらして 10 通り。SSE が同じなら先に試したほう |
| 打ち切り | 中心がすべての値で変わらなくなるか、更新が 300 回に達したら |
| 同点 | 距離が同じなら番号の小さい中心 |
| 空のクラスタ | 点が 1 つも割り当てられなければ、前の中心を残す |
| 足し算の順 | どれも 0 から始めて前から順に足す |
14.4 TODO リストの作成¶
TODO リスト:
- 支出額の列だけを読み込む
- 列ごとに標準化する
- 2 点間の距離の 2 乗を求める
- 各点を最も近い中心に割り当てる
- クラスタごとの平均を新しい中心にする
- SSE を計算する
- 中心が変わらなくなるまで繰り返す
- 初期中心をシードで選び、何通りか試して SSE が最小の結果を採る
- クラスタ数ごとの SSE を並べる（エルボー法）
- クラスタごとの件数と平均支出額をまとめる
ほかの章と同じく、名前は章の題材を含むものにしました（クラスタ割当・卸売標準化 など）。第 9 章の標準化はほかのエージェントが並行して作っているので、ここでは第 14 章専用に書きます。
14.5 支出額の列を読み込む¶
卸売二件 = 「Channel,Region,Fresh,Milk{改行}1,3,10,100{改行}2,3,30,300{改行}」
二件支出額 = (卸売二件を表読込)の卸売支出額抽出
「Channel と Region を除いた列だけを読む」で(二件支出額["列名一覧"])と[「Fresh」, 「Milk」]が一致確認
空欄支出額 = 関数()
(「Channel,Region,Fresh,Milk{改行}1,3,10,{改行}」を表読込)の卸売支出額抽出
ここまで
「空欄があれば失敗する」で空欄支出額を「Milk を数値として読めません: 」でエラー確認
読み込みは第 1 章の 表読込 と 整数取得 をそのまま使います。整数取得 は正規表現で形を確かめるので、空欄は Haskell 版と同じメッセージで止まります。この章の題材に欠損値は無いので、補完はしません。
●(対象表の)卸売支出額抽出とは
列名一覧 = []
全列 = 対象表["列名一覧"]
Iを 0から (全列の要素数) - 1まで繰り返す
もし、(卸売区分列の(全列[I])を配列検索) < 0ならば
列名一覧に全列[I]を配列追加
ここまで
ここまで
行一覧 = 対象表["行一覧"]
支出額一覧 = []
Iを 0から (行一覧の要素数) - 1まで繰り返す
支出額 = {}
Jを 0から (列名一覧の要素数) - 1まで繰り返す
支出額[列名一覧[J]] = 行一覧[I]の(列名一覧[J])を整数取得
ここまで
支出額一覧に支出額を配列追加
ここまで
{"列名一覧": 列名一覧, "支出額一覧": 支出額一覧}で戻る
ここまで
14.6 列ごとに標準化する¶
「列ごとに平均零・標準偏差一にそろえる」で((二件支出額["支出額一覧"])を(二件支出額["列名一覧"])で卸売標準化)と[[-1, -1], [1, 1]]が一致確認
「値がすべて同じ列は標準偏差を一として扱う」で([{"Fresh": 5}, {"Fresh": 5}]を[「Fresh」]で卸売標準化)と[[0], [0]]が一致確認
2 件なら標準化した値は必ず −1 と 1 になるので、誤差なく比べられます。実装は、平均・分散とも 0 から前の行へ順に足し、件数で割ります。Haskell 版の fitStandardizer と同じ順です。
14.7 各点を最も近い中心に割り当てる¶
「各点を最も近い中心に割り当てる」で([[0], [10], [4]]を[[0], [10]]でクラスタ割当)と[0, 1, 0]が一致確認
「距離が同じなら先に並ぶ中心を選ぶ」で([[5]]を[[0], [10]]でクラスタ割当)と[0]が一致確認
●(点一覧を中心一覧で)クラスタ割当とは
ラベル一覧 = []
Iを 0から (点一覧の要素数) - 1まで繰り返す
点 = 点一覧[I]
最近 = 0
最短 = 点と(中心一覧[0])のクラスタ二乗距離
Kを 1から (中心一覧の要素数) - 1まで繰り返す
距離 = 点と(中心一覧[K])のクラスタ二乗距離
もし、距離 < 最短ならば
最近 = K
最短 = 距離
ここまで
ここまで
ラベル一覧に最近を配列追加
ここまで
ラベル一覧で戻る
ここまで
厳密な < で比べるので、同じ距離なら先に並ぶ中心が残ります。この素朴な割り当ては、14.14 節で速くした割り当てが「同じ結果を出しているか」を確かめる基準として、最後まで残します。
14.8 中心を更新する¶
「クラスタごとの平均を新しい中心にする」で([[0], [2], [10]]と[0, 0, 1]を[[0], [0]]でクラスタ中心更新)と[[1], [10]]が一致確認
「点が割り当てられなかったクラスタは前の中心を残す」で([[0], [2]]と[0, 0]を[[0], [7]]でクラスタ中心更新)と[[1], [7]]が一致確認
最初は列ごとの和をなでしこ3 の繰り返しで足していましたが、14.14 節の理由で、Go で書かれた命令 表列合計 に置き換えました。
●(点一覧とラベル一覧を前中心一覧で)クラスタ中心更新とは
所属一覧 = []
Kを 0から (前中心一覧の要素数) - 1まで繰り返す
所属一覧に[]を配列追加
ここまで
Iを 0から (点一覧の要素数) - 1まで繰り返す
所属一覧[ラベル一覧[I]]に点一覧[I]を配列追加
ここまで
新中心一覧 = []
Kを 0から (前中心一覧の要素数) - 1まで繰り返す
所属 = 所属一覧[K]
件数 = 所属の要素数
もし、件数 = 0ならば
新中心一覧に前中心一覧[K]を配列追加
違えば
平均 = []
Jを 0から (所属[0]の要素数) - 1まで繰り返す
平均に((所属のJを表列合計) ÷ 件数)を配列追加
ここまで
新中心一覧に平均を配列追加
ここまで
ここまで
新中心一覧で戻る
ここまで
置き換える前に、表列合計 の足す順を処理系のソース（internal/stdlib/array.go）で確かめました。合計を 0.0 から始め、0 行目から順に total += で足しています。点が現れた順に 0 から足す、自作の繰り返しと同じ順なので、中心の値は最後の桁まで変わりません。
14.9 SSE を計算する¶
●(点一覧とラベル一覧を中心一覧で)クラスタ誤差平方和とは
合算 = 0
Iを 0から (点一覧の要素数) - 1まで繰り返す
合算 = 合算 + (点一覧[I]と(中心一覧[ラベル一覧[I]])のクラスタ二乗距離)
ここまで
合算で戻る
ここまで
440 個の足し算なので、順序が変わると下の桁が動きます。ここも 0 から前の点へ順に足します。
14.10 中心が変わらなくなるまで繰り返す¶
中心がすべての値で等しくなったら止めます。
●(甲一覧と乙一覧が)クラスタ中心同一判定とは
もし、(甲一覧の要素数) ≠ (乙一覧の要素数)ならば
偽で戻る
ここまで
Kを 0から (甲一覧の要素数) - 1まで繰り返す
甲 = 甲一覧[K]
乙 = 乙一覧[K]
Jを 0から (甲の要素数) - 1まで繰り返す
もし、甲[J] ≠ 乙[J]ならば
偽で戻る
ここまで
ここまで
ここまで
真で戻る
ここまで
浮動小数点数の等値比較ですが、収束すると同じ点の同じ順の平均を計算するので、値は「近い」ではなく「同じ」になります（Haskell 版の next == current と同じ考え方です）。比べるのは数値どうしなので、ゆるい比較の ≠ でも困りません。
二塊結果 = 二塊を[[0, 0], [10, 10]]と300でクラスタ中心学習
「離れた二つのかたまりを分ける（ラベル）」で(二塊結果["ラベル一覧"])と[0, 0, 1, 1]が一致確認
「離れた二つのかたまりを分ける（中心）」で(二塊結果["中心一覧"])と[[0.5, 0.5], [10.5, 10.5]]が一致確認
「離れた二つのかたまりを分ける（SSE）」で(二塊結果["SSE"])と2が一致確認
「中心の更新を計算した回数を数える（収束を確かめた回を含む）」で(二塊結果["反復回数"])と2が一致確認
打切結果 = 二塊を[[0, 0], [10, 10]]と0でクラスタ中心学習
「更新の回数の上限に達したらそこで打ち切る」で(打切結果["中心一覧"])と[[0, 0], [10, 10]]が一致確認
結果の辞書には、Haskell 版に無い 反復回数（中心の更新を計算した回数。「変わらなかった」と確かめた最後の 1 回を含む）を足しました。14.14 節で命令の数を見積もるのに使います。
14.11 局所解を再現する¶
Haskell 版と同じ、1 次元に並んだ 3 つのかたまりで局所解を再現します。
悪い結果 = 三塊を[[0], [1], [10]]と300でクラスタ中心学習
良い結果 = 三塊を[[0.5], [10.5], [20.5]]と300でクラスタ中心学習
「初期中心の置き方で結果が変わる（良い置き方）」で(良い結果["SSE"])と1.5が一致確認
「初期中心の置き方で結果が変わる（局所解の SSE）」で(悪い結果["SSE"])と101が一致確認
「初期中心の置き方で結果が変わる（局所解のラベル）」で(悪い結果["ラベル一覧"])と[0, 1, 2, 2, 2, 2]が一致確認
手前のかたまりを 2 つの中心が分け合うと、奥の 4 点を 1 つの中心が抱えたまま動かなくなり、SSE は最良の 1.5 に対して 101 になります。
14.12 初期中心と複数回の試行¶
「シード付きの乱数で並べ替えて先頭から選ぶ」で([[0], [1], [2], [3]]の2を0でクラスタ初期中心選択)と[[3], [0]]が一致確認
「点より多くの中心は選べない」で中心過多選択を「点は 1 件しかないので 2 個の中心を選べません」でエラー確認
「SSE が同じなら先に試した候補を残す」で((二塊を[[[0, 0], [10, 10]], [[10, 10], [0, 0]]]でクラスタ最良選択)["ラベル一覧"])と[0, 0, 1, 1]が一致確認
[[3], [0]] は Haskell 版のテストと同じ値です。第 2 章の 並べ替え が java.util.Random と同じ並びを出すので、ほかの言語版と同じ点が初期中心になります。
エルボー法では、クラスタ数ごとに 10 通り試した最良の結果を [クラスタ数, 結果] の組で並べます。
●(点一覧をクラスタ数一覧とシードと試行回数で)クラスタ数別学習とは
結果 = []
Iを 0から (クラスタ数一覧の要素数) - 1まで繰り返す
クラスタ数 = クラスタ数一覧[I]
結果に[クラスタ数, (点一覧のクラスタ数をシードと試行回数でクラスタ複数試行)]を配列追加
ここまで
結果で戻る
ここまで
Haskell 版の報告は、SSE を求めたあとにクラスタ数 5 の学習をもう一度行います。同じシードから同じ計算をするので結果は同じです。なでしこ3 版では 14.14 節の理由で、エルボー法で求めたクラスタ数 5 の結果を使い回します。
14.13 クラスタごとの特徴をまとめる¶
平均は元の単位（金額）で出し、件数の多い順（同数ならクラスタ番号の小さい順）に並べます。並べ替えは、番号順に並んだものを「件数が真に多いときだけ前へ動かす」挿入ソートにしたので、同数の順が保たれます。
「件数の多い順に並べ、元の単位で平均を出す」で要約一覧と[{"クラスタ": 1, "件数": 2, "平均一覧": [20, 200]}, {"クラスタ": 0, "件数": 1, "平均一覧": [30, 300]}]が一致確認
「件数が同じならクラスタ番号の小さい順に並べる」で(((三件支出額["支出額一覧"])を(三件支出額["列名一覧"])と[2, 0, 1]でクラスタ特徴要約)[0]["クラスタ"])と0が一致確認
平均支出額は整数に丸めて表示します。Haskell 版は printf "%.0f" が 0.5 を偶数側に丸めるため、ちょうど 34708.5 が出たところでほかの言語版と 1 だけ食い違い、「0.5 を足して切り下げる」に書き直しました。なでしこ3 の 四捨五入 は 34708.5 を 34709 にしましたが、振る舞いを Haskell 版とそろえておくため、同じ式で書きます。
●(値の)支出額四捨五入とは
(値 + 0.5)をFLOORで戻る
ここまで
「ちょうど零点五は大きいほうへ丸める」で(34708.5の支出額四捨五入)と34709が一致確認
14.14 実データでクラスタリングする¶
2 億命令で止まる¶
ここまでを Haskell 版と同じ形で書き、実データと同じ形の架空のデータ（440 行・6 列・欠損なし。配布データではない。この節の後半を参照）で報告を作らせると、約 3.2 秒後に止まりました。
[実行時エラー]run14.nako3(35行目): 実行した命令が多すぎます。終わらない繰り返しになっていませんか。
処理系のソース（internal/vm/vm.go）を読むと、実行の既定の設定に MaxInstructions: 200_000_000 があり、1 回の実行で 2 億命令を超えると止まる作りでした。gonako --help にこの上限を変える引数はありません。第 1〜6 章のテスト（学習データの無い環境）は、この上限に当たっていません。
命令の数を数える¶
減らすには、まず数える必要があります。処理系を手元で作り直し、実行の最後に命令の数と関数ごとの内訳を表示するようにしました（計測のためだけの改造で、リポジトリには入れていません）。
最初の書き方（距離も中心の更新もなでしこ3 の繰り返しで書き、毎回すべての点とすべての中心の距離を測る）では、報告に要る 100 回の学習だけで 4 億 7372 万命令でした。上限の 2.4 倍です。更新の回数は 100 回の学習で合わせて 1566 回で、Python の参照実装と同じでした。
命令の数を小さな計測で分けると、なでしこ3 のどこが高いかが分かります。
| 書き方（6 列の点の距離の 2 乗を 1 回求める） | 命令の数 |
|---|---|
Iを 0から (甲の要素数) - 1まで繰り返す で 1 列ずつ足す |
81 |
同じ繰り返しで (甲[I] - 乙[I]) ^ 2 と書く |
87 |
甲を反復 で 対象 と 乙[対象キー] を使う |
172 |
甲[I + 1] のように添字を計算して 6 列ずつ展開する |
84 |
甲[0] 〜 甲[5] と定数の添字で 6 列を展開する |
42 |
（参考）配列合計 で 6 個を足す |
3 |
繰り返しの制御に、足し算そのものと同じくらいの命令がかかっています。 添字に計算を入れると、展開しても減りません。そして Go で書かれた命令（配列合計・表列合計）は、中で何個足しても数命令です。
距離を測る点を減らす¶
命令の大半は距離の計算でした。そこで K-means の結果を変えずに、測る距離の数を減らします。使うのは Hamerly の方法 です。点ごとに次の 2 つを持ち回ります。
- 上界: 所属する中心までの距離（2 乗ではない距離）より大きいとわかっている値
- 下界: ほかのどの中心までの距離よりも小さいとわかっている値
中心が動いたら、三角不等式から、上界には所属する中心が動いた距離を足し、下界からはほかの中心が動いた距離の最大値を引けば、どちらも正しい界のままです。上界 < 下界 なら、所属する中心がいちばん近いままなので、その点は距離を 1 つも測らずに済みます。 「上界 < 所属する中心からいちばん近いほかの中心までの距離の半分」でも同じことが言えます。決めきれない点だけ、所属する中心までの距離を測り直し、それでもだめなら探し直します。
Iを 0から 点数 - 1まで繰り返す
所属 = ラベル一覧[I]
上界 = 上界基準[I] + 累計移動[所属]
境界 = 下界基準[I] - 最大移動累計
もし、半径一覧[所属] > 境界ならば
境界 = 半径一覧[所属]
ここまで
もし、上界 + クラスタ確定幅 × (1 + 上界) < 境界ならば
続ける
ここまで
点 = 点一覧[I]
所属二乗 = 点と(中心一覧[所属])のクラスタ二乗距離
上界 = 所属二乗の平方根
もし、上界 + クラスタ確定幅 × (1 + 上界) < 境界ならば
上界基準[I] = 上界 - 累計移動[所属]
続ける
ここまで
近傍 = 点と中心一覧を所属と所属二乗と(近隣表[所属])でクラスタ近傍探索
ラベル一覧[I] = 近傍[0]
上界基準[I] = 近傍[1] - 累計移動[近傍[0]]
下界基準[I] = 近傍[2] + 最大移動累計
ここまで
なでしこ3 版ならではの工夫が 2 つあります。
- 界を「累計との差」で持つ。上界そのものではなく「上界 − その時点の累計移動」を覚えておけば、距離を測らなかった点の値を毎回書き換えずに済みます。書き換えも 1 命令ずつ数えられるからです
- 探し直しは、所属する中心から近い順に比べて打ち切る。中心ごとにほかの中心を近い順に並べた表（
クラスタ近隣表作成。並べ替えは安定な表数値ソート）を作っておき、距離(所属, 乙) > 2 × 距離(点, 所属)になったところで、それより遠い中心は測らずに打ち切ります
打切距離 = 2 × 所属距離 + クラスタ確定幅 × (1 + 所属距離)
Jを 0から (近隣番号の要素数) - 1まで繰り返す
もし、近隣距離[J] > 打切距離ならば
見積 = 近隣距離[J] - 所属距離
もし、見積 < 下界ならば
下界 = 見積
ここまで
抜ける
ここまで
K = 近隣番号[J]
候補二乗 = 点と(中心一覧[K])のクラスタ二乗距離
候補 = 候補二乗の平方根
# 遠い（候補二乗 > 最短二乗）ことがほとんどなので、その比べ方を先にする。
もし、候補二乗 > 最短二乗ならば
もし、候補 < 下界ならば
下界 = 候補
ここまで
違えば、もし、(候補二乗 < 最短二乗 または K < 最近)ならば
もし、最短 < 下界ならば
下界 = 最短
ここまで
最近 = K
最短二乗 = 候補二乗
最短 = 候補
違えば
もし、候補 < 下界ならば
下界 = 候補
ここまで
ここまで
ここまで
同点の扱いを変えないことが肝心です。比べる順が「前から順」ではなくなったので、< だけでは「番号の小さい中心」を選べません。測った中心どうしは距離の 2 乗で比べ、同じなら番号の小さいほうを選びます。測らずに飛ばすのは、距離が 確かに大きい 中心だけです。
「確かに」のために、クラスタ確定幅（1e-9）の幅を取っています。上界と下界は平方根と足し引きを重ねた浮動小数点数なので、ほぼ同じ距離のところで丸めの誤差が大小を逆転させるかもしれません。幅の内側では「決めきれない」として距離を測り直すので、割り当てを決めるのは常に、測った距離の 2 乗どうしの比べ合いです。
素朴な繰り返しと同じ結果になることを確かめる¶
速くした割り当てが素朴な割り当てと同じであることは、テストで守ります。テストの中に、14.7 節の クラスタ割当 で毎回すべての中心と比べる 比較用素朴学習 を書き、ラベル・中心・SSE・反復回数を丸ごと比べます。
一致件数 = 0
試行件数 = 0
格子設定一覧 = [[6, 3, 1], [6, 5, 2], [3, 4, 3], [2, 6, 4]]
Iを 0から (格子設定一覧の要素数) - 1まで繰り返す
設定 = 格子設定一覧[I]
格子点 = 60と(設定[0])を(設定[1])と(設定[2])で比較用格子点作成
Kを 1から 6まで繰り返す
Sを 0から 1まで繰り返す
初期 = 格子点のKをSでクラスタ初期中心選択
速い結果 = 格子点を初期と300でクラスタ中心学習
素朴結果 = 格子点を初期と300で比較用素朴学習
試行件数 = 試行件数 + 1
もし、(速い結果をJSONエンコード) = (素朴結果をJSONエンコード)ならば
一致件数 = 一致件数 + 1
違えば
「不一致: 設定 {設定} クラスタ数 {K} シード {S}」を表示
ここまで
ここまで
ここまで
ここまで
「格子点で、距離を測る点を減らしても素朴な繰り返しとラベル・中心・SSE・反復回数が一致する」で一致件数と試行件数が一致確認
点は、自作の乱数で 0〜4 程度の整数を座標にした格子点 です。座標が整数なので、2 つの中心までの距離がちょうど同じになる点や、重なった点（同じ座標の点が初期中心に 2 回選ばれ、空のクラスタができる）が現れます。6 列の点（距離を展開した書き方）と、2・3 列の点（繰り返しの書き方）の両方を含めました。
このテストが空振りしていないことは、壊して確かめました。同点で番号の小さい中心を選ぶ条件を消すと、比べた 48 通りの一部が不一致になり、幅を負の値にして「確かに」を崩しても不一致になります。
命令の数の推移¶
報告 1 回（エルボー法の 100 回の学習と、クラスタごとの平均）の命令の数を、手を入れるたびに測りました。数値は架空の 440 件のデータでのものです。
| 手を入れたこと | 命令の数 |
|---|---|
| 最初の書き方（100 回の学習だけ） | 473,717,622 |
Hamerly の方法、中心の更新を 表列合計 に |
331,340,465 |
| 探し直しを近い順にして打ち切る、界を累計との差で持つ | 210,704,449 |
| 6 列の距離を定数の添字で展開する | 168,939,873 |
| 判定を関数にせず式で書く、クラスタ数 5 の結果を使い回す | 142,355,220 |
| 探し直しで「遠い」を先に比べる | 138,721,864 |
最後の内訳は、距離の計算 5406 万、探し直し 3725 万、繰り返しの本体 3393 万、中心の更新 815 万命令です。6 列の距離の展開は、繰り返しの書き方を 6 列のときだけ定数の添字に置き換えたもので、足す順は変えていません（0 + 最初の項 は最初の項そのものなので、距離 = 差 × 差 から始めても値は同じです）。7 列の点と、6 列に 0 を 1 列足した点で値が同じになることもテストにしました。
クラスタ数ごとに見ると、次のようになりました（架空のデータ。初期中心 10 通りの合計）。
| クラスタ数 | 更新の回数 | 命令の数 |
|---|---|---|
| 1 | 20 | 1,469,182 |
| 2 | 88 | 4,360,380 |
| 3 | 142 | 8,366,253 |
| 4 | 178 | 12,192,218 |
| 5 | 179 | 13,019,335 |
| 6 | 191 | 16,409,081 |
| 7 | 194 | 18,698,348 |
| 8 | 199 | 19,852,908 |
| 9 | 187 | 21,145,329 |
| 10 | 188 | 22,808,754 |
更新の回数は Python の参照実装で数えた値で、なでしこ3 版の最初の書き方で数えた合計 1566 回と一致しました。報告 1 回の実行時間は gonako 3.8.8 で約 3.0 秒、同じ規則の Python の参照実装（NumPy を使わない素の Python）で約 3.3 秒でした。
余裕はどれだけあるか¶
実データでも 2 億命令に収まるかは、実データで数えるまで分かりません。代わりに、架空のデータを乱数の種とばらつきの大きさを変えて 5 通り作り、報告 1 回の命令の数を測りました。
| 架空のデータ | 命令の数 |
|---|---|
| 既定（上の表のもの） | 138,721,864 |
| 種 1 | 150,681,889 |
| 種 2・ばらつき大 | 118,112,419 |
| 種 3・ばらつき特大 | 113,475,614 |
| 種 4・ばらつき小 | 129,033,363 |
1.1 億〜1.5 億命令で、上限まで 1.3〜1.8 倍の余裕です。テストファイルは実データの報告とほかの検査を 1 回の実行で行うので、既定の架空のデータを置いたときのテストファイル全体では 1 億 4927 万命令でした。
自作のデータと参照実装で突き合わせる¶
学習データはこの環境に無いので、数値の正しさは次の 2 つで確かめました。
- 自作の 14 行のデータ（値は自分で考えたもの）で報告を組み立てる関数を直接呼ぶテストを書き、期待値は Python の参照実装（
java.util.Randomと同じ乱数、Haskell 版と同じ規則の標準化・K-means・丸め）で求めた - 実データと同じ形の架空のデータ（列名・440 行・欠損なし・整数の支出額。配布データではない）を 5 通り作り、なでしこ3 版と参照実装の報告を比べた
どちらも、報告の文字列が 1 文字も違いませんでした。 架空のデータを ML_DATA_DIR で指して ./tools/check.sh を動かすと、実データのテストは期待値（ほかの言語版の値）との値の不一致だけで落ち、書き間違いのエラーでは止まらないことも確かめました。
実行して結果を表示する¶
学習データのある環境では、次のテストが期待する報告を返すことを確かめます。
「クラスタ数ごとの SSE とクラスタごとの平均支出額がほかの言語版と一致する」で(第14章実行)と実報告が一致確認
期待する出力は、Haskell 版のテストと同じ次の報告です（Haskell 版の記事によれば Java 版・Scala 版・Clojure 版・Elixir 版・PHP 版とも一致しています）。Go 版は初期中心を選ぶ乱数が math/rand なので、SSE とクラスタの内訳が違います。
データ件数: 440（支出額 6 列）
クラスタ数ごとの SSE（初期中心 10 通りの最小値）:
クラスタ数 SSE
1 2640.00
2 1954.18
3 1614.52
4 1334.36
5 1085.27
6 947.20
7 888.22
8 775.24
9 690.81
10 618.17
クラスタ数 5 のクラスタごとの件数と平均支出額:
クラスタ 件数 Fresh Milk Grocery Frozen Detergents_Paper Delicassen
2 265 8909 2967 3804 2248 989 962
1 96 5509 10556 16478 1420 7199 1659
0 65 31117 4260 5374 7225 849 2286
3 10 15965 34709 48537 3055 24875 2943
4 4 52022 31696 18491 29826 2699 19656
実データでの実測について
この章は学習データの無い環境で書いたため、上の報告はまだなでしこ3 版の実装で実測していません。実データで 2 億命令に収まるかどうかも、まだ確かめていません。学習データを配置した環境で npx gulp apps:check:nadesiko3 を実行し、保留が消えて通ることを確かめてから、この注記を外します。自作の 14 行のデータと、実データと同じ形の架空の 440 件のデータ 5 通りでは、Python の参照実装と報告が 1 文字も違わないことを確かめています。
\(k\) = 1 の SSE は、標準化した 6 列の分散がどれも 1 なので 440 × 6 = 2640 に決まります。架空のデータでも 2640.00 になり、ライブラリが無くても使える検算になっています。エルボー法で読むと、SSE の減り方ははっきりした肘を見せず、ほかの言語版と同じくクラスタ数を 5 にしました。
14.15 可視化について¶
Notebook による探索と可視化の節は設けません。クラスタの散布図・エルボー図・クラスタごとの支出額のヒートマップは、Python 版 と Kotlin 版 の「Notebook による探索と可視化」の節を参照してください。
14.16 何が置き換えられて、何が置き換えられないのか¶
何も置き換えられません。 なでしこ3 には K-means のライブラリも行列の命令も無いので、この章は自作だけで完結します。代わりに、実装の正しさを次の 4 つで支えています。
| 支え | この章での形 |
|---|---|
| 決まっている答え | \(k\) = 1 の SSE が件数 × 列数（2640）になる |
| アルゴリズムの性質 | 離れたかたまりが分かれる、局所解が起きる、クラスタ数を増やすと SSE が減る |
| 素朴な実装との一致 | 距離を測る点を減らしても、格子点の 48 通りでラベル・中心・SSE・反復回数が素朴な繰り返しと一致する |
| 参照実装とほかの言語版 | 自作・架空のデータで参照実装と報告が一致する。実データはほかの言語版の値を期待値にする |
3 つ目はなでしこ3 版で加わったものです。命令の上限のために速い書き方を選んだぶん、「速い書き方は素朴な書き方と同じ答えを出す」ことを、テストで支え続ける必要があります。
14.17 品質チェック¶
cd apps/nadesiko3
./tools/check.sh
文法・整形・テスト・テストが最後まで走ったか・関数の網羅の 5 つを確かめます。この章で関数が 19 個増え、関数の網羅は 67 件中 67 件（直接 57 件、間接 10 件）になりました。第 14 章のテストは 50 件（学習データが無いときは実データの 1 件が保留）です。
14.18 まとめ¶
この章では、K-means をなでしこ3 で自作しました。なでしこ3 に固有の論点は次のとおりです。
- gonako は 1 回の実行で 2 億命令を超えると止まる — 処理系の既定の設定で、コマンドから変える引数は無い。素朴な K-means は報告に 4.7 億命令を要した
- 命令を数えてから減らす — 処理系に命令の数と関数ごとの内訳を表示させて測った。繰り返しの制御は足し算と同じくらい高く、Go で書かれた命令は中で何個足しても数命令
- 結果を変えずに距離を測る点を減らす — Hamerly の方法で、距離の上界と下界から「所属が変わらない」とわかる点を飛ばした。界は累計との差で持ち、探し直しは近い順に比べて打ち切る。1.4 億命令になった
- 同点の扱いを守る — 比べる順が変わっても、距離の 2 乗が同じなら番号の小さい中心を選ぶ。飛ばすのは幅（1e-9）を取って「確かに遠い」中心だけ
- 速い書き方は素朴な書き方でテストする — 同点や空のクラスタが起きやすい格子点で、ラベル・中心・SSE・反復回数を丸ごと比べる
表列合計は 0 から前の行へ順に足す — ソースで確かめてから使い、中心の値を最後の桁まで保った二乗は組み込みの命令 — 変数名にすると「関数『二乗』に代入できません」になる
TODO リスト（この章の完了時点）:
- 支出額の列だけを読み込む
- 列ごとに標準化する
- 2 点間の距離の 2 乗を求める
- 各点を最も近い中心に割り当てる
- クラスタごとの平均を新しい中心にする
- SSE を計算する
- 中心が変わらなくなるまで繰り返す
- 初期中心をシードで選び、何通りか試して SSE が最小の結果を採る
- クラスタ数ごとの SSE を並べる（エルボー法）
- クラスタごとの件数と平均支出額をまとめる（実データでの実測は未確認。14.14 節の注記を参照）
なでしこ3 版のほかの章は なでしこ3 版のトップ から辿れます。