kmeans++初期化のアルゴリズム(sklearn)を理解する

## k-means++初期化の概要
### k-means++でやりたいこと
k-meansではデータ構造自体からデータをクラスタリングする。inputとして教師データは不要で、クラスタリングするクラス数と初期クラス中心点を指定する。　　
この初期クラス位置が重要で、k-meansアルゴリズムの性質上、そのとり方によって結果が変わってしまう。k-means++では初期値を良い感じに配置するためのアルゴリズムを与える。
一言でいうと、中心点が、データ点の集合のなかでいい感じに分散するようなアルゴリズムである。

k-means++の初期化についての概念と簡単な実装は[Pythonでの数値計算ライブラリNumPy徹底入門 \- DeepAge](https://deepage.net/features/numpy/)がとてもわかりやすい。

ここでは実際の実装(sklearn)におけるk-means++の初期化アルゴリズムを調べる。

### 初期化では何をしているのか?
やりたいことは、外れ値に強く、かつ各中心点を適度にばらつかせること。
一番最初の中心点はデータ点からランダムに一点を選ぶ。空間からランダムに選ぶのではなく、データ点から選ぶ。
2点目以降の中心点の選び方について図で示す。

![Screenshot from 2020-01-29 10-55-52.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/188904/99e39246-3bcb-e51e-3f3d-372af8147c9c.png)

* 各点から一番近い中心点までの距離を計算する。その距離の2乗に応じて確率的に次の中心点候補を複数選ぶ。(このときすでに中心点に選ばれている点は距離が0なので選ばれる確率も0となる)。確率的に選ぶのは外れ値対策である。
* 次に複数候補点から一つを選ぶ。その基準として、各点における最近傍の中心点までの距離の和を計算する。この和が小さいほど、中心点が各データに近くなるよう配置されていることになるので、和が小さい候補点を次の中心点として採用する。今度は決定論的に最小値をとるものを採用する。

このように２つのステップを踏んで2点目以降を順次きめていく。

距離の総和を最初から使うと、計算量が増えることと、外れ値が選ばれやすくなるのがデメリットと思われる。

それでは具体的な実装をsklearnのコードで確認していく。
[sklearnのk-means++の初期化のソースコード](https://github.com/scikit-learn/scikit-learn/blob/1495f69242646d239d89a5713982946b8ffcf9d9/sklearn/cluster/k_means_.py)

## sklearnでのk-menans++初期化処理の流れ

### 初期化処理の入力
* データX: array, shapeは(サンプルサイズ, 特徴次元) 
* クラスター数: int
* X2乗ノルム: |X|^2, array, shapeは(サンプルサイズ, )  
中心点との2乗距離を計算するときに、  
|x-y|^2 = dot(x, x) - 2 * dot(x, y) + dot(y, y)  
と計算するほうが効率がよいので、先に計算して引数で渡す
* random_state: 最初の点をランダムに生成するのに使用する
* n_local_trials(optional): 各中心点を決めるときの候補の数. defaultは(2+log(k))

### 最初の中心点をデータ点の中からランダムに選ぶ  
特徴空間全体からランダムに選ぶのではなく、データ点の一つを選ぶ

```python
center_id = random_state.randint(n_samples)
centers[0] = X[center_id]
```
### 各点から中心点との距離と、その距離の総和を計算する
どちらも次の候補点選びに使われる

```python
closest_dist_sq = euclidean_distances(
    centers[0, np.newaxis], X, Y_norm_squared=x_squared_norms,
    squared=True)
current_pot = closest_dist_sq.sum()  # potはpotentialのpot
```

### 次の候補点を複数(n_local_trials)選ぶ
やりたいことは最近房の中心点までの距離の2乗に応じて確率的に選択するということ。つまり距離0、距離2、距離4の3点があるとすると、それぞれ、0/(0+2^2+4^4)=0, 2^2/(0+2^2+4^4)=4/20, 4^4/(0+2^2+4^4)=16/20, の確率で選びたい。これを累積和のcumsum関数と乱数を使って実現している. 総和×(0-1乱数)がcumsumのどこにぶつかるか、という感じで、距離が大きいほど選ばれやすくなる.
![ファイル_001.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/188904/d5a440b0-a4a7-bdae-c312-231b79add61d.png)

```python
rand_vals = random_state.random_sample(n_local_trials) * current_pot
candidate_ids = np.searchsorted(stable_cumsum(closest_dist_sq), rand_vals)
```
### 候補点から一つを次の中心点として選ぶ
各点における最近傍の中心点までの距離の和が小さいほど、中心点が各データに近くなるよう配置されていることになるので、和が小さい候補点を次の中心点として採用する。
採用したら前の節の処理に戻り、全ての点を選ぶまで複数選択→決定の処理を繰り返す。
![ファイル_001 (1).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/188904/1a4338a9-c431-29d9-7500-7f511e54a0c5.png)

```python
for trial in range(n_local_trials):
    # Compute potential when including center candidate
    new_dist_sq = np.minimum(closest_dist_sq, distance_to_candidates[trial])
    new_pot = new_dist_sq.sum()

    # Store result if it is the best local trial so far
    if (best_candidate is None) or (new_pot < best_pot):
        best_candidate = candidate_ids[trial]
        best_pot = new_pot
        best_dist_sq = new_dist_sq
```

## ソースコード全体
```python
# scikit-learn/sklearn/cluster/k_means_.py 2019-11-3
def _k_init(X, n_clusters, x_squared_norms, random_state, n_local_trials=None):
    """Init n_clusters seeds according to k-means++

    Parameters
    ----------
    X : array or sparse matrix, shape (n_samples, n_features)
        The data to pick seeds for. To avoid memory copy, the input data
        should be double precision (dtype=np.float64).

    n_clusters : integer
        The number of seeds to choose

    x_squared_norms : array, shape (n_samples,)
        Squared Euclidean norm of each data point.

    random_state : int, RandomState instance
        The generator used to initialize the centers. Use an int to make the
        randomness deterministic.
        See :term:`Glossary <random_state>`.

    n_local_trials : integer, optional
        The number of seeding trials for each center (except the first),
        of which the one reducing inertia the most is greedily chosen.
        Set to None to make the number of trials depend logarithmically
        on the number of seeds (2+log(k)); this is the default.

    Notes
    -----
    Selects initial cluster centers for k-mean clustering in a smart way
    to speed up convergence. see: Arthur, D. and Vassilvitskii, S.
    "k-means++: the advantages of careful seeding". ACM-SIAM symposium
    on Discrete algorithms. 2007

    Version ported from http://www.stanford.edu/~darthur/kMeansppTest.zip,
    which is the implementation used in the aforementioned paper.
    """
    n_samples, n_features = X.shape

    centers = np.empty((n_clusters, n_features), dtype=X.dtype)

    assert x_squared_norms is not None, 'x_squared_norms None in _k_init'

    # Set the number of local seeding trials if none is given
    if n_local_trials is None:
        # This is what Arthur/Vassilvitskii tried, but did not report
        # specific results for other than mentioning in the conclusion
        # that it helped.
        n_local_trials = 2 + int(np.log(n_clusters))

    # Pick first center randomly
    center_id = random_state.randint(n_samples)
    if sp.issparse(X):
        centers[0] = X[center_id].toarray()
    else:
        centers[0] = X[center_id]

    # Initialize list of closest distances and calculate current potential
    closest_dist_sq = euclidean_distances(
        centers[0, np.newaxis], X, Y_norm_squared=x_squared_norms,
        squared=True)
    current_pot = closest_dist_sq.sum()

    # Pick the remaining n_clusters-1 points
    for c in range(1, n_clusters):
        # Choose center candidates by sampling with probability proportional
        # to the squared distance to the closest existing center
        rand_vals = random_state.random_sample(n_local_trials) * current_pot
        candidate_ids = np.searchsorted(stable_cumsum(closest_dist_sq),
                                        rand_vals)
        # XXX: numerical imprecision can result in a candidate_id out of range
        np.clip(candidate_ids, None, closest_dist_sq.size - 1,
                out=candidate_ids)

        # Compute distances to center candidates
        distance_to_candidates = euclidean_distances(
            X[candidate_ids], X, Y_norm_squared=x_squared_norms, squared=True)

        # Decide which candidate is the best
        best_candidate = None
        best_pot = None
        best_dist_sq = None
        for trial in range(n_local_trials):
            # Compute potential when including center candidate
            new_dist_sq = np.minimum(closest_dist_sq,
                                     distance_to_candidates[trial])
            new_pot = new_dist_sq.sum()

            # Store result if it is the best local trial so far
            if (best_candidate is None) or (new_pot < best_pot):
                best_candidate = candidate_ids[trial]
                best_pot = new_pot
                best_dist_sq = new_dist_sq

        # Permanently add best center candidate found in local tries
        if sp.issparse(X):
            centers[c] = X[best_candidate].toarray()
        else:
            centers[c] = X[best_candidate]
        current_pot = best_pot
        closest_dist_sq = best_dist_sq

    return centers
```