0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

機械学習入門 第6回:PCAを「情報を保ったまま次元を減らす方法」として理解する

0
Posted at

前回は、クラスタリングを使って、正解ラベルのないデータを似ているもの同士に分ける方法を見ました。

今回は、同じく教師なし学習でよく使われる 主成分分析 を扱います。英語では Principal Component Analysis、略して PCA と呼ばれます。

PCA は、ざっくり言えば「たくさんある特徴量を、情報をなるべく失わない少数の軸に言い換える方法」です。表の列が多すぎて見通しが悪いとき、2次元にして可視化したいとき、画像のような大きなデータを少し軽く持ちたいときに使われます。

ファインマン風に言えば、PCA は「立体的なものにライトを当てて、いちばん形がよく残る影の向きを探す作業」です。影にすると情報は減りますが、向きをうまく選べば、大事な形はかなり残せます。

目次

  1. PCA で何をしたいのか
  2. この記事で使うデータについて
  3. PCA の中心思想:分散が大きい方向を探す
  4. 中心化と標準化
  5. 共分散行列と固有値分解で PCA を見る
  6. NumPy で PCA を実装する
  7. SVD で PCA を見る
  8. 主成分の数をどう選ぶか
  9. Scikit-learn で PCA を使う
  10. 2次元にして可視化する
  11. PCA で画像を圧縮する
  12. PCA を使うときの注意点
  13. まとめ

1. PCA で何をしたいのか

データ分析では、1つのサンプルに対して多くの特徴量を持つことがあります。

たとえば店舗データなら、売上、来店客数、客単価、広告費、リピート率などが並びます。これらは完全に独立しているとは限りません。来店客数が多い店舗は売上も高くなりやすく、広告費が大きい店舗は来店客数も増えやすい、というように特徴量同士が関連していることがあります。

PCA は、このような関連した特徴量を、互いに重ならない新しい軸に変換します。この新しい軸を 主成分 と呼びます。

PCA の代表的な使い道は、次のようなものです。

  • 高次元データを2次元や3次元に落として可視化する
  • 似た情報を持つ特徴量をまとめ、前処理として使う
  • データの圧縮やノイズの影響を少し抑える
  • 多数の特徴量の背後にある大きな傾向を読む

ここで大切なのは、PCA は分類器でも回帰モデルでもないことです。正解ラベルを予測するのではなく、特徴量の見方を変えて、データの構造を見やすくします。

PCA の目的が見えたところで、この記事で使うデータを用意します。

2. この記事で使うデータについて

この記事では、外部ファイルに依存しないように、説明用の小さな架空データと、Scikit-learn に付属する画像データを使います。

  • PCA の仕組み、NumPy 実装、Scikit-learn 実装では、架空の日本の店舗データ store_jp を使います。
  • 画像圧縮の例では、Scikit-learn 付属の load_digits を使います。これは 8×8 ピクセルの手書き数字画像データで、追加の画像ファイルを用意しなくてもそのまま試せます。

まず、店舗データを作ります。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.preprocessing import StandardScaler


store_jp = pd.DataFrame({
    "店舗": [
        "東京駅前", "新宿西口", "横浜みなと", "名古屋栄",
        "大阪梅田", "京都河原町", "札幌駅前", "仙台中央",
        "福岡天神", "金沢片町", "松山大街道", "那覇国際通り",
    ],
    "月間売上_万円": [980, 920, 860, 760, 810, 700, 620, 590, 680, 430, 390, 520],
    "来店客数": [4200, 3900, 3600, 3100, 3350, 2800, 2500, 2400, 2900, 1700, 1500, 2100],
    "客単価_円": [2330, 2360, 2390, 2450, 2420, 2500, 2480, 2460, 2350, 2530, 2600, 2470],
    "広告費_万円": [85, 78, 70, 62, 67, 54, 45, 43, 52, 28, 24, 36],
    "リピート率_%": [48, 46, 44, 42, 43, 41, 39, 38, 40, 35, 34, 37],
})

store_features = [
    "月間売上_万円",
    "来店客数",
    "客単価_円",
    "広告費_万円",
    "リピート率_%",
]

print(store_jp)

各列の意味は次の通りです。

列名 意味
店舗 説明用に用意した仮の店舗名です。実在店舗の実績ではありません。
月間売上_万円 1か月の売上を万円単位で表した値です。
来店客数 1か月に来店した顧客数です。
客単価_円 1人あたりの平均購入金額です。
広告費_万円 1か月に使った広告費を万円単位で表した値です。
リピート率_% 再来店した顧客の割合です。単位は%です。

このデータは、説明のために作った架空データです。実在の売上データではありません。売上、来店客数、広告費、リピート率はある程度同じ方向に動くようにし、客単価だけは少し違う動きをするようにしています。PCA で「大きな傾向」と「そこから少し外れる傾向」を見やすくするためです。

データの準備ができたので、次は PCA がどのような方向を探しているのかを見ます。

3. PCA の中心思想:分散が大きい方向を探す

PCA の本質は、データを投影したときに、ばらつきが大きく残る方向を探す ことです。

「ばらつきが大きい」とは、サンプル同士の違いがよく見えるということです。反対に、投影した先で点がほとんど同じ場所に重なってしまう方向は、データの違いをあまり残せていません。

中心化されたサンプルを $x^{(i)}$、長さ1の投影方向を $u$ とします。このとき、サンプルを方向 $u$ に投影した値は次のように書けます。

$$
z^{(i)}=u^Tx^{(i)}
$$

PCA は、この $z^{(i)}$ の分散ができるだけ大きくなるような $u$ を選びます。共分散行列を $\Sigma$ とすると、目的は次のように書けます。

$$
\max_{u};u^T\Sigma u
$$

ただし、方向ベクトル $u$ の長さを勝手に大きくすれば値も大きくなってしまうので、次の制約を置きます。

$$
u^Tu=1
$$

この最大化問題の答えは、固有ベクトル の中から見つかることが知られています。固有ベクトル $u$ とは、$\Sigma u=\lambda u$ のように、$\Sigma$ を掛けても向きが変わらず、長さだけが $\lambda$ 倍になるベクトルのことです。この $\lambda$(固有値)は、ちょうどその方向にデータを投影したときの分散の大きさに対応します。つまり、一番大きい固有値を持つ固有ベクトルを選べば、そのまま分散が最大になる方向が手に入ります。こうして、共分散行列の最大固有値に対応する固有ベクトルが第1主成分になります。

第2主成分は、第1主成分と直交する方向の中で、次に分散が大きい方向です。直交する方向を選ぶのは、すでに第1主成分が拾った情報とできるだけ重ならないようにするためです。

ここでいう「情報」とは、主に分散のことです。もちろん、すべての情報が分散だけで測れるわけではありません。それでも PCA では、「大きく変わっている方向には、データを見分ける手がかりが多い」と考えます。

この考え方をコードにする前に、PCA でほぼ必ず必要になる中心化と標準化を確認します。

4. 中心化と標準化

PCA では、まず各特徴量の平均を 0 にそろえます。これを 中心化 と呼びます。中心化をしないと、データの重心が原点からずれたままになり、分散が本当に大きい方向ではなく、原点からの距離そのものが計算に影響してしまいます。

特徴量 $x_j$ の平均を $\mu_j$ とすると、中心化は次のように書けます。

$$
x_j'=x_j-\mu_j
$$

さらに、特徴量の単位やスケールが違う場合は、標準偏差で割ってスケールもそろえます。これを 標準化 と呼びます。

$$
x_j'=\frac{x_j-\mu_j}{\sigma_j}
$$

今回の店舗データでは、月間売上_万円、来店客数、客単価_円、リピート率_% のように単位がばらばらです。標準化しないまま PCA を行うと、値の大きい 来店客数 の影響が強くなりすぎます。

そこで、まず数値特徴量を標準化します。

X_store = store_jp[store_features]

scaler = StandardScaler()
X_store_scaled = scaler.fit_transform(X_store)

scaled_df = pd.DataFrame(
    X_store_scaled,
    columns=store_features,
    index=store_jp["店舗"],
)

print(scaled_df.round(2))

StandardScaler は、各列の平均がおおよそ 0、標準偏差がおおよそ 1 になるように変換します。これで、単位の違いではなく、相対的な高低を比べやすくなります。

前処理ができたので、次は共分散行列から主成分を求めます。

5. 共分散行列と固有値分解で PCA を見る

中心化または標準化済みのデータ行列を $X\in\mathbb{R}^{m\times n}$ とします。ここで、$m$ はサンプル数、$n$ は特徴量数です。

各行が1つのサンプル、各列が1つの特徴量だとすると、共分散行列は次のように書けます。

$$
\Sigma=\frac{1}{m-1}X^TX
$$

共分散行列は、「特徴量同士がどの方向に一緒に動くか」をまとめた行列です。対角成分はそれぞれの特徴量自身の分散、対角以外の成分は2つの特徴量が一緒に増減する度合い(共分散)を表します。PCA では、この共分散行列を固有値分解します。

$$
\Sigma U=U\Lambda
$$

ここで、$U$ の各列が主成分の方向、$\Lambda$ の対角成分が固有値です。固有値は、その主成分がどれくらいの分散を説明しているかを表します。

前 $K$ 個の主成分方向を集めた行列を $U_{\text{reduce}}$ とすると、低次元表現は次のようになります。

$$
Z=XU_{\text{reduce}}
$$

この形なら、$m\times n$ のデータが $m\times K$ のデータに変わります。

流れをまとめると、次のようになります。

数式の流れが見えたところで、NumPy だけで実装してみます。

6. NumPy で PCA を実装する

共分散行列は実対称行列なので、NumPy では np.linalg.eigh を使うと安定して固有値分解できます。np.linalg.eig でも計算できますが、対称行列には eigh が向いています。

cov_matrix = np.cov(X_store_scaled, rowvar=False)

eigenvalues, eigenvectors = np.linalg.eigh(cov_matrix)

order = np.argsort(eigenvalues)[::-1]
eigenvalues = eigenvalues[order]
eigenvectors = eigenvectors[:, order]

explained_ratio = eigenvalues / eigenvalues.sum()
cumulative_ratio = np.cumsum(explained_ratio)

eigen_table = pd.DataFrame({
    "固有値": eigenvalues,
    "説明分散比": explained_ratio,
    "累積説明分散比": cumulative_ratio,
}, index=[f"PC{i}" for i in range(1, len(eigenvalues) + 1)])

print(eigen_table.round(3))

説明分散比 は、各主成分が全体の分散のうち何割を説明しているかを表します(主成分分析の文脈では 寄与率 と呼ばれることもあります)。累積説明分散比 は、第1主成分から順に足し上げた割合で、累積寄与率 と呼ばれることもあります。

ここでは、まず2次元に落としてみます。

n_components = 2
U_reduce = eigenvectors[:, :n_components]
Z_np = X_store_scaled @ U_reduce

pca_np_result = pd.DataFrame(Z_np, columns=["PC1", "PC2"])
pca_np_result.insert(0, "店舗", store_jp["店舗"])

print(pca_np_result.round(3))

主成分が元の特徴量とどう関係しているかも見ておきます。このように、各主成分がどの特徴量をどれくらい反映しているかを表す値を 負荷量 と呼びます。

loading_np = pd.DataFrame(
    U_reduce,
    index=store_features,
    columns=["PC1", "PC2"],
)

print(loading_np.round(3))

この表は、各主成分がどの特徴量を強く見ているかを読むためのものです。実際に計算すると、PC1 は 月間売上_万円・来店客数・広告費_万円・リピート率_% がどれも $-0.45$ 前後とほぼ同じ大きさで並び、客単価_円 だけ符号が逆の $+0.41$ になります。つまり PC1 は「店舗規模や集客力が動くと、客単価だけは逆向きに動く」という軸だと読めます。一方 PC2 は 客単価_円 が $-0.90$ と抜けて大きく、ほかの特徴量は $-0.1$〜$-0.3$程度にとどまります。つまり PC2 は、ほぼ 客単価_円 だけを映す軸になっています。

ただし、固有ベクトルの符号は反転しても同じ方向を表します。ある環境では PC1 の符号がすべて逆に見えることがありますが、その場合も解釈の軸そのものは同じです。

共分散行列から主成分を求める方法が見えたので、次は PCA でよく使われるもう1つの見方である SVD を確認します。

7. SVD で PCA を見る

PCA は、共分散行列を作って固有値分解する代わりに、データ行列そのものを 特異値分解 して求めることもできます。英語では Singular Value Decomposition、略して SVD と呼ばれます。

標準化済みデータ行列 $X$ の SVD は、次のように書けます。

$$
X=USV^T
$$

ここで、$U$ の列は左特異ベクトル、$S$ は特異値を大きい順に対角に並べた行列、$V$ の列は右特異ベクトルです。この中で、$V$ の列が PCA の主成分方向に対応します。NumPy の np.linalg.svd は $V^T$ を返すので、先頭の行を取り出して転置すると、主成分方向として使えます。

u, singular_values, v_t = np.linalg.svd(X_store_scaled, full_matrices=False)

V_k = v_t[:2].T
Z_svd = X_store_scaled @ V_k

explained_variance_svd = singular_values ** 2 / (X_store_scaled.shape[0] - 1)
explained_ratio_svd = explained_variance_svd / explained_variance_svd.sum()

print("SVD から計算した説明分散比:")
print(np.round(explained_ratio_svd, 3))

loading_svd = pd.DataFrame(V_k, index=store_features, columns=["PC1", "PC2"])
print(loading_svd.round(3))

SVD で得られる主成分方向も、固有値分解で得られる方向と同じです。ただし、ここでも符号は反転することがあります。PC1 全体がプラスからマイナスに反転しても、同じ直線上の反対向きを見ているだけです。

SVD は、大きなデータや数値計算の実装でよく使われます。Scikit-learn の PCA も、内部では状況に応じて SVD 系の計算を使います。

主成分を求める方法がわかったところで、次は何個の主成分を残すかを考えます。

8. 主成分の数をどう選ぶか

PCA では、主成分を多く残せば元の情報に近くなります。その代わり、次元削減や圧縮の効果は小さくなります。

反対に、主成分を少なくすればデータは軽くなりますが、情報の損失は大きくなります。

よく使われる考え方は、累積説明分散比 を見ることです。固有値を大きい順に並べたものを $\lambda_1,\lambda_2,\ldots,\lambda_n$ とすると、前 $K$ 個の主成分が説明する割合は次のように書けます。

$$
R_K=\frac{\sum_{i=1}^{K}\lambda_i}{\sum_{i=1}^{n}\lambda_i}
$$

たとえば $R_K\ge0.90$ なら、前 $K$ 個の主成分で分散の約90%を説明できる、という意味です。

target_ratio = 0.90
n_required = np.searchsorted(cumulative_ratio, target_ratio) + 1

print(f"{target_ratio:.0%} 以上を説明するために必要な主成分数: {n_required}")

plt.plot(
    range(1, len(cumulative_ratio) + 1),
    cumulative_ratio,
    marker="o",
)
plt.axhline(target_ratio, color="gray", linestyle="--")
plt.xlabel("主成分の数")
plt.ylabel("累積説明分散比")
plt.title("主成分数と累積説明分散比")
plt.ylim(0, 1.05)
plt.show()

今回の店舗データでは、第1主成分だけで90%以上を説明できる結果になるかもしれません。ただし、散布図として見たい場合は横軸と縦軸が必要なので、可視化のために2主成分を残すこともよくあります。目的が「圧縮」なのか「2次元で眺めること」なのかで、選ぶ主成分数は変わります。

もう1つの見方は、投影してから元の空間へ戻したときの誤差を見ることです。2次元に落とした表現から、標準化後のデータを近似的に復元してみます。

X_reconstructed_scaled = Z_np @ U_reduce.T

relative_error = (
    np.sum((X_store_scaled - X_reconstructed_scaled) ** 2)
    / np.sum(X_store_scaled ** 2)
)

print("2主成分での相対再構成誤差:", round(relative_error, 3))
print("2主成分で残せた情報の目安:", round(1 - relative_error, 3))

相対再構成誤差が 0.10 なら、かなり大ざっぱには約90%の情報を残したと見られます。99%を残したいなら、誤差は 0.01 程度まで小さくしたい、という関係です。

なお、SVD から説明分散比を計算する場合は、特異値そのものではなく、特異値の2乗を使います。分散に対応するのは $s_i^2$ だからです。

主成分数の選び方が見えたところで、次は Scikit-learn で同じことを短く書きます。

9. Scikit-learn で PCA を使う

Scikit-learn では、PCA は PCA クラスとして用意されています。

必要なライブラリが入っていない場合は、先にインストールします。

pip install numpy pandas matplotlib scikit-learn

2次元へ落とす場合は、n_components=2 と指定します。

from sklearn.decomposition import PCA


pca = PCA(n_components=2)
Z_sklearn = pca.fit_transform(X_store_scaled)

sklearn_result = pd.DataFrame(Z_sklearn, columns=["PC1", "PC2"])
sklearn_result.insert(0, "店舗", store_jp["店舗"])

print(sklearn_result.round(3))
print("説明分散比:", np.round(pca.explained_variance_ratio_, 3))
print("累積説明分散比:", np.round(pca.explained_variance_ratio_.cumsum(), 3))

fit_transform は、主成分方向を学習し、そのまま低次元データへ変換します。

主成分方向は components_ で確認できます。Scikit-learn では、components_ の各行が1つの主成分です。見やすいように転置して表にします。

loadings = pd.DataFrame(
    pca.components_.T,
    index=store_features,
    columns=["PC1", "PC2"],
)

print(loadings.round(3))

累積説明分散比を指定して、自動で主成分数を選ぶこともできます。たとえば、90%以上の分散を残したい場合は次のように書きます。

pca_90 = PCA(n_components=0.90, svd_solver="full")
Z_90 = pca_90.fit_transform(X_store_scaled)

print("選ばれた主成分数:", pca_90.n_components_)
print("累積説明分散比:", round(pca_90.explained_variance_ratio_.sum(), 3))
print("変換後の形:", Z_90.shape)

n_components=0.90 は「90%を保つために必要な主成分数を選ぶ」という意味です。整数を指定したときは「その個数にする」、0から1の小数を指定したときは「その割合以上を保つ」という読み方になります。

Scikit-learn で短く書けるようになったので、次は2次元に落とした結果を図で見ます。

10. 2次元にして可視化する

PCA は、高次元データを2次元に落として散布図で見るときによく使われます。

plt.scatter(Z_sklearn[:, 0], Z_sklearn[:, 1], s=80)

for store_name, pc1, pc2 in zip(store_jp["店舗"], Z_sklearn[:, 0], Z_sklearn[:, 1]):
    plt.text(pc1 + 0.03, pc2 + 0.03, store_name, fontsize=9)

plt.axhline(0, color="gray", linewidth=0.8)
plt.axvline(0, color="gray", linewidth=0.8)
plt.xlabel("PC1")
plt.ylabel("PC2")
plt.title("PCA による店舗データの2次元表示")
plt.show()

図の横軸と縦軸は、元の 月間売上_万円 や 来店客数 そのものではありません。複数の特徴量を混ぜて作った新しい軸です。

そのため、PCA の散布図を見るときは、次の2つをセットで確認します。

  • 各サンプルが PC1・PC2 上でどこにあるか
  • components_ や loadings を見て、PC1・PC2 がどの特徴量を強く反映しているか

PCA の散布図は、クラスタリング結果のようにグループ名を自動で付けるものではありません。点の位置を見やすくするための地図だと考えると、使いどころがつかみやすくなります。

2次元への可視化ができたので、次は PCA を圧縮として使う例を見ます。

11. PCA で画像を圧縮する

PCA は画像圧縮の考え方にも使えます。画像は、ピクセルを並べると大きな数値ベクトルとして扱えます。PCA で重要な方向だけを残せば、少ない数で画像を近似できます。

Qiita の記事では、読者が同じ画像ファイルを持っているとは限りません。そこでここでは、Scikit-learn に付属する load_digits を使います。これは 8×8 ピクセルの手書き数字画像なので、本物の写真ほど大きくはありませんが、「画像を数値行列として扱い、主成分で近似する」という流れは同じです。

from sklearn.datasets import load_digits


digits = load_digits()
X_digits = digits.data
images = digits.images

print("画像枚数:", X_digits.shape[0])
print("1枚あたりの次元数:", X_digits.shape[1])

各画像は 8×8 ピクセルなので、1枚あたり 64 次元のベクトルです。まず 16 個の主成分だけを残して、元の画像を復元してみます。

pca_image = PCA(n_components=16)
Z_digits = pca_image.fit_transform(X_digits)
X_digits_reconstructed = pca_image.inverse_transform(Z_digits)

sample_index = 0
original_image = images[sample_index]
reconstructed_image = X_digits_reconstructed[sample_index].reshape(8, 8)

fig, axes = plt.subplots(1, 2, figsize=(6, 3))

axes[0].imshow(original_image, cmap="gray_r")
axes[0].set_title("元の画像")
axes[0].axis("off")

axes[1].imshow(reconstructed_image, cmap="gray_r")
axes[1].set_title("PCAで復元")
axes[1].axis("off")

plt.tight_layout()
plt.show()

print("元の次元数:", X_digits.shape[1])
print("圧縮後の次元数:", Z_digits.shape[1])
print("累積説明分散比:", round(pca_image.explained_variance_ratio_.sum(), 3))

inverse_transform は、PCA で低次元にしたデータを元の特徴量空間へ戻すメソッドです。完全に元通りになるわけではありません。捨てた主成分が持っていた情報はすでに失われているため、残した主成分だけを使って、近い画像を作り直していると考えます。主成分数が少ないほど、輪郭はわかっても細部がぼやけた画像になりやすくなります。

95%の分散を残すために必要な主成分数も確認できます。

pca_image_95 = PCA(n_components=0.95, svd_solver="full")
Z_digits_95 = pca_image_95.fit_transform(X_digits)

print("95%を保つための主成分数:", pca_image_95.n_components_)
print("変換後の形:", Z_digits_95.shape)

主成分数を増やすほど復元画像は元に近づきますが、圧縮率は下がります。主成分数を減らすほど軽くなりますが、輪郭や細部は失われやすくなります。

画像圧縮の例まで見たところで、最後に PCA を使うときの注意点を整理します。

12. PCA を使うときの注意点

PCA は便利ですが、何にでも効く万能な方法ではありません。

まず、PCA は線形変換です。曲がった構造をそのまま1本の直線や平面で表すのは得意ではありません。非線形な構造を見たい場合は、t-SNE や UMAP、カーネル PCA など別の方法を検討することがあります。

次に、標準化の有無で結果が大きく変わります。単位が違う特徴量を一緒に扱うなら、基本的には標準化してから PCA を行います。ただし、すべての特徴量が同じ単位で、元の分散の大きさ自体に意味がある場合は、中心化だけにすることもあります。

また、主成分の符号には本質的な意味がありません。PC1 の値が全体的に反転しても、表している軸は同じです。大切なのは、どの特徴量が同じ向きに効いているか、どの特徴量が反対向きに効いているかを読むことです。

最後に、説明分散比が高いことと、予測精度が高いことは同じではありません。PCA は特徴量のばらつきをよく残す方法であって、目的変数をよく予測する方向を直接探しているわけではありません。教師あり学習の前処理として使う場合は、最終的なモデルの評価も必ず確認します。

13. まとめ

今回の内容をまとめます。

  • PCA は、高次元データを少数の主成分に変換する教師なし学習の方法である。
  • 主成分は、データを投影したときに分散が大きく残る方向である。
  • 第1主成分は最も分散が大きい方向、第2主成分はそれと直交する中で次に分散が大きい方向である。
  • 単位やスケールが違う特徴量では、PCA の前に標準化することが重要である。
  • 共分散行列の固有値分解から主成分方向を求められる。
  • SVD を使っても PCA を求められ、実装ではこちらの見方もよく使われる。
  • 主成分の数は、累積説明分散比や再構成誤差を見ながら選ぶ。
  • Scikit-learn では PCA を使うと、次元削減、説明分散比の確認、再構成を短く書ける。
  • PCA は可視化や圧縮に便利だが、線形手法であり、説明分散比と予測性能は別物である。

PCA は、たくさんの特徴量を「少ないけれど意味のある軸」に言い換えるための方法です。

第5回のクラスタリングが「似ているものを分ける」方法だったのに対して、PCA は「データを見る座標軸を作り直す」方法です。この違いを押さえておくと、教師なし学習の道具をかなり整理して使えるようになります。

0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?