SVDでお手軽irisデータセットの可視化と分類

# SVD、もしかして需要がない?
昨日の朝に、[特異値分解(SVD)](https://qiita.com/gyu-don/items/2e4fd7e945d8c349afb5)に関する記事を書きました。ですが、あまりにもLGTMが少なくて、「もしかして特異値分解って需要がない?」という気持ちになりました。

今回は、SVDのよさを知ってほしくて、記事を書きます。
SVDは次元が大きいデータを整理して次元を圧縮できる、とても素晴らしいアルゴリズムなのです。

例として、機械学習の練習でおなじみの、アヤメの種類を特徴から分類するIrisデータセットを使って次元削減して、2次元で図示してみましょう。
**この記事は、機械学習よりは、次元削減の方法や、それによって図示がしやすくなることに重みを置いています。**

## Irisデータセットをscikit-learnから取ってくる

このデータセットでは、4つの次元の特徴量と、0, 1, 2のラベルが用意されています。4次元というのは人間には難しく、なかなか図示もしづらいですね。


```python
from itertools import combinations
import numpy as np
import pandas as pd
from sklearn.datasets import load_iris
from sklearn.model_selection import train_test_split
from sklearn.svm import LinearSVC
import matplotlib.pyplot as plt

iris = load_iris(as_frame=True)
iris.frame
```

<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>sepal length (cm)</th>
      <th>sepal width (cm)</th>
      <th>petal length (cm)</th>
      <th>petal width (cm)</th>
      <th>target</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>5.1</td>
      <td>3.5</td>
      <td>1.4</td>
      <td>0.2</td>
      <td>0</td>
    </tr>
    <tr>
      <th>1</th>
      <td>4.9</td>
      <td>3.0</td>
      <td>1.4</td>
      <td>0.2</td>
      <td>0</td>
    </tr>
    <tr>
      <th>2</th>
      <td>4.7</td>
      <td>3.2</td>
      <td>1.3</td>
      <td>0.2</td>
      <td>0</td>
    </tr>
    <tr>
      <th>3</th>
      <td>4.6</td>
      <td>3.1</td>
      <td>1.5</td>
      <td>0.2</td>
      <td>0</td>
    </tr>
    <tr>
      <th>4</th>
      <td>5.0</td>
      <td>3.6</td>
      <td>1.4</td>
      <td>0.2</td>
      <td>0</td>
    </tr>
    <tr>
      <th>...</th>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
    </tr>
    <tr>
      <th>145</th>
      <td>6.7</td>
      <td>3.0</td>
      <td>5.2</td>
      <td>2.3</td>
      <td>2</td>
    </tr>
    <tr>
      <th>146</th>
      <td>6.3</td>
      <td>2.5</td>
      <td>5.0</td>
      <td>1.9</td>
      <td>2</td>
    </tr>
    <tr>
      <th>147</th>
      <td>6.5</td>
      <td>3.0</td>
      <td>5.2</td>
      <td>2.0</td>
      <td>2</td>
    </tr>
    <tr>
      <th>148</th>
      <td>6.2</td>
      <td>3.4</td>
      <td>5.4</td>
      <td>2.3</td>
      <td>2</td>
    </tr>
    <tr>
      <th>149</th>
      <td>5.9</td>
      <td>3.0</td>
      <td>5.1</td>
      <td>1.8</td>
      <td>2</td>
    </tr>
  </tbody>
</table>
</div>


試しに、各2次元について、グラフを描いてみましょう。


```python
axes = list(combinations(range(4), 2))
plt_cols = 3
plt_rows = (len(axes) - 1) // plt_cols + 1
targets = [0, 1, 2]

plt.figure(figsize=(12, 8))
for i, (x, y) in enumerate(axes):
    plt.subplot(plt_rows, plt_cols, i + 1)
    plt.xlabel(iris.frame.columns[x])
    plt.ylabel(iris.frame.columns[y])
    for t in targets:
        plt.plot(iris.frame[iris.frame.target == t].iloc[:, x], iris.frame[iris.frame.target == t].iloc[:, y], ".")
```

![output_4_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/88256/6faec302-cd60-3bbc-5ad1-32071712648c.png)
    

各軸ごとにプロットしてみました。おぼろげに何かは見えるけど、4次元を直感的に捉えるのはあまりにも難しすぎます。

# SVDをやってみる
SVDをやるだけなら、train, testという概念は特になく、単なる次元削減という側面が強いです。情報を落とすことになるので、単に機械学習をしたいのであれば、SVDをしなくてもいいかもしれません。けれど、次元圧縮と機械学習の両方を同時にやってみましょう。

そのために、データセットをtrainデータとtestデータに分けます。そして、trainデータに対してSVDをやってみます。

```python
x_train, x_test, y_train, y_test = train_test_split(iris.data, iris.target, test_size=0.2, random_state=1)
```

trainデータはデータ数が114、特徴量の数が4です。


```python
print(x_train.shape) # => (120, 4)
```

SVD自体は、教師あり機械学習ではないので、ラベルのデータは使いません。単に `x_train` を次元削減します。


```python
u, s, vt = np.linalg.svd(x_train, full_matrices=False)
```

112x4の行列である `x_train` を、SVDにより3つの行列の積に分割しました(`s`は対角行列なので対角成分だけが返ってきています)

`u * s @ vt` 、あるいは同じ意味ですが `np.dot(u * s, vt)` が、`x_train`と一致します。


```python
print(np.allclose(x_train, u * s @ vt)) # => True
```

あるいは、式変形により、`x_train`, `vt`, `s`から`u`を導く式も作れます。(今回はうまく合いましたが、割り算を使っているので、`s`の値によっては、浮動小数点誤差が大きく合わないかもしれません。しかし、この後、次元削減により小さな値は無視することにします)
この式は、trainデータでSVDした結果をtestデータに使うときに必要になります。


```python
print(np.allclose(x_train @ vt.T / s, u)) # => True
```

ところで、`s`の値を見てみましょう。


```python
print(s) # => [85.99229833 15.66745575  3.29372488  1.58329849]
```

`s`は、数字が大きい順になっています。`u`と`vt`は、`x_train`を行列として見たとき、行列を長さを変えずに変換するもので、`s`の部分が長さの変換に相当します。そして、`s`が小さい部分というのは、あまり変換に寄与しないということです。

**この値を見るに、最初の2つだけで十分だと思いませんか?**

ということで、次元の削減ができます！
`s`のデータを削って、それと次元が合うように `u`, `vt`も削ります。


```python
x_train_2d = u[:, :2] * s[:2] @ vt[:2, :2]
for t in targets:
    plt.plot(x_train_2d[y_train == t][:, 0], x_train_2d[y_train == t][:, 1], ".")
```

![output_19_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/88256/1dc347be-80ad-1db2-8b3d-b3af619bebaa.png)

    


## 次元削減したデータをざっくりと線形分離してしまう
情報を削減しているので、単に機械学習をするなら、元のデータの方が精度はいいかもしれません。けれど、2次元平面で図示できるというのは圧倒的に分かりやすいです。少々精度は落ちても、SVM(サポートベクターマシン)で線形分離してみましょう。

ここでは、setosa(青)とversicolor(橙)、versicolor(橙)とvirginica(緑)を分類する2本の直線を引きます。


```python
clf1 = LinearSVC()
clf1.fit(x_train_2d[y_train.isin([0, 1])], y_train[y_train.isin([0, 1])])
clf2 = LinearSVC(max_iter=2000)
clf2.fit(x_train_2d[y_train.isin([1, 2])], y_train[y_train.isin([1, 2])])
```

はい、できました。図示しましょう。

```python
# train-data
for t in targets:
    plt.plot(x_train_2d[y_train == t][:, 0], x_train_2d[y_train == t][:, 1], ".")
# line
x0, x1 = plt.xlim()
y0, y1 = plt.ylim()
plt.plot([x0, -(clf1.coef_[0, 1] * y1 + clf1.intercept_[0]) / clf1.coef_[0, 0]], [-(clf1.coef_[0, 0] * x0 + clf1.intercept_[0]) / clf1.coef_[0, 1], y1], color="black")
plt.plot([-(clf2.coef_[0, 1] * y0 + clf2.intercept_[0]) / clf2.coef_[0, 0], x1], [y0, -(clf2.coef_[0, 0] * x1 + clf2.intercept_[0]) / clf2.coef_[0, 1]], color="black")
```

![output_22_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/88256/c4d68592-9f61-0b96-dcd1-4d9d39a6d02c.png)
    
所詮は2次元の線形SVMですが、多次元で計算するよりも圧倒的に図として見やすいですし、「何かよさそうな特徴量を2個選びました」よりは、「SVDで寄与の大きい次元を選びました」の方が、選んだ根拠としても強いです。さて、ここにtestデータを載せてみたら、どうなるでしょう?

## testデータを載せてみる

testデータ用の`u`をテストデータと`vt`, `s`から作り、さらに、testデータを次元圧縮するのですが、計算してみると、一度の計算にまとめられることが分かります。

```python
# testデータの計算
# u_test_2d = x_test.to_numpy() @ vt.T[:, :2] / s[:2]
# x_test_2d = u_test_2d * s[:2] @ vt[:2, :2]
# と、2段階でできるが、上の2つをまとめると以下のようになる
x_test_2d = x_test.to_numpy() @ vt.T[:, :2] @ vt[:2, :2]
```

さて、これをプロットしましょう。

```python
# train-data
for t in targets:
    plt.plot(x_train_2d[y_train == t][:, 0], x_train_2d[y_train == t][:, 1], ".", alpha=0.2)
# test-data
colors = plt.rcParams['axes.prop_cycle'].by_key()['color'] # point colors
for t in targets:
    plt.plot(x_test_2d[y_test == t][:, 0], x_test_2d[y_test == t][:, 1], "*", color=colors[t])

# line
x0, x1 = plt.xlim()
y0, y1 = plt.ylim()
plt.plot([x0, -(clf1.coef_[0, 1] * y1 + clf1.intercept_[0]) / clf1.coef_[0, 0]], [-(clf1.coef_[0, 0] * x0 + clf1.intercept_[0]) / clf1.coef_[0, 1], y1], color="black")
plt.plot([-(clf2.coef_[0, 1] * y0 + clf2.intercept_[0]) / clf2.coef_[0, 0], x1], [y0, -(clf2.coef_[0, 0] * x1 + clf2.intercept_[0]) / clf2.coef_[0, 1]], color="black")
```

![output_25_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/88256/050fa6d8-602b-0346-0e87-082a453f50e3.png)

trainデータの色を薄くして、testデータは星でプロットしてみました。
どうですか? なかなか綺麗じゃないですか？