「はじめてのパターン認識」（平井有三 著）という本の「第９章 部分空間法」で説明されているCLAFIC法のRでの実行例を紹介します。
学習アルゴリズム等に関しては以下のページを参考にしてください。
(1)はじパタ第９章関連
はじめてのパターン認識 9章（前半）
はじパタ9章で飛ばされがちな部分の説明2(特異値分解、CLAFIC法)
メモ
(2)部分空間法のアルゴリズム
例題ではじめる部分空間法- パターン認識へのいざない -
CLAFIC法の実行例
1.準備
- データは、$28\times28=784$画素の手書き数字データを使用します。データは以下のページからダウンロードが可能です。（kaggle）
- 私がこれから使うデータは次の2つのデータです。（名前は私が勝手に指定しています。）
（１）学習データ$train$(42,000行785列)
データ数は42,000で、1列目にlabel（クラス）データ、残りの784列がpixel（画素）データが並べられています。
（２）テストデータ$test$(28,000行784列)
データ数は$28,000$ですが、今回は20しか使いません。ちなみにこちらはlabel（クラス）についてのデータはありません。
データの整理
- 目的変数と説明変数に分ける。
y_train<-train$label
x_train<-data.matrix(train[,-1])
x_test<-data.matrix(test)
- データベクトル（説明変数）の長さが１になるように揃える。
for(i in 1:42000){
x_train[i,]<-x_train[i,]/sqrt(x_train[i,]%*%x_train[i,])
}
for(i in 1:28000){
x_test[i,]<-x_test[i,]/sqrt(x_test[i,]%*%x_test[i,])
}
※上記のベクトルの長さを1にするプログラムは実行可能だが、書き方としてはあまり良くないらしい...わかる方がいたら教えていただけると助かります。
- データベクトルの大きさを1に揃えたので、データ$train$を再度定義する。
train<-data.frame(label=y_train,x_train)
データ(説明変数)を画像表示する関数
- x_dataには784画素のデータを入れる。
- rangeで画像表示するデータ数を指定する。
- (range/5)行５列で表示される。
view_train <- function(x_data, range) {
par(mfrow=c(length(range)/5, 5))
par(mar=c(0,0,0,0))
for (i in range) {
m <- matrix(x_data[i,], 28, 28)
image(m[,28:1])
}
}
データラベル(目的変数)を表示する関数
- y_dataにはlabel(クラス)データを入れる。
- rangeで表示するラベルデータ数を指定する。
- (range/5)行５列で表示される。
view_label <- function(y_data, range) {
matrix(y_data[range], length(range)/5, 5, byrow = TRUE)
}
20個の手書き数字データを可視化
range <- 1:20
view_train(x_train, range)
view_label(y_train, range)
データtrainをクラス別に分割
- X_**$i$**にクラス$i$のデータ行列が保存されている。
X_0<-data.matrix(subset(train,train$label==0)[-1])
X_1<-data.matrix(subset(train,train$label==1)[-1])
X_2<-data.matrix(subset(train,train$label==2)[-1])
X_3<-data.matrix(subset(train,train$label==3)[-1])
X_4<-data.matrix(subset(train,train$label==4)[-1])
X_5<-data.matrix(subset(train,train$label==5)[-1])
X_6<-data.matrix(subset(train,train$label==6)[-1])
X_7<-data.matrix(subset(train,train$label==7)[-1])
X_8<-data.matrix(subset(train,train$label==8)[-1])
X_9<-data.matrix(subset(train,train$label==9)[-1])
2.CLAFIC法の実装
相関行列
- クラスごとに相関行列$Q_i$を求める。
$$
Q_i=\frac{1}{N_i}X_i^TX_i
$$
($N_i$：クラス$i$のデータ数、$X_i$：クラス$i$のデータ行列)
q0<-1/nrow(X_0)*t(X_0)%*%X_0
- 今回、行列$Q_i$は$784\times784$型の行列である。
- 各クラスの相関行列を同様に求める。(q$i$に保存）
相関行列に関する固有値問題
$$
Q_{i}u_{ij}=\lambda_{ij}u_{ij}
$$
($\lambda_{ij}$：クラス$i$の$j$番目の固有値、$u_{ij}$：クラス$i$の$j$番目の固有値に対応する固有ベクトル)
- eigen()を使って、固有値と固有ベクトルを求める。
- valuesに固有値、vectorsに固有ベクトルが保存されている。
eig0<-eigen(q0)
str(eig0)
- $Q_0$の固有値をlm_0、固有ベクトルをu_0に保存する。
lm_0<-eig0$values
u_0<-eig0$vectors
- 各クラスの相関行列にも同様の操作を行う。
固有ベクトルが正規直交基底か？
- 固有ベクトルの内積を計算してみる。
u_q0[,1]%*%u_q0[,1];u_q0[,1]%*%u_q0[,2]
u_q3[,4]%*%u_q3[,4];u_q3[,2]%*%u_q3[,5]
- 固有ベクトルの大きさが１で、異なる固有ベクトルの内積が０になることから、正規直交基底であることがわかる。
正規直交基底の可視化
view_train(t(u_q3), 1:10)
忠実度の導入
-
忠実度を$\kappa$とする。
-
固有値の累積値を用いる。
$$
a(d_i)=\frac{\sum_{j=1}^{d_i}\lambda_{ij}}{\sum_{All}\lambda_{ij}}
$$
($\sum_{All}\lambda_{ij}$：相関行列$Q_i$の固有値の総和) -
次式を満たす$d_i$を$i$番目のクラスの部分空間の次元として採用する。
$$
a(d_i-1)\leq\kappa\leq a(d_i)
$$
忠実度を定めて次元を求める！
- $\kappa=0.95$とする。
- 固有値は降順で並べられている。
- cumsum()は累積和を求めてくれる関数である。
kpa<-0.95
a0<-cumsum(lm_q0)/sum(lm_q0)
for(i in 1:784){
adi<-a0[i]
if(kpa>adi){d0<-i}
}
d0
a0[53];a0[54]
- 上記の結果からクラス０の部分空間の次元を５４として採用する。
- 他クラスの部分空間も同様に調べると、以下のようになる。
射影行列
- クラス$i$の射影行列を$P_i$とする。
- 射影行列の性質を確認しておきましょう。
射影行列の性質
(1)べき等律：$P_i^2=P_i$
(2)対称行列：$P_i^T=P_i$
(3)要素行列(固有ベクトルの外積)$u_{ij}u_{ij}^T$のランクは１
- 射影行列は部分空間の次元数の要素行列からなる。
- 要素行列は固有ベクトルの外積（テンソル積）で求められる。
- ベクトルの外積はouter()で計算可能。
u01u01<-outer(u_q0[,1],u_q0[,1]) # 要素行列
qr(u01u01)$rank # 要素行列のランク１
- 上記のことから射影行列の性質(3)を満たすことがわかる。
- クラス０の射影行列$P_0$を求める。（先程求めた次元数がここで扱う固有ベクトル数である。）
p0<-matrix(0,nrow=784,ncol=784)
for(i in 1:54){
m<-outer(u_q0[,i],u_q0[,i])
p0<-p0+m
}
- 同様に他クラスの射影行列も求められる。
sum(t(p0)-p0) # （２）転置行列との差は０
sum(p0%*%p0-p0) # （１）2乗との差も-4×10^(-14)なので０
識別規則
$$
^\forall j\neq i\quad,\quad \vec{x}^TP_j\vec{x}<\vec{x}^TP_i\vec{x}\Rightarrow \vec{x}\in C_i
$$
識別規則のイメージ
$$
\vartheta_{i}^2=||\vec{x}||^2-\sum_{j=1}^{d_i}(\vec{x}^T\vec{u_{ij}})^2=||\vec{x}||^2-\vec{x}^TP_i\vec{x}
$$
- $\vartheta_{i}^2$が最小となるクラス$C_i$に分類する。
$\Rightarrow||\vec{x}||^2$が全クラス共通であるから、$\vec{x}^TP_i\vec{x}$が最大となるクラス$C_i$に分類すればよい！
学習データの一つを識別してみよう！
- 学習データx_train[1,]を識別してみる。
- $l$に各射影行列による長さを保存する。
l<-c()
l[1]<-x_train[1,]%*%p0%*%x_train[1,];l[2]<-x_train[1,]%*%p1%*%x_train[1,]
l[3]<-x_train[1,]%*%p2%*%x_train[1,];l[4]<-x_train[1,]%*%p3%*%x_train[1,]
l[5]<-x_train[1,]%*%p4%*%x_train[1,];l[6]<-x_train[1,]%*%p5%*%x_train[1,]
l[7]<-x_train[1,]%*%p6%*%x_train[1,];l[8]<-x_train[1,]%*%p7%*%x_train[1,]
l[9]<-x_train[1,]%*%p8%*%x_train[1,];l[10]<-x_train[1,]%*%p9%*%x_train[1,]
y_train[1];which.max(l)-1
- x_train[1,]のラベルy_train[1]は１で、$which.max(l)−1$の結果が１になってるから良さそう…
識別関数を定義してみよう！
- dataにはには784画素のデータを入れる。
- rangeで識別するデータ数を指定する。
- 下記の関数はdataを識別した結果（クラス）をベクトルで返してくれる。
classfun<-function(data,range){
y_pred<-c()
for(i in range){
l<-c()
l[1]<-data[i,]%*%p0%*%data[i,]
l[2]<-data[i,]%*%p1%*%data[i,]
l[3]<-data[i,]%*%p2%*%data[i,]
l[4]<-data[i,]%*%p3%*%data[i,]
l[5]<-data[i,]%*%p4%*%data[i,]
l[6]<-data[i,]%*%p5%*%data[i,]
l[7]<-data[i,]%*%p6%*%data[i,]
l[8]<-data[i,]%*%p7%*%data[i,]
l[9]<-data[i,]%*%p8%*%data[i,]
l[10]<-data[i,]%*%p9%*%data[i,]
y_pred[i]<-which.max(l)-1
}
return(y_pred)
}
（注意）
list()やfor文をうまく使えば関数はもう少し簡単に書くことができると思います。
再代入誤り率を求める！
- 42,000個のx_trainデータに識別関数を用いる。
- y_predに識別したクラスを保存する。
range<-1:(42000)
y_pred<-classfun(x_train,range)
クロス表
tab<-table(y_train[range],y_pred)
tab
再代入誤り率
1-sum(diag(tab))/sum(tab)
- 正答率をもう少し高くできたらな...
テストデータを識別してみよう！
- データx_testの２０個を識別関数で識別する。
range1<-1:20
y_test<-classfun(x_test,range1)
テストデータの可視化
view_train(x_test, range1)
テストデータの識別されたクラスの表示
view_label(y_test, range1)
3.まとめ
再代入誤り率をみても、１００個の手書き数字を識別しようとしたとき、約８個の手書き文字を誤識別してしまうので、識別能力はまだ低いのかなと感じます。まだ工夫の余地がありそうです...
共分散行列を用いた投影距離法も実行してましたが、似たような結果が出ました。それはまた機会があったら、実行例を載せようと思います。