固有値問題版フィッシャーの線形判別分析

機械学習の分類アルゴリズムとして先駆的な手法であり*1、例として用いたirisデータセットを世に広めたフィッシャー御大の線形判別分析（LDA）の説明用のコード*2を探したのですが、教科書の説明と対応のよい固有値の計算をしていたのが見つからなかったので書きました。

概説

フィッシャーのLDAは、複数の群によるデータから、データから群を予測する分類器を構成する手法です*3。多変量を一変数に変換する写像（ $a$ ）をつくり、一次元の値からどの群に含まれるかを非確率的に予測します。予測するデータの写像の値を、群ごとの平均値の写像の値と比較し、もっとも近い平均値の写像の値の群を予測値とします。

写像をどう作るのかが問題になるわけですが、データ全体の共分散 $S$ を、群ごとの平均値の共分散（ $B$ ; Between-class co-variance）と、群ごとの写像の共分散（ $W$ ; Within-class co-variance）に、標本分散の計算と似たような方法で分離します。

$nS = \sum^g_{j=1}\sum^{n_j}_{i=1}(x_{i,j}-\bar{x})(x_{i,j}-\bar{x})^T = nW + nB$

$W = \frac{1}{n} \sum^G_{j=1}\sum^{n_j}_{i=1}(x_{i,j}-\hat{\mu}_j)(x_{i,j}-\hat{\mu}_j)^T$

$B = \frac{1}{n} \sum^G_{j=1} n_j (\hat{\mu}_j - \bar{x})(\hat{\mu}_j - \bar{x})^T$

$n$ が全体のサンプルサイズ、 $G$ が群の集合、 $n_j$ が群 $j$ のサンプルサイズ、 $x_{i,j}$ が観測値のベクトル、 $\mu$ が群ごとの平均ベクトル、 $\bar{x}$ がサンプル全体の平均ベクトルです。

次に、データ $x_i$ の写像 $a^Tx_i$ の分散 $\mbox{VAR}(a^Tx_i)$ を、平均値の写像分散 $a^TBa$ と、群ごとの写像の分散 $a^TWa$ に分けて考え、

$\max_{a} \frac{a^TBa}{a^TWa}$ ，

を満たす $a$ を、構成する分類器の写像として使います。ある群の平均値と別の群の平均値が、群の中の散らばりと比較して、なるべく離れるようにするわけです。この最適値を求める計算は、 $W^{-1}B$ の最大固有値に対応する固有ベクトルを求める固有値問題として考えることができます*4。

$a$ は要素の比が重要で、大きさは任意です。 $a^TWa=1$ と置くことができます。さらに $\alpha = W^\frac{1}{2}a$ と置いて、最大化問題を書き直します。 $W^{\frac{1}{2}} (W^{\frac{1}{2}})^T = W$ です。

$\max_{\alpha} \alpha^T (W^{-\frac{1}{2}})^T B W^{-\frac{1}{2}} \alpha \quad \mbox{s.t.} \ \alpha^T \alpha = 1$

$B$ は実対称行列なので、直交行列 $V$ と固有値が並んだ対角行列 $\Lambda$ を使って表すことができます。

$B = V \Lambda V^T$

$\beta = V^T W^{-\frac{1}{2}} \alpha$ と置くと

$\alpha^T (W^{-\frac{1}{2}})^T B W^{-\frac{1}{2}} \alpha = \alpha^T (W^{-\frac{1}{2}})^T V \Lambda V^T W^{-\frac{1}{2}} \alpha = \beta^T \Lambda \beta = \sum^n_{i=1}\lambda_i \beta_i^2$

となります。 $\lambda_i$ は固有値になります。 $\lambda_1 \geq \lambda_2 \geq \cdots$ と置き、（順番は任意にとれるので） $\Lambda$ の対角成分とします。

直交行列の逆行列が転置行列になることに注意すると、

$\sum^n_{i=1}\beta_i^2 = \beta^T \beta = \alpha^T (W^{-\frac{1}{2}})^T V V^T W^{-\frac{1}{2}} \alpha = \alpha ^T \alpha = 1$

となります。この制約の下で $\sum^n_{i=1}\lambda_i \beta_i^2$ の最大化を考えると、 $\beta_1^2$ の係数 $\lambda_1$ が一番大きいので、 $\beta_1=1$ で $\beta_i=0\ (i>1)$ です。よって、

$V^T W^{-\frac{1}{2}} \alpha = \beta$
$\alpha = (W^{-\frac{1}{2}})^T V \beta$

となります。 $(W^{-\frac{1}{2}})^TV\beta$ は $(W^{-\frac{1}{2}})^TV$ の最初の行があらわすベクトルになるので、 $\alpha$ は最大の固有値 $\lambda_1$ に対応する固有ベクトルになります。

$\beta_i = 1$ で $\beta_{-i} = 0$ と置くと、 $\alpha$ は任意の固有ベクトルとなります。このとき $W^{-1}Ba$ を整理すると、

$W^{-1}Ba = W^{-1}BW^{-\frac{1}{2}}\alpha = W^{-\frac{1}{2}} (W^{-\frac{1}{2}})^T V \Lambda V^T W^{-\frac{1}{2}} \alpha = W^{-\frac{1}{2}} ((W^{-\frac{1}{2}})^T V \Lambda V^T W^{-\frac{1}{2}}) \alpha$

となります。 $\alpha$ が $(W^{-\frac{1}{2}})^T V \Lambda V^T W^{-\frac{1}{2}}$ の固有ベクトルであることに注意して、

$W^{-1}Ba = W^{-\frac{1}{2}} \lambda_i \alpha = \lambda_i (W^{-\frac{1}{2}} \alpha) = \lambda_i a$

となり、 $a$ が $W^{-1}B$ の固有値 $\lambda_i$ に対応する固有ベクトルになることが分かります。 $W^{-1}B$ と $(W^{-\frac{1}{2}})^T B W^{-\frac{1}{2}}$ の行列のランクは等しいので、この二つの行列の固有値はすべて一致します。 $a$ が $W^{-1}B$ の最大の固有値に対応する固有ベクトルとなるとき、 $\alpha$ も $(W^{-\frac{1}{2}})^T B W^{-\frac{1}{2}}$ の最大の固有値に対応する固有ベクトルになっており、求める最大化条件が満たされていることが分かります。

ソースコード

Irisデータセットを使って実際に分類してみましょう。訓練データと検証データは分けた方がよいので、行番号が奇数のデータを訓練に、偶数のデータを検証に使います。

# 3種すべてを分類（1つは減らせる）
G <- c("setosa", "virginica", "versicolor")
# がくと花弁の長さで分類（2つ以外にするとプロットがエラーになる）
column <- c("Sepal.Length", "Petal.Length")

# 分類器（写像と平均値）の学習
training <- subset(iris[seq(1, nrow(iris), 2), ], Species %in% G)
X <- training[, column] 
n <- nrow(X)
mu <- apply(X, 2, mean)
Mu <- list()
B <- W <- matrix(0, ncol(X), ncol(X))
for(g in G){
  X_s <- subset(training, Species == g)[, column]
  mu_s <- apply(X_s, 2, mean)
  Mu[[g]] <- mu_s
  B <- B + (mu_s - mu) %*% t(mu_s - mu) * nrow(X_s) / n
  X_c <- apply(X_s, 1, \(x) x - mu_s)
  W <- W + X_c %*% t(X_c) / n
}
e <- eigen(solve(W) %*% B)
a <- e$vectors[, 1]

# 検証データの分類
test_data <- subset(iris[seq(2, nrow(iris), 2), ], Species %in% G)
test_data$Species <- factor(test_data$Species, labels = G)
Species_id <- as.numeric(test_data[, "Species"])
X <- test_data[, column] 
z <- c(as.matrix(X) %*% a)
z_d <- matrix(NA, length(z), length(G))
z_mu <- sapply(1:length(G), \(i) a %*% Mu[[G[i]]])
p_g <- G[t(apply(replicate(length(G), z), 1, \(z){
  z_d <- sqrt((z - z_mu)^2)
  which(min(z_d)==z_d)
}))]

# 精度の計算
accuracy <- sum(p_g == test_data$Species)/length(p_g)

# プロット（2変数のときのみ動く）
pch <- c(21, 23, 24)[1:length(G)]
bg <- c("white", "red", "blue")[1:length(G)]
plot(as.formula(sprintf("%s ~ %s", column[2], column[1])), data = test_data, pch = pch[Species_id], bg = bg[Species_id], main = "Fisher's linear discriminant")

G_o <- G[order(z_mu)]
for(i in 1:(length(G)-1)){
  c <- c(a %*% Mu[[G_o[i]]] + a %*% Mu[[G_o[i + 1]]])/2
  curve({ -(a[1]/a[2])*x + c/a[2] }, 1, 9, add = T)
}

legend("bottomright", G, pch = pch, pt.bg = bg, bty = "n")

*1:深層学習が一般化していますし、それよりは簡便な手法であるランダムフォレストやSVM（Rで機械学習（SVM） - 餡子付゛録゛）やQDAも、LDAよりは柔軟な予測器を構成できます。

*2:計算や図示はMASS::ldaやklaR::partimatのような既存パッケージの関数で間に合います。

*3:係数の標準誤差などが必要な計量分析には使えない一方、ロジットモデルなどと比較してデータが完全分離する場合でも一意な値が得られる利点があります。

*4:詳しくは8.3 Fisher’s linear discriminant rule | Multivariate Statisticsを参照してください。導出までの計算は教科書的な線形代数の操作ですが、実対称行列の性質などの復習機会になりました。