Math＆Julia #04｜多変量解析 » 主成分分析 » 次元削減 » IRIS

はじめに

Irisデータセットについて主成分分析（Principal Component Analysis, PCA）を実装します。Irisデータセットの特徴量は 4次元ですが、これを 2次元と 3次元に削減して可視化します。

次元削減の手法として 2つの行列分解（固有値分解，特異値分解）を採用します。前者は教科書的に考え方を理解するのに適しており、後者は大規模なデータ処理に適しています。

本稿で重点を置くのは、主成分軸がどのようなものかを概観することと、2つの行列分解のそれぞれについて主成分得点の計算方法を確認することの 2点です。分析にとって重要な寄与率や主成分負荷量などの指標は扱いません。

主成分軸の概観

主成分分析の核心部分には「高次元空間におけるデータの本質的な構造を低次元空間で効率的に表現する」という目標があります。この目標がどのようにして達成されるのかを、主成分軸を中心にして概観してみます。

図 1　主成分軸

中心化データ

データの分散を計算するための前処理として、特徴量ごとの平均が0になるように中心化を行います。前処理については、Zスコアによる標準化を行うこともありますが本稿では採用しません。

主成分軸

データの分散が最も大きくなる方向を主成分軸とします。この新しい座標軸は、データの特徴を最もよく表しています。これを主成分軸$\,1$として、以降は直前の主成分軸に直交しつつ次に分散が大きくなる方向を主成分軸とする、というふうに続きます。

射影

射影は、中心化データを主成分軸上に垂直に投影する操作です。この操作では、 元の高次元空間のデータ点を低次元の主成分軸上に写します。射影を計算するには、中心化データと主成分係数（主成分軸の成分）の内積をとります。あとで出てくる数式や実装で「Z = ⋯」となっている部分がそれです。

主成分得点

射影によって得られた新しい座標値を主成分得点と呼びます。主成分得点を用いて散布図を作成すると、データの構造を可視化することができます。

【計算方法】 固有値分解を用いる方法

図２ 計算過程で登場する行列

図 2は特徴量を$X\in\Reals^{8×4}$として作成しました。

1.　中心化データ

$$\tag{1}X_c=X-\bold{1}\pmb{μ}^\mathsf{T}$$

• $X\in\Reals^{n×d}$は特徴量（$n$はサンプル数、$d$は特徴量の数）
• $\pmb{μ}$は特徴量ごとの平均ベクトル
• $\bold{1}$は$\pmb{μ}$を$X$の形状に変換するために使用
• $X_{c}\in\Reals^{n×d}$は中心化データ行列（中心化された特徴量）

2.　共分散行列

中心化データ行列$X_c$から共分散行列$Σ$を求めます。共分散行列はデータの分散と相関を表しており、主成分軸の方向を決定するための基礎となります。

$$\tag{2}Σ=\frac{1}{n-1}{X_c}^\mathsf{T}{X_c}$$

• $Σ\in\Reals^{d×d}$は共分散行列（対角成分は各特徴量の分散、非対角成分は共分散）

実装ではStatistics.cov関数を使います。この関数はデフォルト（dims=1）で列ごとの共分散を計算します。

3.　固有値分解

共分散行列$Σ$を固有値と固有ベクトルに分解して、主成分軸を決定します。

$$\tag{3}Σ=VΛV^\mathsf{T}$$

• $V\in\Reals^{d×d}$は固有ベクトル行列（列は主成分軸で、その成分は主成分係数）
• $Λ$は固有値行列（対角成分は固有値、非対角成分は0）

実装ではLinearAlgebra.eigen関数を使います。

4.　固有ベクトル行列をソート

固有ベクトル行列$V$の列（主成分軸）を分散が大きい順（降順）にソートします。

実装ではBase.sortperm関数を使ってソートします。この関数はソート結果そのものではなく、並び順を指示するインデックスを返します。

5.　主成分得点

中心化データ行列$X_{c}$と固有ベクトル行列$V_{k}$の内積を計算します。

$$\tag{4}Z=X_{c}V_{k}$$

• $V_k\in\Reals^{d×k}$は固有ベクトル行列で、分散が大きい上位 $k$ 個
• $Z\in\Reals^{n×k}$は次元削減された主成分得点行列（先に行ったソートによって、1列めから順に第1主成分、第2主成分、⋯ と呼ぶことができる）

【計算方法】 特異値分解を用いる方法

1.　中心化データ

固有値分解の場合と同じ方法で、データを中心化します。

2.　特異値分解（SVD）

中心化データ行列$X_c$を特異値分解します。

$$\tag{5}X_c=UΣV^\mathsf{T}$$

• $U\in\Reals^{n×n}$は左特異ベクトル行列
• $Σ\in\Reals^{n×d}$は特異値行列
• $V\in\Reals^{d×d}$は右特異ベクトル行列（固有値分解で得られた固有ベクトル行列$V$と一致する）

実装ではLinearAlgebra.svd関数を使います。

3.　主成分得点

式$(5)$を変形すると次のようになります。$k$は分散の大きい上位$k$個を表します。

$$\tag{6}U_{k}Σ_{k}=X_{c}V_{k}$$

式$(4)$と式$(6)$の右辺を比べると、式$(6)$の左辺と右辺のどちらでも主成分得点を求められることが分かります。実は左辺の $U_{k}Σ_{k}$ で実装すると、誤差の蓄積が無いうえに計算効率やメモリ効率が高いことが知られています。しかし今回は、固有値分解に合わせて敢えて右辺の $X_{c}V_{k}$ で実装します。

実装

Julia

using RDatasets, Statistics, LinearAlgebra, CairoMakie

# 主成分得点を計算（固有値分解）
function pca_eig(X, k=3)
    Xc   = X .- mean(X, dims=1)     # 特徴量を中心化
    Σ    = cov(Xc)                  # 共分散行列
    λ, V = eigen(Σ)                 # 固有値分解
    desc = sortperm(λ, rev=true)    # 主成分軸を分散の降順に並べ替えるためのインデックス
    Z    = Xc * V[:, desc[1:k]]     # 主成分得点（k次元に削減）
end

# 主成分得点を計算（特異値分解）
function pca_svd(X, k=3)
    Xc  = X .- mean(X, dims=1)      # 特徴量を中心化
    UΣV = svd(Xc)                   # 特異値分解
    Z   = Xc * UΣV.V[:, 1:k]        # 主成分得点（k次元に削減）
end

# 主成分得点を可視化
function vis_score(pc, species, labels)
    color = [:blue3,  :green,  :firebrick3]
    shape = [:circle, :xcross, :diamond]
    m = Dict(l => (color=color[i], shape=shape[i]) for (i, l) in enumerate(labels))
    fig = Figure(size=(600,800))
    ax1 = Axis( fig[1,1], xlabel="PC1", ylabel="PC2", aspect=1)
    ax2 = Axis3(fig[2,1], xlabel="PC1", ylabel="PC2", zlabel="PC3")
    for l in labels
        i = species .== l  # 品種名の中からラベルと一致するものについて、すべてのインデックスを抽出
        scatter!(ax1, pc[i,1], pc[i,2],          label=l, color=m[l].color, marker=m[l].shape)
        scatter!(ax2, pc[i,1], pc[i,2], pc[i,3], label=l, color=m[l].color, marker=m[l].shape)
    end
    axislegend(ax2, position=:rt)
    display(fig)
    fig
end

ハンズオン

初期設定

Julia

iris = dataset("datasets", "iris")  # Irisデータセットを取得
X = Matrix(iris[:, 1:4])            # 特徴量を取得
labels = unique(iris.Species)       # 品種名を重複排除してラベル一覧を作成

固有値分解を用いる方法

Julia

fig = vis_score(pca_eig(X), iris.Species, labels)    # 主成分得点を可視化（固有値分解）
save("04-handson-eig.png", fig)

図 3　主成分分析（固有値分解）

特異値分解を用いる方法

Julia

fig = vis_score(pca_svd(X), iris.Species, labels)    # 主成分得点を可視化（特異値分解）
save("04-handson-svd.png", fig)

図 4　主成分分析（特異値分解）

符号の任意性（sign indeterminacy）

図 3と図 4を比べると固有値分解と特異値分解では、グラフの上下や左右が反転していることが分かります。これは、固有値分解や特異値分解の結果はベクトルの向き（符号）が一意に定まらず、その選択が計算アルゴリズムの実装依存になっていることから生じます。しかし、この現象には本質的な意味はなく、重要なのはクラスターの構造やデータ点同士の相対的な位置関係です。両者は数学的に等価で、分析上の解釈の仕方もまったく同じです。

おわりに

行列分解（固有値分解，特異値分解）を用いることで、僅かなコードでIrisデータセットに次元削減を施して可視化することができました。もし、理論の細部に深入りしたくなければ、多変量解析パッケージ（MultivariateStats.jlなど）の利用を検討するとよいでしょう。