1
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?

@arairuca(Iruca)

最急降下法を用いた主成分分析の実装

1
Last updated at Posted at 2026-05-28

はじめに

本稿では，最急降下法（Gradient Descent）を用いて主成分分析（Principal Component Analysis, PCA）を実現する手法について述べる．主成分分析とは，多次元データを分散が最大となる方向へ射影することにより，データの本質的な構造を低次元で表現する手法である．

通常，主成分分析は固有値分解または特異値分解（SVD）によって解析的に求められるが，本稿では最急降下法による数値最適化を用いてその第1主成分方向を推定する．これにより，最急降下法の適用範囲の広さと，主成分分析の幾何学的直感を同時に理解することを目的とする．

本稿の実装では，2次元データを対象とし，最急降下法によって分散を最大化する方向ベクトル $(a, b)$ を推定する．

(本稿は、ClaudeCodeを用いて作成した。)

以下、本稿で登場するプログラムの改良後のプログラムで作成した図である。

分析結果の図示

主成分分析.png

座標変換後

主成分分析で回転したデータ.png

分散の逆数の最小化（最急降下法）

主成分分析と最急降下法.png

分散の最大化

分散の推移.png

理論

主成分分析の原理

$N$ 個の2次元データ点 ${(x_i, y_i)},$$({i=1,2,\cdot \cdot \cdot ,N})$ が与えられたとする．データの平均を

$$
\mu_x = \frac{1}{N}\sum_{i=1}^{N} x_i, \quad \mu_y = \frac{1}{N}\sum_{i=1}^{N} y_i
$$

と定義する．

主成分分析の目的は，平均中心化されたデータ $(\boldsymbol{x}_i - \boldsymbol{\mu})$ を射影したときに分散が最大となる単位ベクトル $\boldsymbol{b}$ を求めることである．方向を表すパラメータ $(a, b)$ に対し，単位ベクトルを

\boldsymbol{b} = \frac{1}{\sqrt{a^2 + b^2}} \begin{pmatrix} a \\ b \end{pmatrix}

と定義する．

このとき，各データ点の射影値 $d_i$ は

$$
d_i = \boldsymbol{b}^{\top} (\boldsymbol{x}_i - \boldsymbol{\mu}) = \frac{a(x_i - \mu_x) + b(y_i - \mu_y)}{\sqrt{a^2 + b^2}}
$$

と表される．射影後のデータの分散は

$$
\sigma^2 = \frac{1}{N}\sum_{i=1}^{N} d_i^2 - \left( \frac{1}{N}\sum_{i=1}^{N} d_i \right)^2
$$

である．第1主成分方向は，この $\sigma^2$ を最大化するベクトル $\boldsymbol{b}$ に相当する．

誤差関数の定義

最急降下法は最小化問題として定式化されるため，分散の逆数を誤差関数として定義する．

$$
E(a, b) = \frac{1}{\sigma^2(a, b)}
$$

$E(a, b)$ を最小化することは，分散 $\sigma^2$ を最大化することと等価である．

最急降下法による最適化

最急降下法は，誤差関数の勾配方向に沿ってパラメータを逐次更新する手法である．更新則は以下のように記述される．

$$
a \leftarrow a - \alpha \frac{\partial E}{\partial a}, \quad b \leftarrow b - \alpha \frac{\partial E}{\partial b}
$$

ただし $\alpha > 0$ は学習率である．

数値微分

本実装では，解析的な偏微分式を用いず，数値微分によって勾配を近似する．前進差分を用いると，

$$
\frac{\partial E}{\partial a} \approx \frac{E(a + h,, b) - E(a,, b)}{h}, \quad \frac{\partial E}{\partial b} \approx \frac{E(a,, b + h) - E(a,, b)}{h}
$$

ここで $h$ は十分小さな微小量（本実装では $h = 10^{-5}$）である．

推定結果の回転変換

最急降下法によって推定された主成分方向 $(a, b)$ から，角度 $\theta$ を

\theta = \arctan(\frac{b}{a})

と求め，回転行列

R = \begin{pmatrix}
     \cos(-\theta) & -\sin(-\theta) \\
     \sin(-\theta) & \cos(-\theta) \\ 
    \end{pmatrix}

を用いてデータを回転することで，第1主成分を横軸に対応させた座標系へのデータ変換が得られる．

プログラム

全体コード

# 2次元主成分分析を行うプログラム（最急降下法）

import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib

# サンプルデータ
num = 10
x = np.linspace(-2, 2, num)
y = (2 * x + 2) + np.random.normal(0, 1, num)

mu_x = np.mean(x)
mu_y = np.mean(y)

h = 1.0e-5


def error(a, b):
    norm = (a**2 + b**2) ** 0.5
    # データを主成分の方向に投影
    projections = (a * (x - mu_x) + b * (y - mu_y)) / norm
    variance = np.mean(projections**2) - np.mean(projections) ** 2
    # 分散の逆数を誤差と定義する
    return 1.0 / variance


# 最急降下法の設定
iterations = 10000
alpha = 0.1
a, b = -2.0, 2.0

a_history = []
b_history = []
error_history = []
variance_history = []

# 最急降下法の実行
for _ in range(iterations):
    err = error(a, b)
    grad_a = (error(a + h, b) - err) / h
    grad_b = (error(a, b + h) - err) / h
    a -= alpha * grad_a
    b -= alpha * grad_b

    a_history.append(a)
    b_history.append(b)
    error_history.append(error(a, b))
    variance_history.append(1.0 / error(a, b))

誤差関数

誤差関数 error(a, b) では，データ点の射影値をNumPyのブロードキャストにより一括計算し，その分散の逆数を返す．

def error(a, b):
    norm = (a**2 + b**2) ** 0.5
    projections = (a * (x - mu_x) + b * (y - mu_y)) / norm
    variance = np.mean(projections**2) - np.mean(projections) ** 2
    return 1.0 / variance

最急降下法の実行

各イテレーションにおいて，数値微分により勾配を計算し，パラメータ $(a, b)$ を更新する．更新後の誤差と分散をそれぞれ履歴として記録する．

for _ in range(iterations):
    err = error(a, b)
    grad_a = (error(a + h, b) - err) / h
    grad_b = (error(a, b + h) - err) / h
    a -= alpha * grad_a
    b -= alpha * grad_b

    a_history.append(a)
    b_history.append(b)
    error_history.append(error(a, b))
    variance_history.append(1.0 / error(a, b))

可視化

誤差曲面と最急降下法の軌跡

パラメータ空間 $(a, b)$ 上の誤差曲面をコンター図で示し，最急降下法によるパラメータの収束軌跡を重ねて描画する．

n = 10
a_grid = np.linspace(-10, 10, n)
b_grid = np.linspace(-10, 10, n)
A, B = np.meshgrid(a_grid, b_grid)
Errors = np.vectorize(error)(A, B)

plt.contourf(A, B, Errors, cmap="jet", alpha=0.5)
plt.scatter(a_history, b_history, c=error_history, cmap="jet", label="最急降下法")
plt.colorbar(label="誤差")
plt.legend()
plt.xlabel("a")
plt.ylabel("b")
plt.show()

分散の推移

イテレーションに伴う分散の変化を描画することで，最急降下法が分散を単調に増加させながら収束する様子を確認する．

plt.plot(np.arange(iterations), variance_history)
plt.xlabel("iteration")
plt.ylabel("分散")
plt.title("分散の推移")
plt.show()

推定された主成分方向

元データと，推定された第1主成分方向を単位ベクトルとして重ねて描画する．

norm = (a**2 + b**2) ** 0.5
plt.scatter(x, y, label="データ", color="blue")
plt.quiver(
    0, 0, a / norm, b / norm,
    color="red", angles="xy", scale_units="xy", scale=1,
    label="主成分分析の方向",
)
plt.legend()
plt.xlabel("x")
plt.ylabel("y")
plt.show()

主成分座標への変換

推定した主成分方向を基準に，回転行列 $R$ を用いてデータを変換する．変換後のデータでは，第1主成分方向が $x'$ 軸と一致する．

theta = np.arctan2(b, a)
R = np.array([
    [np.cos(-theta), -np.sin(-theta)],
    [np.sin(-theta),  np.cos(-theta)],
])
centered = np.vstack([x - mu_x, y - mu_y])
rotated = R @ centered

plt.scatter(rotated[0], rotated[1])
plt.xlabel("x'")
plt.ylabel("y'")
plt.title("主成分分析で回転したデータ")
plt.show()

実行結果

誤差曲面と最急降下法の軌跡

パラメータ空間上の誤差曲面（コンター図）に，最急降下法による更新軌跡を重ねた図を以下に示す．初期値 $(a, b) = (-2, 2)$ から出発し，誤差の極小点へと収束する様子が確認できる．

線形回帰2乗誤差と最急降下法.png

分散の推移

イテレーションが進むにつれ，射影分散が増加し，やがて収束していることが確認できる．これは，誤差関数 $E(a, b) = 1/\sigma^2$ の最小化が分散の最大化に対応することと整合的である．

分散の推移.png

推定された主成分方向

散布図上に，最急降下法によって推定された第1主成分方向（赤矢印）を重ねた図を示す．矢印はデータの広がりが最大となる方向を指しており，線形トレンドと良く一致していることが確認できる．

主成分分析.png

主成分座標への変換

推定した主成分方向を基準に回転したデータを示す．回転後の $x'$ 軸方向にデータの分散が集中しており，第1主成分の抽出が正しく行われていることが確認できる．

主成分分析で回転したデータ.png

改良プログラム

上記のコードをさらに改良したプログラムを以下に示す。
具体的には、サンプルデータの点と基準線との偏差が分かるようにした。

# 2次元主成分分析を行うプログラム（最急降下法）

import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib

# サンプルデータ
num = 10
x = np.linspace(-2, 2, num)
y = (2 * x + 2) + np.random.normal(0, 1, num)

mu_x = np.mean(x)
mu_y = np.mean(y)

h = 1.0e-5


def error(a, b):
    norm = (a**2 + b**2) ** 0.5
    # データを主成分の方向に投影（ベクトル化）
    projections = (a * (x - mu_x) + b * (y - mu_y)) / norm
    variance = np.mean(projections**2) - np.mean(projections) ** 2
    # 分散の逆数を誤差と定義する
    return 1.0 / variance


# 最急降下法の設定
iterations = 10000
alpha = 0.1
a, b = -2.0, 2.0

a_history = []
b_history = []
error_history = []
variance_history = []

# 最急降下法の実行
for _ in range(iterations):
    err = error(a, b)
    grad_a = (error(a + h, b) - err) / h
    grad_b = (error(a, b + h) - err) / h
    a -= alpha * grad_a
    b -= alpha * grad_b

    a_history.append(a)
    b_history.append(b)
    error_history.append(error(a, b))
    variance_history.append(1.0 / error(a, b))

# -------------------------------------------------------
# グラフ 1：誤差曲面のコンター図と最急降下法の軌跡
# -------------------------------------------------------
n_grid = 60
a_grid = np.linspace(-10, 10, n_grid)
b_grid = np.linspace(-10, 10, n_grid)
A, B = np.meshgrid(a_grid, b_grid)
Errors = np.vectorize(error)(A, B)

fig, ax = plt.subplots(figsize=(7, 6))
cf = ax.contourf(A, B, Errors, levels=30, cmap="jet", alpha=0.6)
plt.colorbar(cf, ax=ax, label="誤差（分散の逆数）")
ax.scatter(a_history, b_history, c=error_history, cmap="jet",
           s=5, zorder=3, label="最急降下法の軌跡")
ax.scatter(a_history[0], b_history[0], color="white", s=80,
           edgecolors="black", zorder=5, label="開始点")
ax.scatter(a_history[-1], b_history[-1], color="black", s=80,
           edgecolors="white", zorder=5, label="収束点")
ax.legend(fontsize=9)
ax.set_xlabel("a")
ax.set_ylabel("b")
ax.set_title("誤差曲面と最急降下法の軌跡")
plt.tight_layout()
plt.savefig("主成分分析と最急降下法.png", dpi=150)
plt.show()

# -------------------------------------------------------
# グラフ 2：分散の推移
# -------------------------------------------------------
fig, ax = plt.subplots(figsize=(7, 4))
ax.plot(np.arange(iterations), variance_history, color="royalblue")
ax.set_xlabel("iteration")
ax.set_ylabel("分散")
ax.set_title("最急降下法による分散の推移")
plt.tight_layout()
plt.savefig("分散の推移.png", dpi=150)
plt.show()

# -------------------------------------------------------
# 主成分・射影の計算（グラフ 3・4 共通）
# -------------------------------------------------------
norm = (a**2 + b**2) ** 0.5
pc_dir = np.array([a / norm, b / norm])      # 主成分方向の単位ベクトル
mean   = np.array([mu_x, mu_y])

# 各データ点を主成分直線へ射影
centered = np.column_stack([x - mu_x, y - mu_y])   # shape (num, 2)
t_proj   = centered @ pc_dir                         # 各点の射影スカラー
proj_pts = mean + np.outer(t_proj, pc_dir)           # 射影点の座標

# 主成分直線の描画範囲（余白 0.5 追加）
margin = 0.5
t_line = np.array([t_proj.min() - margin, t_proj.max() + margin])
line_pts = mean + np.outer(t_line, pc_dir)

# -------------------------------------------------------
# グラフ 3：データ点・主成分直線・偏差（垂線）の可視化
# -------------------------------------------------------
fig, ax = plt.subplots(figsize=(7, 7))
ax.set_aspect("equal")

# 主成分直線
ax.plot(line_pts[:, 0], line_pts[:, 1],
        color="red", linewidth=2, zorder=2, label="第1主成分")

# 各点 → 射影点への垂線と射影点
for i in range(num):
    ax.plot([x[i], proj_pts[i, 0]], [y[i], proj_pts[i, 1]],
            color="gray", linewidth=1.0, linestyle="--", alpha=0.8, zorder=3)
    ax.scatter(proj_pts[i, 0], proj_pts[i, 1],
               color="orange", s=50, zorder=4)

# データ点・平均点
ax.scatter(x, y, color="steelblue", s=80, zorder=5, label="データ点")
ax.scatter(mu_x, mu_y, color="green", s=150, marker="*",
           zorder=6, label=f"平均点 ({mu_x:.2f}, {mu_y:.2f})")

# 凡例用ダミー
ax.plot([], [], color="gray", linestyle="--", label="偏差（垂線方向）")
ax.scatter([], [], color="orange", s=50, label="射影点")

ax.legend(fontsize=9)
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_title("主成分分析：データ点と第1主成分への射影・偏差")
plt.tight_layout()
plt.savefig("主成分分析.png", dpi=150)
plt.show()

# -------------------------------------------------------
# グラフ 4：回転後のデータ（主成分座標系）＋偏差の可視化
# -------------------------------------------------------
theta = np.arctan2(b, a)
R = np.array([
    [np.cos(-theta), -np.sin(-theta)],
    [np.sin(-theta),  np.cos(-theta)],
])
rotated = R @ centered.T   # shape (2, num)
rx, ry  = rotated[0], rotated[1]

fig, ax = plt.subplots(figsize=(8, 5))

# x' 軸（第1主成分方向）
ax.axhline(0, color="red", linewidth=1.5, linestyle="-", label="第1主成分軸（x'）", zorder=2)

# 各点から x' 軸への垂線（これが PCA が最小化する偏差）
for i in range(num):
    ax.plot([rx[i], rx[i]], [0, ry[i]],
            color="gray", linewidth=1.0, linestyle="--", alpha=0.8, zorder=3)
    ax.scatter(rx[i], 0, color="orange", s=50, zorder=4)

# 回転後データ点
ax.scatter(rx, ry, color="steelblue", s=80, zorder=5, label="回転後データ点")

# 凡例用ダミー
ax.plot([], [], color="gray", linestyle="--", label="垂直偏差（最小化対象）")
ax.scatter([], [], color="orange", s=50, label="x' 軸上の射影点")

ax.legend(fontsize=9)
ax.set_xlabel("x'（第1主成分方向）")
ax.set_ylabel("y'（第2主成分方向）")
ax.set_title("主成分座標系でのデータ分布：垂直偏差の可視化")
plt.tight_layout()
plt.savefig("主成分分析で回転したデータ.png", dpi=150)
plt.show()

まとめ

本稿では，最急降下法を用いて2次元データに対する主成分分析を実装した．分散の逆数を誤差関数として定義し，数値微分により勾配を近似することで，解析的な固有値分解を用いることなく第1主成分方向の推定が可能であることを示した．

最急降下法を主成分分析に適用することで，以下の知見が得られた．

  • 誤差関数 $E(a, b) = 1/\sigma^2$ の勾配降下によって，分散が最大となる方向への収束が確認された．
  • 分散の推移グラフより，最急降下法は単調に分散を増加させながら収束することが観察された．
  • 推定された主成分方向は，データのトレンドと良く一致した．

本手法は2次元への適用を示したが，同様のアプローチは高次元データへも拡張可能である．また，固有値分解との比較や，確率的勾配降下法（SGD）との組み合わせも今後の課題として挙げられる．

参考

1
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
Sign upLogin
1
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?