「次元の呪い」って何？Pythonで次元削減をやさしく学ぶ（PCA・LDA・Kernel PCA編）

「次元の呪い」──これは、機械学習においてデータの次元（特徴量の数）が多すぎることで、モデルの性能が低下したり、計算負荷が増大したりする現象を指す。こうした問題を軽減するために有効なアプローチの1つが次元削減であり、大きく「特徴量選択」と「特徴量抽出」の2種類に分けられる。

[前回の記事](https://qiita.com/Djay/items/0d2eee691533facd8f96)では、不要な特徴量を取り除く「特徴量選択」について解説した。今回は、次元の呪いの本質と、特徴量抽出による次元削減の代表的な手法（PCA・LDA・Kernel PCA）について、Pythonコードとともにわかりやすく紹介する。

本記事は、[Python機械学習プログラミング 達人データサイエンティストによる理論と実践](https://amzn.to/4ktOMjQ)の第5章を参考にしている。

本記事では以下の流れで説明する。
- [1 特徴量抽出](#1-特徴量抽出)
- [2 主成分分析（PCA）](#2-主成分分析pca)
    - [2.1 PCAの実装](#21-pcaの実装)
    - [2.2 主成分数の決定方法](#22-主成分数の決定方法)
- [3 線形判別分析（LDA）](#3-線形判別分析lda)
    - [3.1 LDAの実装](#31-ldaの実装)
- [4 カーネル主成分分析（Kernel PCA）](#4-カーネル主成分分析kernel-pca)
    - [4.1 Kernel PCAの実装](#41-kernel-pcaの実装)
- [5 まとめ](#5-まとめ)


# 1 特徴量抽出
特徴量選択アルゴリズムでは、元の特徴量は変換されることなくそのままの形で維持されていた。一方、特徴量抽出では、データが新しい特徴量空間に変換または射影される。特に、特徴量抽出は正則化されていないモデルを扱う時、**次元の呪い**（特徴量の次元（変数の数）が増えることで、データ解析や学習が困難になること）を減らすことで予測性能を向上させる目的でも利用できる。すなわち、モデルが高次元データに対して過学習しやすい場合、特徴量抽出を通じてよりコンパクトな表現にすることで、予測精度の向上が期待できる。

# 2 主成分分析（PCA）
**主成分分析（PCA）**は様々な分野で広く使われている教師なし線形変換法。PCAの目的は、高次元データにおいて分散が最大となる方向を見つけ出し、元の次元と同じかそれよりも低い次元の新しい部分空間へ射影すること。
変換後の第1主成分の分散が最大となり、第2主成分はそれに直交するような形になる。PCAはデータのスケーリングに対して敏感なので、PCAの前に特徴量を標準化する必要がある。

主成分分析のプロセスは以下のようになる。
1.　d次元のデータセットを標準化する
2.　標準化したデータセットの共分散行列（covariance matrix）を作成する
2つの特徴量$x_j$と $x_k$の間の共分散は

$$
\sigma_{jk} = \frac{1}{n} \sum_{i=1}^{n} (x_{j}^{(i)} - \mu_j)(x_{k}^{(i)} - \mu_k)
$$

ゆえに、3つの特徴量からなる共分散行列は以下のようにかける。

$$
\Sigma = 
\begin{bmatrix}
\sigma_{1}^2 & \sigma_{12} & \sigma_{13} \\\\
\sigma_{21} & \sigma_{2}^2 & \sigma_{23} \\\\
\sigma_{31} & \sigma_{32} & \sigma_{3}^2
\end{bmatrix}
$$

3.　共分散行列を固有ベクトルと固有値に分解する
固有ベクトル$\nu$はスカラー$\lambda$（固有値）を用いて以下のように記述される。
$$
\Sigma \nu = \lambda\nu
$$

4.　固有値を降順でソートすることで、対応する固有ベクトルをランク付けする
5.　削減後の次元（k）を設定して、最も大きいk個の固有値に対応するk個の固有ベクトルを選択する(k<d)
6.　上位k個の固有ベクトルから射影行列$W$を作成する
7.　射影行列$W$を使って、d次元の入力データセット$X$を変換し、新しいk次元の特徴量部分空間を取得する
$$
\mathbf{X'} = \mathbf{X}\mathbf{W}
$$


## 2.1 PCAの実装
ここでは、[UCI Machine Learning Repository](https://archive.ics.uci.edu/dataset/109/wine)のWineデータセットを使う。

<details><summary>コードを表示</summary>

```python
# UCI Machine Learning RepositoryからWineデータセットを取得
import pandas as pd
import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler

df_wine = pd.read_csv('https://archive.ics.uci.edu/'
                      'ml/machine-learning-databases/wine/wine.data',
                      header=None)

df_wine.columns = ['Class label', 'Alcohol', 'Malic acid', 'Ash',
                   'Alcalinity of ash', 'Magnesium', 'Total phenols',
                   'Flavanoids', 'Nonflavanoid phenols', 'Proanthocyanins',
                   'Color intensity', 'Hue', 'OD280/OD315 of diluted wines',
                   'Proline']

print('Class labels', np.unique(df_wine['Class label']))
# Class labels [1 2 3]
先頭の5行を表示　
df_wine.head().iloc[:, :5]

# データセットを訓練データとテストデータに分割
X, y = df_wine.iloc[:, 1:].values, df_wine.iloc[:, 0].values
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, stratify=y, random_state=0)

# 特徴量の標準化
sc = StandardScaler()
X_train_std = sc.fit_transform(X_train)
X_test_std = sc.transform(X_test)

```
</details>

<details><summary>データを表示</summary>

| Class label | Alcohol | Malic acid | Ash | Alcalinity of ash |
|-------------|---------|------------|-----|--------------------|
| 1           | 14.23   | 1.71       | 2.43| 15.6               |
| 1           | 13.20   | 1.78       | 2.14| 11.2               |
| 1           | 13.16   | 2.36       | 2.67| 18.6               |
| 1           | 14.37   | 1.95       | 2.50| 16.8               |
| 1           | 13.24   | 2.59       | 2.87| 21.0               |

</details>

Wineデータに対して、`scikit-learn`の`PCA`クラスを用いる。訓練データを使って、PCAのモデルを適合させた上で、訓練データとテストデータを同じモデルパラメータに基づいて変換する。PCAを使って、 13次元から2次元に特徴量を削減し、返還後のデータを使って、ロジスティック回帰で分類する。

```python
from sklearn.linear_model import LogisticRegression
from sklearn.decomposition import PCA
import matplotlib.pyplot as plt

# PCAのインスタンスを生成
pca = PCA(n_components=2)
# 訓練データをPCAで変換
X_train_pca = pca.fit_transform(X_train_std)
X_test_pca = pca.transform(X_test_std)

# ロジスティック回帰のインスタンスを生成
lr = LogisticRegression(multi_class='ovr', random_state=1, solver='lbfgs')

# 変換した訓練データを使って学習
lr.fit(X_train_pca, y_train)
```

<details><summary>決定領域図を出力するコード</summary>

```python
# 決定領域を描画
x_min, x_max = X_train_pca[:, 0].min() - 1, X_train_pca[:, 0].max() + 1
y_min, y_max = X_train_pca[:, 1].min() - 1, X_train_pca[:, 1].max() + 1
xx, yy = np.meshgrid(np.arange(x_min, x_max, 0.02),
                     np.arange(y_min, y_max, 0.02))

# グリッド全点を予測
Z = lr.predict(np.array([xx.ravel(), yy.ravel()]).T)
Z = Z.reshape(xx.shape)

# 背景の決定領域を描画
plt.contourf(xx, yy, Z, alpha=0.3, cmap=plt.cm.RdYlBu)

# 訓練データを散布図で描画
plt.scatter(X_train_pca[y_train == 1, 0], X_train_pca[y_train == 1, 1],
            c='r', marker='o', label='Class 1')
plt.scatter(X_train_pca[y_train == 2, 0], X_train_pca[y_train == 2, 1],
            c='g', marker='s', label='Class 2')
plt.scatter(X_train_pca[y_train == 3, 0], X_train_pca[y_train == 3, 1],
            c='b', marker='x', label='Class 3')

plt.xlabel('Principal Component 1')
plt.ylabel('Principal Component 2')
plt.title('PCA of Wine Dataset with Decision Regions')
plt.legend(loc='upper left')
plt.show()

```
</details>

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4072443/32bb8569-c1bc-4b7c-a4bc-fde6df91e3c8.png" width="70%">

各主成分の分散説明率にアクセスしたい場合は、`n_components`パラメータをNoneに設定する。

```python
pca = PCA(n_components=None)
X_train_pca = pca.fit_transform(X_train_std)
# 分散説明率を取得
pca.explained_variance_ratio_
# array([0.36951469, 0.18434927, 0.11815159, 0.07334252, 0.06422108,
#        0.05051724, 0.03954654, 0.02643918, 0.02389319, 0.01629614,
#        0.01380021, 0.01172226, 0.00820609])
```

## 2.2 主成分数の決定方法
PCAを実行した後は、何個の主成分を用いて以降の分析を行うかを決める必要がある。その方法として、本記事では以下の3つの手法を紹介する。
1. 累積寄与率：一般的には累積寄与率が80~95%に達するまでの主成分を用いる
1. スクリープロット：主成分ごとの固有値または寄与率をプロットして、「肘」にあたる部分を探す
1. タスク依存のパフォーマンス：実務においてはおそらく最も有効な方法で、主成分の数を変えたいくつかのモデルを作成し、予測精度やクラスタリングのパフォーマンス評価によって決める

1つ目の方法は以下のようなコードから、累積寄与率を計算し、累積寄与率が80~95%に達するまでの主成分を用いる方法である。これは、「全体の分散の80〜95%を説明できる」ことを意味し、元のデータの大部分の情報を保持できる。今回の場合であれが、第4~7主成分までの用いるのが良い。

```python
# 寄与率の確認
pca_full = PCA()
pca_full.fit(X_train_std)
print('Explained variance ratio:', pca_full.explained_variance_ratio_)

# Explained variance ratio: [0.36951469 0.18434927 0.11815159 0.07334252 0.06422108 0.05051724
#  0.03954654 0.02643918 0.02389319 0.01629614 0.01380021 0.01172226
#  0.00820609]

plt.plot(np.cumsum(pca_full.explained_variance_ratio_))
plt.xlabel('Number of components')
plt.ylabel('Cumulative explained variance')
plt.title('Explained Variance by PCA Components')
plt.grid(True)
plt.show()
```
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4072443/f1368ab9-1811-4a27-b796-82931e6130a8.png" width="60%">

次に、スクリープロット法に用いるグラフを以下に記述する。この手法では、グラフが急激に下がった後に緩やかになる点（Elbow）で、重要な情報が含まれている次元数を見極める。今回の例では、第4~8主成分までの用いるのが良い。

```python
plt.plot(range(1, len(pca.explained_variance_ratio_)+1), pca.explained_variance_ratio_, marker='o')
plt.xlabel('Component number')
plt.ylabel('Explained variance ratio')
plt.title('Scree Plot')
plt.show()
```

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4072443/857534c6-15ef-4630-8e51-93e315384618.png" width="60%">

最後に、全ての特徴量を用いた場合と、第4主成分までを用いてモデルを訓練した場合の予測精度の違いを比較する。
```python
from sklearn.metrics import accuracy_score

pca = PCA(n_components=4)
X_train_pca = pca.fit_transform(X_train_std)
X_test_pca = pca.transform(X_test_std)

# PCA前のロジスティック回帰
lr_orig = LogisticRegression(multi_class='ovr', random_state=1, solver='lbfgs')
lr_orig.fit(X_train_std, y_train)
y_pred_orig = lr_orig.predict(X_test_std)

# PCA後のロジスティック回帰の予測
lr = LogisticRegression(multi_class='ovr', random_state=1, solver='lbfgs')
lr.fit(X_train_pca, y_train)
y_pred_pca = lr.predict(X_test_pca)

print('Original accuracy:', accuracy_score(y_test, y_pred_orig))
# Original accuracy: 1.0
print('PCA reduced accuracy:', accuracy_score(y_test, y_pred_pca))
# PCA reduced accuracy: 0.9259259259259259
```

以上から、第4主成分まで含めた分析結果は、13個の特徴量をもつ分析結果とそれほど変わらない予測精度をもつ（第7主成分まで含めると、予測結果は同じになる）。

# 3. 線形判別分析（LDA）
**線形判別分析**（Linear Discriminant Analysis:LDA）は特徴量抽出手法の1つ。LDAは正則化されていないモデルで次元の呪いによる過学習を抑制するのに利用できる。LDAの基本的な考え方はPCAと似ているが、PCAはデータセットにおいて分散が最大となる直交成分軸を見つけ出そうとするのに対して、LDAはクラスの分離を最適化する特徴量部分空間を見つけ出そうとする。
PCAは教師なし、LDAは教師あり線形変換であるため、LDAの方が分類タスクの特徴量抽出としては優れている。
一応、LDAはデータが以下の条件を満たしていることを前提条件としている。しかし、次元削減の観点においては条件を満たしていなくても、それなりにうまくいく。

- データが正規分布に従っている
- クラスの共分散が全く同じである
- 訓練データが互いに独立している

LDAの手順を以下にまとめる。手順4~7はPCAの手順とほぼ同じ。

1.　d次元の（d個の特徴量を持つ）データを標準化する
2.　クラスごとにd次元の平均ベクトルを計算する
各平均ベクトル$m_i$は、クラスiのデータ点に関する平均特徴量の値$\mu_i$を格納する。

$$
m_i = 
\begin{bmatrix}
\mu_{i, j1} \\\\
\mu_{i, j2} \\\\
\vdots \\\\
\mu_{i, jn} \\\\
\end{bmatrix}^\top, \quad i \in \{1, 2, \dots, c\}
$$

3.　平均ベクトルを使って、クラス間変動行列$S_B$とクラス内変動行列$S_W$を生成する
$$
S_W = \sum_{i=1}^{c}S_i
$$
$S_i$は変動行列を表し、
$$
S_i = \sum_{\mathbf{x} \in D_i} (\mathbf{x} - m_i)^T (\mathbf{x} - m_i)
$$

クラス間変動行列$S_B$は、全てのクラス$c$のデータ点を対象として計算される全体平均$m$を用いて、
$$
S_B = \sum_{i=1}^c n_i(m_i - m)^T (m_i - m)
$$

4.　行列$S_W^{-1}S_B$の固有ベクトルと対応する固有値を計算
5.　固有値を降順でソートすることで、対応する固有ベクトルをランク付けする
6.　最も大きいk個の固有値に対する固有ベクトルを選択し、変換行列$W$を生成する
7.　変換行列$W$を使ってデータ点を新しい特徴量部分空間へ射影する

## 3.1 LDAの実装
同じようにWineデータに対して、`scikit-learn`の`LinearDiscriminantAnalysis`を使って次元削減を行う。その後、削減した次元を元にロジスティック回帰を行う。

```python
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA
# LDAのインスタンスを生成
lda = LDA(n_components=2)
X_train_lad = lda.fit_transform(X_train_std, y_train)
# ロジスティック回帰のインスタンスを生成
lr = LogisticRegression(multi_class='ovr', random_state=1, solver='lbfgs')
lr = lr.fit(X_train_lad, y_train)
```

<details><summary>決定領域図を出力するコード</summary>

```python
# ロジスティック回帰のインスタンスを生成
lr = LogisticRegression(multi_class='ovr', random_state=1, solver='lbfgs')
lr = lr.fit(X_train_lad, y_train)
# 決定領域を描画
x_min, x_max = X_train_lad[:, 0].min() - 1, X_train_lad[:, 0].max() + 1
y_min, y_max = X_train_lad[:, 1].min() - 1, X_train_lad[:, 1].max() + 1
xx, yy = np.meshgrid(np.arange(x_min, x_max, 0.02),
                     np.arange(y_min, y_max, 0.02))
# グリッド全点を予測
Z = lr.predict(np.array([xx.ravel(), yy.ravel()]).T)
Z = Z.reshape(xx.shape)
# 背景の決定領域を描画
plt.contourf(xx, yy, Z, alpha=0.3, cmap=plt.cm.RdYlBu)
# 訓練データを散布図で描画
plt.scatter(X_train_lad[y_train == 1, 0], X_train_lad[y_train == 1, 1],
            c='r', marker='o', label='Class 1')
plt.scatter(X_train_lad[y_train == 2, 0], X_train_lad[y_train == 2, 1],
            c='g', marker='s', label='Class 2')
plt.scatter(X_train_lad[y_train == 3, 0], X_train_lad[y_train == 3, 1],
            c='b', marker='x', label='Class 3')
plt.xlabel('LDA Component 1')
plt.ylabel('LDA Component 2')
plt.title('LDA of Wine Dataset with Decision Regions')
plt.legend(loc='upper left')
plt.show()

```
</details>

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4072443/915cb1ea-50d3-466b-b3d2-f4b070093fd7.png" width="70%">


# 4. カーネル主成分分析（Kernel PCA）
現実の問題では、非線形問題を扱うことが多いが、前述のPCAやLDAは線形手法であるため、非線形問題の次元削減手法としてはベストではない。そこで、非線形変換法として、**カーネルPCA**がある。
まず、非線形問題を解くには、より高い次元の新しい特徴量空間へ射影し、そこで線形分離（PCA）を行うといったアプローチを取る。その際に用いられるのが、別記事で解説したカーネルトリックである。カーネルトリックを用いることで、高次元への写像にかかる計算コストを削減することができる。
データ点$\mathbf{x}$をより高次元の空間へ射影するための、非線形射影関数を$\phi$とする。高次元空間でPCA線形変換を行うには、高次元空間においてカーネル行列$K$の固有ベクトルを計算する必要があるが、以下のようなカーネルトリックを行うことで、その固有ベクトルを明示的に計算する必要がなくなる。
$$
K(x^{(i)}, x^{(j)}) =  \phi(x^{(i)})^T \phi(x^{(j)})
$$

[カーネルトリックの種類や詳細はこちらを参照](https://qiita.com/Djay/items/3a5cbf409d1518a5142b#32-%E3%82%AB%E3%83%BC%E3%83%8D%E3%83%ABsvm)

# 4.1 Kernel PCAの実装
ここでは、以下のような2つの半月形のデータを用いて、カーネルPCAを実装する。
<details><summary>コードを表示</summary>

```python
# 2つの半月形データを生成
from sklearn.datasets import make_moons
X, y = make_moons(n_samples=100, noise=0, random_state=0)
plt.scatter(X[y == 0, 0], X[y == 0, 1], color='red', marker='o', label='Class 0')
plt.scatter(X[y == 1, 0], X[y == 1, 1], color='blue', marker='s', label='Class 1')
plt.xlabel('Feature 1')
plt.ylabel('Feature 2')
plt.title('Two Moons Dataset')
plt.legend(loc='upper left')
plt.show()
```
</details>
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4072443/5b6b8681-5b47-47c1-a520-5d5797f479b0.png" width="50%">

まず、通常のPCAが非線形データに有効でないことを確認する。左図から、2つの主成分により新しい部分空間に射影した結果、うまく線形分離が出来ていないことがわかる。右図では、第一主成分だけ取り出しているが、同様にして線形分離が出来ていないことを示している。
<details><summary>コードを表示</summary>

```python
# 通常のPCAを実行
pca = PCA(n_components=2)
X_pca = pca.fit_transform(X)
fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(8, 4))
# PCAの結果を描画
ax[0].scatter(X_pca[y == 0, 0], X_pca[y == 0, 1], color='red', marker='o', label='Class 0')
ax[0].scatter(X_pca[y == 1, 0], X_pca[y == 1, 1], color='blue', marker='s', label='Class 1')
ax[0].set_xlabel('Principal Component 1')
ax[0].set_ylabel('Principal Component 2')
plt.legend(loc='upper left')
ax[0].set_title('PCA of Two Moons Dataset')
# 最初の主成分だけをプロット
ax[1].scatter(X_pca[y == 0, 0], np.zeros((50, 1)), color='red', marker='o', alpha=0.2, label='Class 0')
ax[1].scatter(X_pca[y == 1, 0], np.zeros((50, 1)), color='blue', marker='s', alpha=0.2, label='Class 1')
ax[1].set_xlabel('Principal Component 1')
ax[1].set_yticks([])
ax[1].set_title('PCA of Two Moons Dataset (PC1 only)')
plt.legend(loc='upper left')
plt.show()
```
</details>
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4072443/1f6f5e2c-179e-49f5-9cde-db1ee87a26e2.png" width="80%">

そこで、RBFカーネルを使ったカーネルPCAで次元削減を実行する。図を確認すると、元の非線形データの2つのクラスがうまく線形分離可能にな新しい部分空間にデータが射影されていることがわかる。
```python
# カーネルPCAを実行
from sklearn.decomposition import KernelPCA
# カーネルPCAのインスタンスを生成
kpca = KernelPCA(n_components=2, kernel='rbf', gamma=15)
X_kpca = kpca.fit_transform(X)
# カーネルPCAの結果を描画
fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(8, 4))
ax[0].scatter(X_kpca[y == 0, 0], X_kpca[y == 0, 1], color='red', marker='o', label='Class 0')
ax[0].scatter(X_kpca[y == 1, 0], X_kpca[y == 1, 1], color='blue', marker='s', label='Class 1')
ax[0].set_xlabel('Principal Component 1')
ax[0].set_ylabel('Principal Component 2')
plt.legend(loc='upper left')
ax[0].set_title('Kernel PCA of Two Moons Dataset')
# 最初の主成分だけをプロット
ax[1].scatter(X_kpca[y == 0, 0], np.zeros((50, 1)), color='red', marker='o', alpha=0.2, label='Class 0')
ax[1].scatter(X_kpca[y == 1, 0], np.zeros((50, 1)), color='blue', marker='s', alpha=0.2, label='Class 1')
ax[1].set_xlabel('Principal Component 1')
ax[1].set_yticks([])
ax[1].set_title('Kernel PCA of Two Moons Dataset (PC1 only)')
plt.legend(loc='upper left')
plt.show()
```

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4072443/fde8cb76-fccb-4afc-91ef-b3a219dfb3ff.png" width="80%">

# 5. まとめ
本記事では、主成分分析（PCA）、線形判別分析（LDA）、カーネル主成分分析（Kernel PCA）といった代表的な特徴量抽出手法について、理論とPythonによる実装の両面から解説した。これらの手法は、モデルの予測精度を向上させるだけでなく、計算コストや過学習リスクを下げる目的でも極めて有用である。

実務においては、次のような観点から特徴量抽出の活用が求められる：

- 可視化や異常検知の前処理：高次元データを2〜3次元に圧縮することで、分布の傾向や外れ値を視覚的に把握しやすくなる
- ノイズ除去と精度向上：ノイズを含んだ不要な特徴量を圧縮・統合することで、学習モデルの汎化性能が改善されるケースが多い
- 学習スピードの高速化：特徴量数を削減することで、学習・推論の計算時間が短縮され、特に大規模データを扱う現場では効率向上に直結する

特に実務では、「どの次元削減手法が最適か」はデータの構造や目的に強く依存する。まずは少数の手法を試してみて、モデルの性能や解釈性とのバランスを確認しながら進めることが重要である。

最後に、特徴量抽出はあくまで前処理手法の1つに過ぎない。モデルの選定やチューニングと組み合わせてこそ、その効果が最大化される。多次元のデータを扱う際は、特徴量抽出の活用を「当たり前の選択肢」として意識しておくとよいだろう。


```python
Python version: 3.10.4
numpy version: 2.2.5
matplotlib version: 3.10.1
scikit-learn version: 1.6.1
pandas version: 2.2.3
```