100個の変数で記述されたデータがあるとしましょう。たとえば、ある製品の品質を100種類のセンサーで測定しているような状況です。100次元すべてをそのまま分析するのは容易ではありません。しかし、センサーの多くは互いに相関しており、実は「本質的に重要な方向」はずっと少ないのではないか — この直感を数学的に定式化し、データの「最も重要な方向」を体系的に見つけるのが主成分分析（PCA: Principal Component Analysis）です。
PCAは、データサイエンスにおける最も基本的かつ広く使われる次元削減手法です。その応用は驚くほど広範囲にわたります。
- 画像処理・コンピュータビジョン: 顔画像の「固有顔（eigenfaces）」による顔認識。数千ピクセルの画像を数十個の主成分で表現できます
- 金融工学: 金利の期間構造やポートフォリオのリスク要因分析。何百もの金融商品の価格変動を数本の主成分（レベル・傾き・曲率）で説明できます
- 遺伝学・バイオインフォマティクス: 数万個の遺伝子発現データから人種や疾患に関連する主要な変動パターンを抽出
- 気象学・地球科学: 気温や気圧の空間分布データからテレコネクションパターン（EOFとも呼ばれる）を同定
本記事では、PCAの理論を分散最大化と近似誤差最小化の2つの視点から厳密に導出し、それが固有値問題に帰着することを示します。さらにPythonで実装し、理論の内容を直感的に確認します。
本記事の内容
- 主成分分析の直感的な理解 — データの分散が最大となる方向を探す
- 分散最大化の定式化とラグランジュ乗数法による導出
- 近似誤差最小化からの導出（分散最大化との等価性）
- 寄与率と累積寄与率の意味
- Pythonでの実装と可視化
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- 固有値と固有ベクトル — PCAの核心は固有値分解にあります
- 分散共分散行列 — データの散らばり方を行列で表現する方法
- ラグランジュ乗数法 — 等式制約付き最適化問題の解法
主成分分析とは — 直感的な理解
データの「広がり」が情報を持つ
PCAのアイデアを理解するために、2次元のデータを想像してみましょう。たとえば、ある学校の生徒たちの「数学の得点」と「物理の得点」を散布図にプロットしたとします。数学ができる生徒は物理もできる傾向があるため、データは右肩上がりの楕円形に分布するでしょう。
この楕円の「長軸」の方向こそ、データが最も大きく広がっている方向です。長軸方向に沿ってデータを見れば、生徒間の差（成績の個人差）を最もよく捉えることができます。逆に「短軸」の方向にはデータの広がりが小さく、情報量が少ないと言えます。
PCAは、この「楕円の長軸・短軸」を見つける操作を、任意の次元に一般化したものです。$p$ 次元のデータに対して、$p$ 本の軸を「分散が大きい順」に並べ替えます。最も分散が大きい軸を第1主成分、次に大きい軸を第2主成分、…と呼びます。
射影という操作
もう少し数学的に言うと、PCAは「高次元空間のデータを、ある方向に射影する」という操作です。$p$ 次元の各データ点 $\bm{x}_i$ を単位ベクトル $\bm{w}$ に射影すると、スカラー値 $z_i = \bm{w}^\top \bm{x}_i$ が得られます。この射影値 $z_i$ は、データ点 $\bm{x}_i$ を方向 $\bm{w}$ 上の「影」に落としたもの、と考えることができます。
PCAの第1主成分は、射影値の分散 $\text{Var}(z) = \text{Var}(\bm{w}^\top \bm{x})$ を最大にする方向 $\bm{w}$ です。分散が大きいほど、データ点同士の射影上での距離が大きくなり、識別能力が高くなるからです。
ではこの最大化問題を数学的に定式化し、解を求めてみましょう。
分散最大化による定式化
問題設定
$n$ 個のデータ $\bm{x}_1, \bm{x}_2, \dots, \bm{x}_n \in \mathbb{R}^p$ が与えられているとします。ここでデータは中心化されている（平均がゼロベクトル）と仮定します。つまり次が成り立っています。
$$ \begin{equation} \bar{\bm{x}} = \frac{1}{n}\sum_{i=1}^{n}\bm{x}_i = \bm{0} \end{equation} $$
中心化されていない場合は、事前にすべてのデータから平均を引けばよいだけです。
データ行列 $\bm{X}$ を、各データを行として並べた $n \times p$ 行列とします。
$$ \begin{equation} \bm{X} = \begin{pmatrix} \bm{x}_1^\top \\ \bm{x}_2^\top \\ \vdots \\ \bm{x}_n^\top \end{pmatrix} \end{equation} $$
このとき、分散共分散行列（以下、共分散行列）は次のように書けます。
$$ \begin{equation} \bm{S} = \frac{1}{n}\bm{X}^\top\bm{X} \end{equation} $$
$\bm{S}$ は $p \times p$ の対称半正定値行列です。
第1主成分の定式化
単位ベクトル $\bm{w}_1 \in \mathbb{R}^p$（$\|\bm{w}_1\| = 1$）にデータを射影したとき、射影値は $z_i = \bm{w}_1^\top\bm{x}_i$ です。データは中心化されているので、射影値の平均もゼロになります。射影値の分散は次のように計算されます。
$$ \text{Var}(z) = \frac{1}{n}\sum_{i=1}^n z_i^2 = \frac{1}{n}\sum_{i=1}^n (\bm{w}_1^\top\bm{x}_i)^2 $$
この式を行列の言葉で書き直しましょう。$\bm{w}_1^\top\bm{x}_i$ はスカラーなので、$(\bm{w}_1^\top\bm{x}_i)^2 = \bm{w}_1^\top\bm{x}_i\bm{x}_i^\top\bm{w}_1$ と変形できます。和を取ると次のようになります。
$$ \text{Var}(z) = \frac{1}{n}\sum_{i=1}^n \bm{w}_1^\top\bm{x}_i\bm{x}_i^\top\bm{w}_1 = \bm{w}_1^\top\left(\frac{1}{n}\sum_{i=1}^n \bm{x}_i\bm{x}_i^\top\right)\bm{w}_1 = \bm{w}_1^\top\bm{S}\bm{w}_1 $$
したがって、第1主成分は次の制約付き最適化問題の解です。
$$ \begin{equation} \max_{\bm{w}_1} \quad \bm{w}_1^\top\bm{S}\bm{w}_1 \quad \text{subject to} \quad \bm{w}_1^\top\bm{w}_1 = 1 \end{equation} $$
制約 $\bm{w}_1^\top\bm{w}_1 = 1$ がなければ、$\bm{w}_1$ のノルムを大きくするだけでいくらでも分散を大きくできてしまうため、方向だけを決めるために単位ベクトルに制約します。
この問題を解く準備が整いました。次にラグランジュ乗数法を使って最適解を求めます。
ラグランジュ乗数法による導出
第1主成分の導出
ラグランジュ乗数法を適用するために、ラグランジュ関数を構成します。
$$ \begin{equation} \mathcal{L}(\bm{w}_1, \lambda) = \bm{w}_1^\top\bm{S}\bm{w}_1 – \lambda(\bm{w}_1^\top\bm{w}_1 – 1) \end{equation} $$
$\mathcal{L}$ を $\bm{w}_1$ で微分してゼロとおきます。$\bm{S}$ は対称行列なので、$\frac{\partial}{\partial \bm{w}_1}\bm{w}_1^\top\bm{S}\bm{w}_1 = 2\bm{S}\bm{w}_1$ であり、$\frac{\partial}{\partial \bm{w}_1}\bm{w}_1^\top\bm{w}_1 = 2\bm{w}_1$ です。したがって次の条件が得られます。
$$ \frac{\partial \mathcal{L}}{\partial \bm{w}_1} = 2\bm{S}\bm{w}_1 – 2\lambda\bm{w}_1 = \bm{0} $$
整理すると次の式になります。
$$ \begin{equation} \bm{S}\bm{w}_1 = \lambda\bm{w}_1 \end{equation} $$
これは共分散行列 $\bm{S}$ の固有値方程式に他なりません。$\bm{w}_1$ は $\bm{S}$ の固有ベクトルであり、$\lambda$ は対応する固有値です。
では、どの固有ベクトルを選べばよいでしょうか。最適化の目的関数（分散）に $\bm{S}\bm{w}_1 = \lambda\bm{w}_1$ を代入してみましょう。
$$ \bm{w}_1^\top\bm{S}\bm{w}_1 = \bm{w}_1^\top(\lambda\bm{w}_1) = \lambda\bm{w}_1^\top\bm{w}_1 = \lambda $$
つまり、射影値の分散は固有値 $\lambda$ そのものに等しくなります。分散を最大化したいのですから、最大固有値に対応する固有ベクトルを選べばよいことがわかります。
この結果をまとめると、第1主成分の方向は、共分散行列の最大固有値に対応する固有ベクトルであるということです。
第2主成分以降の導出
第2主成分 $\bm{w}_2$ は、第1主成分 $\bm{w}_1$ と直交するという追加条件のもとで、分散を最大化する方向です。
$$ \begin{equation} \max_{\bm{w}_2} \quad \bm{w}_2^\top\bm{S}\bm{w}_2 \quad \text{subject to} \quad \bm{w}_2^\top\bm{w}_2 = 1, \quad \bm{w}_2^\top\bm{w}_1 = 0 \end{equation} $$
ラグランジュ関数は次のようになります。
$$ \mathcal{L}(\bm{w}_2, \lambda, \mu) = \bm{w}_2^\top\bm{S}\bm{w}_2 – \lambda(\bm{w}_2^\top\bm{w}_2 – 1) – \mu(\bm{w}_2^\top\bm{w}_1) $$
$\bm{w}_2$ で微分してゼロとおくと次の式が得られます。
$$ 2\bm{S}\bm{w}_2 – 2\lambda\bm{w}_2 – \mu\bm{w}_1 = \bm{0} $$
この式の両辺に左から $\bm{w}_1^\top$ を掛けます。$\bm{w}_1^\top\bm{S}\bm{w}_2 = \bm{w}_2^\top\bm{S}\bm{w}_1 = \lambda_1\bm{w}_2^\top\bm{w}_1 = 0$（$\bm{S}$ の対称性と直交条件を使用）、$\bm{w}_1^\top\bm{w}_2 = 0$（直交条件）、$\bm{w}_1^\top\bm{w}_1 = 1$（正規化条件）より次のようになります。
$$ 0 – 0 – \mu = 0 \quad \Longrightarrow \quad \mu = 0 $$
$\mu = 0$ を代入すると、第2主成分も固有値方程式 $\bm{S}\bm{w}_2 = \lambda\bm{w}_2$ を満たすことがわかります。分散 $\lambda$ を最大化するには、$\bm{w}_1$ と直交する固有ベクトルの中で最大の固有値を選べばよく、それは2番目に大きい固有値に対応する固有ベクトルです。
帰納的に同じ議論を繰り返せば、第 $k$ 主成分は共分散行列の第 $k$ 番目に大きい固有値に対応する固有ベクトルであることが示されます。
このように、PCAの本質は共分散行列の固有値分解（スペクトル分解）に帰着するのです。対称行列のスペクトル分解が存在すること、および固有ベクトルが互いに直交することは線形代数のスペクトル定理が保証しています。
ここまで「分散最大化」の観点からPCAを導出しました。実は、まったく別の観点 —「近似誤差の最小化」— からも同じ結果が導けます。この同値性を見ていきましょう。
近似誤差最小化からの導出
低次元空間への射影と再構成
PCAのもう一つの重要な解釈は、最良の低次元近似を与えるという見方です。$p$ 次元のデータを $k$ 次元（$k < p$）の部分空間に射影し、元のデータをその射影から再構成したとき、再構成誤差を最小にする部分空間を見つける問題として定式化できます。
データ点 $\bm{x}_i$ を $k$ 本の直交基底 $\bm{w}_1, \bm{w}_2, \dots, \bm{w}_k$ が張る部分空間に射影すると、射影点は次のようになります。
$$ \begin{equation} \hat{\bm{x}}_i = \sum_{j=1}^{k}(\bm{w}_j^\top\bm{x}_i)\bm{w}_j \end{equation} $$
再構成誤差は、元のデータ点と射影点の二乗距離の総和で測ります。
$$ \begin{equation} J = \frac{1}{n}\sum_{i=1}^{n}\|\bm{x}_i – \hat{\bm{x}}_i\|^2 \end{equation} $$
再構成誤差の展開
直交基底を $p$ 本まで拡張して $\bm{w}_1, \bm{w}_2, \dots, \bm{w}_p$ を $\mathbb{R}^p$ の正規直交基底とします。すると任意のデータ点は次のように展開できます。
$$ \bm{x}_i = \sum_{j=1}^{p}(\bm{w}_j^\top\bm{x}_i)\bm{w}_j $$
射影から再構成した点は最初の $k$ 成分だけを使うので、再構成誤差は残りの $p – k$ 成分で書けます。
$$ \bm{x}_i – \hat{\bm{x}}_i = \sum_{j=k+1}^{p}(\bm{w}_j^\top\bm{x}_i)\bm{w}_j $$
したがって、直交基底の直交性を利用すると次のようになります。
$$ \|\bm{x}_i – \hat{\bm{x}}_i\|^2 = \sum_{j=k+1}^{p}(\bm{w}_j^\top\bm{x}_i)^2 $$
再構成誤差の平均は次のように整理できます。
$$ J = \frac{1}{n}\sum_{i=1}^{n}\sum_{j=k+1}^{p}(\bm{w}_j^\top\bm{x}_i)^2 = \sum_{j=k+1}^{p}\bm{w}_j^\top\bm{S}\bm{w}_j $$
ここで $\bm{S}$ は共分散行列です。
分散最大化との等価性
全分散（共分散行列のトレース）は正規直交基底の取り方によらず一定なので、次の関係式が成り立ちます。
$$ \text{tr}(\bm{S}) = \sum_{j=1}^{p}\bm{w}_j^\top\bm{S}\bm{w}_j = \sum_{j=1}^{k}\bm{w}_j^\top\bm{S}\bm{w}_j + \sum_{j=k+1}^{p}\bm{w}_j^\top\bm{S}\bm{w}_j $$
右辺の第1項は $k$ 次元部分空間に射影したときに保存される分散、第2項は再構成誤差 $J$ です。$\text{tr}(\bm{S})$ は定数なので、再構成誤差 $J$ を最小化することは、保存される分散 $\sum_{j=1}^{k}\bm{w}_j^\top\bm{S}\bm{w}_j$ を最大化することと完全に等価です。
$$ \begin{equation} \min J = \min \sum_{j=k+1}^{p}\bm{w}_j^\top\bm{S}\bm{w}_j = \text{tr}(\bm{S}) – \max \sum_{j=1}^{k}\bm{w}_j^\top\bm{S}\bm{w}_j \end{equation} $$
保存される分散を最大化する直交基底は、先ほど導出したように共分散行列の上位 $k$ 個の固有ベクトルです。したがって、近似誤差最小化の観点からも同じ解が得られます。
具体的に計算すると、再構成誤差の最小値は「捨てた」固有値の和で与えられます。
$$ \begin{equation} J_{\min} = \sum_{j=k+1}^{p}\lambda_j \end{equation} $$
ここで $\lambda_1 \geq \lambda_2 \geq \dots \geq \lambda_p$ は共分散行列の固有値を大きい順に並べたものです。
この結果は直感的にも理解できます。大きな固有値に対応する方向にデータは大きく広がっているので、その方向を保存すれば多くの情報を保持でき、小さな固有値の方向を捨てても情報の損失は小さいのです。
では、どの程度の情報が保持されるかを定量的に評価する指標を見ていきましょう。
寄与率と累積寄与率
寄与率の定義
第 $k$ 主成分の寄与率（contribution ratio, proportion of variance explained）は、その主成分が全分散のうちどれだけの割合を説明するかを示す指標で、次のように定義されます。
$$ \begin{equation} r_k = \frac{\lambda_k}{\sum_{j=1}^{p}\lambda_j} = \frac{\lambda_k}{\text{tr}(\bm{S})} \end{equation} $$
寄与率は0以上1以下の値を取り、すべての寄与率の和は1になります。第1主成分の寄与率が最も大きく、後の主成分ほど小さくなります。
直感的には、寄与率はその主成分が「データの散らばりの何パーセントを担っているか」を表しています。第1主成分の寄与率が0.8（80%）であれば、1つの方向だけでデータの散らばりの80%を捉えているということです。
累積寄与率
上位 $k$ 個の主成分による累積寄与率は、次のように定義されます。
$$ \begin{equation} R_k = \sum_{j=1}^{k}r_j = \frac{\sum_{j=1}^{k}\lambda_j}{\sum_{j=1}^{p}\lambda_j} \end{equation} $$
累積寄与率は、$k$ 個の主成分でデータの全分散のうちどれだけを説明できるかを示します。先ほど導出した再構成誤差の観点では、累積寄与率 $R_k$ は「情報の保持率」であり、$1 – R_k$ が「情報の損失率」に対応します。
主成分数の選択
PCAの実用上の重要な問題は、何個の主成分を保持するかの決定です。代表的な基準として以下があります。
- 累積寄与率の閾値: $R_k \geq 0.80$ や $R_k \geq 0.90$ を満たす最小の $k$ を選ぶ。どの閾値を使うかはデータの性質や分析の目的による
- スクリープロット: 固有値を大きい順にプロットし、急激に減少する「肘（elbow）」の位置を $k$ とする
- カイザー基準: 固有値が1以上（データを標準化している場合）の主成分のみを保持する
いずれの基準も絶対的なものではなく、データの特性や分析目的に応じて判断する必要があります。
以上の理論を踏まえ、PythonでPCAを実装して理論を確認しましょう。
Pythonによる実装と可視化
共分散行列の固有値分解によるPCA
まず、PCAの理論をそのまま実装します。(1) データの中心化、(2) 共分散行列の計算、(3) 固有値分解、(4) 主成分スコアの計算、という手順を踏みます。ここでは合成データを使って、PCAの動作を確認します。
import numpy as np
import matplotlib.pyplot as plt
# --- 合成データの生成 ---
np.random.seed(42)
n = 200 # サンプル数
# 相関のある2次元データ
mean = [3, 5]
# 共分散行列: 対角成分が分散、非対角が共分散
cov = [[2.5, 1.8],
[1.8, 1.5]]
X = np.random.multivariate_normal(mean, cov, n)
# --- Step 1: 中心化 ---
X_mean = X.mean(axis=0)
X_centered = X - X_mean
# --- Step 2: 共分散行列の計算 ---
S = X_centered.T @ X_centered / n
# --- Step 3: 固有値分解 ---
eigenvalues, eigenvectors = np.linalg.eigh(S)
# eigh は対称行列用で固有値を昇順に返すので、降順に並べ替え
idx = np.argsort(eigenvalues)[::-1]
eigenvalues = eigenvalues[idx]
eigenvectors = eigenvectors[:, idx]
# --- Step 4: 主成分スコアの計算 ---
# 主成分スコア = 中心化データ × 固有ベクトル行列
PC_scores = X_centered @ eigenvectors
# --- 寄与率 ---
contribution = eigenvalues / eigenvalues.sum()
cumulative = np.cumsum(contribution)
print("=== PCA 結果 ===")
print(f"固有値: λ1 = {eigenvalues[0]:.4f}, λ2 = {eigenvalues[1]:.4f}")
print(f"寄与率: r1 = {contribution[0]:.4f}, r2 = {contribution[1]:.4f}")
print(f"第1主成分の方向: w1 = ({eigenvectors[0,0]:.4f}, {eigenvectors[1,0]:.4f})")
print(f"第2主成分の方向: w2 = ({eigenvectors[0,1]:.4f}, {eigenvectors[1,1]:.4f})")
print(f"直交性の確認: w1·w2 = {eigenvectors[:,0] @ eigenvectors[:,1]:.2e}")
このコードの出力から、いくつかの重要なことが確認できます。
- 固有値は分散を表す: 第1固有値は第2固有値よりもはるかに大きく、データの変動の大部分が第1主成分の方向に集中していることがわかります
- 寄与率: 第1主成分だけで全分散の90%以上を説明できます。これは、2つの変数が強く相関しているため、データが本質的に1次元的であることを反映しています
- 直交性: 2つの固有ベクトルの内積はほぼゼロ（数値誤差の範囲内）であり、互いに直交していることが確認できます
主成分方向の可視化
次に、元のデータ空間上に主成分方向を可視化し、PCAが何をしているかを幾何学的に確認します。
import numpy as np
import matplotlib.pyplot as plt
# --- 合成データ（再生成） ---
np.random.seed(42)
n = 200
mean = [3, 5]
cov = [[2.5, 1.8], [1.8, 1.5]]
X = np.random.multivariate_normal(mean, cov, n)
X_mean = X.mean(axis=0)
X_centered = X - X_mean
S = X_centered.T @ X_centered / n
eigenvalues, eigenvectors = np.linalg.eigh(S)
idx = np.argsort(eigenvalues)[::-1]
eigenvalues = eigenvalues[idx]
eigenvectors = eigenvectors[:, idx]
PC_scores = X_centered @ eigenvectors
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
# (a) 元のデータと主成分方向
ax = axes[0]
ax.scatter(X_centered[:, 0], X_centered[:, 1], alpha=0.4, s=15, color="steelblue")
scale = 2.5
for k in range(2):
ax.annotate("", xy=scale*np.sqrt(eigenvalues[k])*eigenvectors[:, k],
xytext=(0, 0),
arrowprops=dict(arrowstyle="->", color=["red", "blue"][k], lw=2.5))
ax.text(scale*np.sqrt(eigenvalues[k])*eigenvectors[0, k] + 0.15,
scale*np.sqrt(eigenvalues[k])*eigenvectors[1, k] + 0.15,
f"PC{k+1} (λ={eigenvalues[k]:.2f})", fontsize=10,
color=["red", "blue"][k], fontweight="bold")
ax.set_xlabel("$x_1$ (centered)", fontsize=11)
ax.set_ylabel("$x_2$ (centered)", fontsize=11)
ax.set_title("(a) Data with principal components", fontsize=12)
ax.set_aspect("equal")
ax.grid(True, alpha=0.3)
ax.set_xlim(-5, 5)
ax.set_ylim(-5, 5)
# (b) 主成分スコア空間
ax = axes[1]
ax.scatter(PC_scores[:, 0], PC_scores[:, 1], alpha=0.4, s=15, color="steelblue")
ax.set_xlabel("PC1 score", fontsize=11)
ax.set_ylabel("PC2 score", fontsize=11)
ax.set_title("(b) Principal component scores", fontsize=12)
ax.set_aspect("equal")
ax.grid(True, alpha=0.3)
ax.set_xlim(-5, 5)
ax.set_ylim(-5, 5)
# (c) PC1のみで再構成
ax = axes[2]
X_reconstructed = PC_scores[:, 0:1] @ eigenvectors[:, 0:1].T
ax.scatter(X_centered[:, 0], X_centered[:, 1], alpha=0.2, s=10, color="gray",
label="Original")
ax.scatter(X_reconstructed[:, 0], X_reconstructed[:, 1], alpha=0.4, s=10,
color="red", label="Reconstructed (PC1)")
for i in range(0, n, 10):
ax.plot([X_centered[i, 0], X_reconstructed[i, 0]],
[X_centered[i, 1], X_reconstructed[i, 1]],
"k-", alpha=0.15, linewidth=0.5)
ax.set_xlabel("$x_1$ (centered)", fontsize=11)
ax.set_ylabel("$x_2$ (centered)", fontsize=11)
ax.set_title("(c) Reconstruction from PC1 only", fontsize=12)
ax.set_aspect("equal")
ax.grid(True, alpha=0.3)
ax.legend(fontsize=9)
ax.set_xlim(-5, 5)
ax.set_ylim(-5, 5)
plt.tight_layout()
plt.savefig("pca_2d_visualization.png", dpi=150, bbox_inches="tight")
plt.show()
この3つのグラフから、PCAの幾何学的な意味が明確に読み取れます。
-
図(a): 楕円形に広がったデータに対して、2本の主成分方向が描かれています。赤い矢印（PC1）は楕円の長軸方向、つまりデータが最も大きく広がっている方向を指しています。青い矢印（PC2）はそれに直交する短軸方向です。矢印の長さは $\sqrt{\lambda_k}$ に比例しており、分散の大きさを視覚的に表しています
-
図(b): データを主成分座標系で表示すると、PC1軸方向に大きく広がり、PC2軸方向にはあまり広がっていないことがわかります。これは元のデータの楕円を「軸に沿った」座標系に回転させた結果です。主成分座標系では相関がゼロになっている（無相関化されている）ことに注目してください
-
図(c): PC1のみで再構成した点（赤）は、元のデータ点（灰色）をPC1方向の直線上に射影した点です。黒い線分は再構成誤差（元の点と再構成点の距離）を表しています。この誤差は第2主成分方向の成分に対応し、PC2の固有値 $\lambda_2$ が小さいほど再構成精度は高くなります
高次元データへの適用 — アイリスデータセット
PCAの真価は高次元データに適用したときに発揮されます。scikit-learnのアイリスデータセット（4次元）に対してPCAを適用し、2次元に削減して可視化します。
import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import load_iris
# --- データの読み込みと中心化 ---
iris = load_iris()
X = iris.data # (150, 4): 4つの特徴量
y = iris.target # クラスラベル (0, 1, 2)
feature_names = iris.feature_names
target_names = iris.target_names
X_mean = X.mean(axis=0)
X_centered = X - X_mean
n = X.shape[0]
# --- 共分散行列の固有値分解 ---
S = X_centered.T @ X_centered / n
eigenvalues, eigenvectors = np.linalg.eigh(S)
idx = np.argsort(eigenvalues)[::-1]
eigenvalues = eigenvalues[idx]
eigenvectors = eigenvectors[:, idx]
# --- 主成分スコア ---
PC_scores = X_centered @ eigenvectors
# --- 寄与率 ---
contribution = eigenvalues / eigenvalues.sum()
cumulative = np.cumsum(contribution)
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
# (a) スクリープロット
ax = axes[0]
ax.bar(range(1, 5), contribution, color="steelblue", alpha=0.7, label="Individual")
ax.plot(range(1, 5), cumulative, "ro-", markersize=8, label="Cumulative")
for i, (c, cum) in enumerate(zip(contribution, cumulative)):
ax.text(i + 1, c + 0.01, f"{c:.1%}", ha="center", fontsize=9)
ax.text(i + 1, cum + 0.02, f"{cum:.1%}", ha="center", fontsize=9, color="red")
ax.set_xlabel("Principal Component", fontsize=11)
ax.set_ylabel("Proportion of Variance", fontsize=11)
ax.set_title("(a) Scree Plot", fontsize=12)
ax.set_xticks(range(1, 5))
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3, axis="y")
ax.set_ylim(0, 1.15)
# (b) PC1 vs PC2
ax = axes[1]
colors = ["#e41a1c", "#377eb8", "#4daf4a"]
for k in range(3):
mask = y == k
ax.scatter(PC_scores[mask, 0], PC_scores[mask, 1],
c=colors[k], alpha=0.6, s=30, label=target_names[k])
ax.set_xlabel(f"PC1 ({contribution[0]:.1%})", fontsize=11)
ax.set_ylabel(f"PC2 ({contribution[1]:.1%})", fontsize=11)
ax.set_title("(b) PC1 vs PC2", fontsize=12)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
# (c) 因子負荷量（loadings）
ax = axes[2]
loadings = eigenvectors[:, :2] * np.sqrt(eigenvalues[:2])
for i, name in enumerate(feature_names):
ax.annotate("", xy=(loadings[i, 0], loadings[i, 1]), xytext=(0, 0),
arrowprops=dict(arrowstyle="->", color="red", lw=1.5))
ax.text(loadings[i, 0] * 1.12, loadings[i, 1] * 1.12,
name.replace(" (cm)", ""), fontsize=9, ha="center")
circle = plt.Circle((0, 0), 1, color="gray", fill=False, linestyle="--", alpha=0.5)
ax.add_patch(circle)
ax.set_xlabel(f"PC1 loading", fontsize=11)
ax.set_ylabel(f"PC2 loading", fontsize=11)
ax.set_title("(c) Loading plot (biplot arrows)", fontsize=12)
ax.set_xlim(-1.5, 1.5)
ax.set_ylim(-1.5, 1.5)
ax.set_aspect("equal")
ax.grid(True, alpha=0.3)
ax.axhline(0, color="k", linewidth=0.5)
ax.axvline(0, color="k", linewidth=0.5)
plt.tight_layout()
plt.savefig("pca_iris.png", dpi=150, bbox_inches="tight")
plt.show()
print("=== 固有値と寄与率 ===")
for i in range(4):
print(f"PC{i+1}: λ = {eigenvalues[i]:.4f}, "
f"寄与率 = {contribution[i]:.4f}, "
f"累積 = {cumulative[i]:.4f}")
アイリスデータセットへのPCA適用結果から、以下の重要な知見が得られます。
-
スクリープロット（図a）: 第1主成分だけで全分散の約72%を説明し、第2主成分を加えると約95%に達します。つまり、4次元データの情報の大部分は2次元で捉えられるということです。第3・第4主成分は合わせても5%程度しか寄与しておらず、これらの方向にはほとんど情報が含まれていません
-
主成分スコアプロット（図b）: 2次元に削減してもなお、3つの品種が視覚的に分離されていることが確認できます。特にsetosaは他の2品種から明確に分離されており、versicolorとvirginicaは一部重なっていますが大まかな分離は維持されています。PCAは教師なしの手法（クラスラベルを使わない）であるにもかかわらず、品種の違いをよく捉えています
-
因子負荷量プロット（図c）: 各原変数が主成分にどれだけ寄与しているかを矢印で表しています。petal lengthとpetal widthは互いに近い方向を向いており、PC1に強く寄与しています。sepal widthはPC2方向に寄与しています。このプロットから、PC1は花弁のサイズ、PC2はがく片の幅に対応する軸と解釈できます
再構成誤差の定量的検証
最後に、理論で導いた「再構成誤差 = 捨てた固有値の和」という関係を数値的に確認します。
import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import load_iris
# --- データ準備 ---
iris = load_iris()
X = iris.data
X_mean = X.mean(axis=0)
X_centered = X - X_mean
n, p = X_centered.shape
S = X_centered.T @ X_centered / n
eigenvalues, eigenvectors = np.linalg.eigh(S)
idx = np.argsort(eigenvalues)[::-1]
eigenvalues = eigenvalues[idx]
eigenvectors = eigenvectors[:, idx]
# --- 各 k について再構成誤差を計算 ---
reconstruction_errors = []
theoretical_errors = []
for k in range(p + 1): # k = 0, 1, 2, 3, 4
if k == 0:
X_reconstructed = np.zeros_like(X_centered)
else:
W_k = eigenvectors[:, :k]
PC_scores_k = X_centered @ W_k
X_reconstructed = PC_scores_k @ W_k.T
# 実際の再構成誤差
error = np.mean(np.sum((X_centered - X_reconstructed) ** 2, axis=1))
reconstruction_errors.append(error)
# 理論的な再構成誤差（捨てた固有値の和）
theoretical_error = np.sum(eigenvalues[k:])
theoretical_errors.append(theoretical_error)
fig, ax = plt.subplots(figsize=(8, 5))
ks = range(p + 1)
ax.plot(ks, reconstruction_errors, "bo-", markersize=8, linewidth=2,
label="Actual reconstruction error")
ax.plot(ks, theoretical_errors, "r^--", markersize=8, linewidth=1.5,
label=r"Theoretical: $\sum_{j>k}\lambda_j$")
for k in ks:
ax.text(k + 0.08, reconstruction_errors[k] + 0.02,
f"{reconstruction_errors[k]:.4f}", fontsize=9)
ax.set_xlabel("Number of components $k$", fontsize=12)
ax.set_ylabel("Mean reconstruction error", fontsize=12)
ax.set_title("Reconstruction error vs number of principal components", fontsize=13)
ax.set_xticks(ks)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("pca_reconstruction_error.png", dpi=150, bbox_inches="tight")
plt.show()
print("=== 再構成誤差の比較 ===")
for k in ks:
print(f"k={k}: 実際 = {reconstruction_errors[k]:.6f}, "
f"理論 = {theoretical_errors[k]:.6f}, "
f"差 = {abs(reconstruction_errors[k] - theoretical_errors[k]):.2e}")
この検証結果から、以下のことが確認できます。
-
理論値と実測値が完全に一致: 全ての $k$ について、実際に計算した再構成誤差（青丸）と理論的な予測値 $\sum_{j>k}\lambda_j$（赤三角）がぴたりと重なっています。数値的な差は浮動小数点演算の精度の範囲内であり、導出の正しさが実証されています
-
$k$ を増やすと再構成誤差は単調に減少: $k = 0$（全く主成分を使わない）のとき誤差は全分散 $\text{tr}(\bm{S})$ に等しく、$k = p$（全ての主成分を使う）のとき誤差はゼロになります。固有値の大きい主成分を先に加えるため、最初の数個の主成分で誤差が大幅に減少します
-
次元削減の効果: アイリスデータでは $k = 2$ でほとんどの分散が説明されるため、4次元のうち2次元を捨てても再構成誤差は非常に小さくなります。これが、先ほどの2次元プロットで品種の分離がよく見えた理由です
まとめ
本記事では、主成分分析（PCA）の理論を2つの視点から導出しました。
- 分散最大化: 射影値の分散を最大にする方向を求める制約付き最適化問題を、ラグランジュ乗数法で解くと、共分散行列の固有値方程式 $\bm{S}\bm{w} = \lambda\bm{w}$ に帰着します。最大固有値の固有ベクトルが第1主成分です
- 近似誤差最小化: 低次元部分空間への射影による再構成誤差を最小化する問題も、同じ固有値問題に帰着します。再構成誤差の最小値は捨てた固有値の和で与えられます
- 2つの視点は等価: 全分散は一定なので、保存される分散の最大化と再構成誤差の最小化は等価な問題です
- 寄与率: 固有値は各主成分の分散を表し、寄与率 $r_k = \lambda_k / \text{tr}(\bm{S})$ は各主成分の説明力を示します
- Pythonによる実装: 共分散行列の固有値分解を直接計算することでPCAを実装し、理論を数値的に検証しました
PCAは線形の手法であり、非線形な構造を持つデータに対しては限界があります。次のステップとして、PCAの確率的な拡張である因子分析や、非線形次元削減手法について学ぶとよいでしょう。
次のステップとして、以下の記事も参考にしてください。
- 因子分析の理論 — PCAと似ているが「潜在変数モデル」として定式化される手法
- 固有値と固有ベクトル — PCAの数学的基盤
- 分散共分散行列 — データの散らばり方を表現する行列