PCAとその発展(カーネルPCA、欠測値対応、latent variable-PCA,GP-LVM)

  • 参考
  • Rmdで書いてみる
  • PCAとその周辺のために、Rのパッケージkernlab,pcaMethods,dmtを使って、処理をなぞってみる
  • 基本的には固有値分解という、よいライブラリが準備されていて解釈も簡単な線形代数計算に帰着させよう、という手法で、一番基本的ないわゆるPCAを基本に、カーネル関数活用して非線形対応するという方針の下、複数のカーネル関数を用意してバリエーションを持たせ、欠測値対応するという方針の下、補間のアルゴリズムをオプションメニュー化して、バリエーションを持たせ、latent-variableでデータセットに複雑な構造を想定しようという方針の下、そのモデルをメニュー化してバリエーションを持たせる、という、形になっているのだが、その「メニュー化・バリエーション化」の点を理解するには、「実行関数の基本形はこうで、それの実行オプションはこう」という確認の仕方をするのは、よいやり方なので、その方針に沿って、整理してみたもの


kindle

---
title: "PCAとその拡張版、カーネルPCA、latent variable-PCA、欠測値対応PCA、GP-LVNM"
author: "ryamada"
date: "Thursday, April 09, 2015"
output: html_document
---



# はじめに〜基本となるPCAとその諸工夫について

一般的に、PCA(Principal Component Analysis)は、多次元データが持つ粗密・相互関連の構造を固有値分解を利用して、重要な軸とそうでない軸とに分配しなおして、データ構造を理解する方法である。

もっとも基本的なPCA(basic-PCAと呼ぶことにする。通常は、PCAと言えば、この手法のことである)は、サンプル数Ns x 観測項目数 Np の行列から、サンプル同士の遠近関係を分散共分散行列Mで表し、その固有値・固有ベクトルを求めることによって実行する。

このbasic-PCAを基本として、いくつかの工夫がなされる。

工夫には大きく分けて3方向の工夫がある。

1. basic-PCAは線形代数計算のみでできているので、「線形」な構造のみが見出されるが、「非線形」な構造を取り出すための工夫をする。カーネルPCAがこれに相当する

2. basic-PCAは線形代数演算をしているが、欠測値があるとその演算が成立しないので、欠測値に対応するための工夫をする。これが欠測値対応PCAに相当する。
3. (与えられたデータに)潜む(latentな)変数があり、その影響の仕方が異なる亜群によって全体が構成されていると考えると、basic-PCAでは描けない構造を見出すことができる。これがlatent variable-PCAに相当する。

これらの工夫を単独で行うこともあれば、組み合わせて行うこともある。

以下、Rのパッケージを利用して、それぞれについて確認する。

# 使用パッケージの導入

basic-PCAは基本関数prcomp()でも実行できるが、欠測値PCAのパッケージpcaMethodsを使うと、欠測値対応手法と同列で、欠測値に不対応のbasic-PCAが実行できる(prcomp()関数と同じ結果になる))ので、この文書では、pcaMethodsパッケージの関数を使うものとする。

欠測値対応PCAには複数の手法があるが、それをpcaMethodsパッケージでは手法オプション選択として選べるので、それを用いる。

カーネルPCAにはkernlabパッケージを用いる。

latent variable-PCAには、dmtパッケージを用いる。

以下のコマンドでそれらを読み込むことができる。依存関係のあるパッケージの必要なパッケージについては、コメントアウトした行で、それを示してある。

```{r}
# install.packages("kernlab","knitr","digest","dmt")
# source("http://bioconductor.org/biocLite.R")
# biocLite("pcaMethods")
library(knitr)
library(kernlab)
library(pcaMethods)
library(dmt)
```

# データセットの準備

データセットはlatent variable-PCA以外はirisを使う。
latent variable-PCAでは、構造が見出しやすいデータ構成のデータセットを使う方がわかりやすいので、dmtパッケージに付属のデータセットを利用する。

このデータは150サンプルが、4つの量的変数の値を持ち、品種というカテゴリ情報を持つデータである。

4変数データをdata.mとし、品種をtypeとする。

一部、値が欠損しているデータをmissing.mとして作成する。

120サンプルの品種情報を、training.typeとする。

```{r}
data(iris)
Ns <- nrow(iris) # サンプル数
Np <- ncol(iris)-1 # 変数の数
data.m <- as.matrix(iris[,1:Np]) # 行列オブジェクトに換える
type <- iris[[5]]

# rate.m : missing rate
rate.m <- 0.05
missing.m <- data.m
missing.cell <- sample(1:length(missing.m),length(missing.m)*rate.m)
missing.m[missing.cell] <- NA
```

# いわゆるPCA (basic-PCA)

サンプルが、分散共分散行列で指定されるNp次元正規分布であるとしたときに、分散の分配軸の取り方を工夫して、正規直交基底を取り直した結果が得られる。

データを正規化し、その分散共分散行列を固有値分解する。

ちなみに、分散共分散行列とは、すべての正規化されたサンプルベクトルのペアの内積を取ったものである。ここで内積云々と書くのは、カーネルPCAでこの内積を別のものに替えるというステップを入れるのだが、それとの対応がわかりやすいようにするためであるが、今は、気にしなくてもよい。

## パッケージpcaMethodsの関数pca()を用いてでbasic-PCA

basic-PCAを実施するのは、簡単である。
そのための関数を実行すればよい。
pca()関数のmethodオプションで"svd"を指定する。

```{r}
# library(pcaMethods)
basic.out <- pca(data.m, method="svd", nPcs=Np)
plot(as.data.frame(scores(basic.out)), col=as.integer(type) + 1,main="pca() basic pca")
#plotPcs(basic.out, col=as.integer(type) + 1)
```

## 手作業でbasic-PCA

手計算によって、上記の処理を確認しておく。

データの各変数の平均を0に補正する。
正規化した行列をMとすれば、PCAとは、$M M^{T}$を固有値分解の結果である。

各軸の値の正負が逆転している(かもしれない)点を除けば、同じ結果が得られる。

```{r}
preped.m <- prep(data.m)
apply(preped.m,2,mean)
# apply(preped.m,2,sd) # 標準偏差が1になるように補正することもあるが、prep()関数のデフォルトでは、それは実行しない
K <- preped.m %*% t(preped.m)
basic.out2 <- eigen(K)
plot(as.data.frame(basic.out2[[2]])[,1:Np], col=as.integer(type) + 1,main="manual basic pca")
```

# kernel-PCAで非線形に対応

## kernlabパッケージのkpca()関数でkernel-PCA

まずは、RにあるカーネルPCA用の関数で実行してみる。

```{r}
#library(kernlab)
gamma <- 0.1 # 色々な値を指定できるが、ここでは0.1を使う
kernel.out <- kpca(data.m,kernel="rbfdot",kpar=list(sigma=gamma),features=Np)
plot(as.data.frame(pcv(kernel.out)), col=as.integer(type) + 1)
```

## 手作業でkernel-PCA

前項のbasic-PCAではデータから作成した正方行列$K=M M^{T}$の固有値分解をした。
行列Mはデータを標準化したものであった。

カーネルPCRでは、$K=M M^T$の作り方を変えて、データからMと同サイズ(サンプル数xサンプル数)の正方行列を作り、それを固有値分解してやる。

その心は以下のように説明される。

$K=M M^T$の各要素は、サンプルごとのデータベクトルから得られるスカラー値($f(x_i,x_j)$)になっている。
$M M^T$では、この関数f()が、ベクトルの内積であると言い換えられる。

そこで、何か都合のよい関数$g(x_i,x_j)$を持ってきて、別の$K_g=(g(x_i,x_j))$を作ってやって、それを固有値分解することにする。

ただし、$K=M M^T$を固有値分解することが好都合だったので、結果が好都合になるような制約を課すことにする、というのはよい案だろうから、そのようにすることにする。

そのよい案の一つが次のような$g(x_i,x_j)$の取り方であり、そのようにしてできた$Kg=(g(x_i,x_j))$を少し変形(標準化する)したものである。

$$
g(x_i,x_j) = exp(-\gamma ||x_i-x_j||^2)
$$

ちなみに、basic-PCAで用いた内積は、2つのベクトルが一致しているときに大きくなり、反対を向いているときに小さくなるという性質があり、おおざっぱに言って、2つのベクトルが表す点の遠近情報(近ければ大きく、遠ければ小さい)になっている。
カーネルPCAで使う、この$g$でも、2つのベクトルが近い点を表していれば大きな値をとり、遠い点であれば小さな値を取る、という関係にある。

この$g$を使ってできた$Kg$を少し変形する、というのは、標準化の処理であって列和も行和も1にする。

計算自体は行列演算であり、以下のようにする。

Kの要素をすべて1/Ns倍して、そのうえで
$$
Kg' = 1/Ns Kg\\
Kg'' = Kg'-A Kg' - Kg' A + A Kg' A
$$
ただし、$A$はすべての成分が$1/Ns$ (Ns:要素数)という行列である。


手作業でのカーネルpCAを実行する。

```{r}
# ベクトル間の距離を出す
D <- as.matrix(dist(preped.m))
gamma <- 0.1 # kpca()関数を使ったときと同じ値にする
Kg <- exp(-gamma*D^2)/Ns
A <- matrix(1/Ns,Ns,Ns)
Kg. <- Kg-A%*%Kg-Kg%*%A+A%*%Kg%*%A
```

Rのkpca()関数を使った場合と、手作業のkernel pcaの結果を並べて図示する。

```{r}
plot(as.data.frame(pcv(kernel.out)), col=as.integer(type) + 1,main="kpc()")
```
```{r}
kernel.out2 <- eigen(Kg.)
plot(as.data.frame(kernel.out2[[2]])[,1:Np], col=as.integer(type) + 1,main="manual kpc")
```

念のため、

$Kg$$K=MM^T$とが、行和、列和ともに1であるという条件を満足していることを図示しておく。

```{r}
plot(c(apply(K,1,sum),apply(K,2,sum),apply(Kg.,1,sum),apply(Kg.,2,sum)),ylim = c(-0.1,0.1),ylab="row-sum,col-sum")
```


上記で使った関数$g(x_i,x_j) = exp(-\gamma ||x_i-x_j||^2)$をガウシアン・カーネルと呼ぶ。
kpca()関数で指定したkernel="rbfdot"が、ガウシアン・カーネル選択する引数である。

これ以外の関数を使っても、うまくデータ点同士の遠近関係を与える関数を選べば、同様の手続きでkernel pcaができる。

実際、kpca()関数のヘルプ記事を見れば、8種類のカーネル関数が選べることがわかる。
```{}
rbfdot Radial Basis kernel function "Gaussian"
polydot Polynomial kernel function
vanilladot Linear kernel function
tanhdot Hyperbolic tangent kernel function
laplacedot Laplacian kernel function
besseldot Bessel kernel function
anovadot ANOVA RBF kernel function
splinedot Spline kernel
```

## カーネル・トリック

カーネルPCAではカーネル・関数を使って、固有値分解するべき正方行列を作った。

これを、以下のように考えると、トリックを使って処理を簡便化しているとみなせるので、「カーネル・トリック」と呼ぶ。

カーネル関数は、ベクトル間の遠近関係をスカラーで表す関数であるが、そのような遠近関係がうまく実現されている空間があるとする。
与えられたデータは変数の数の次元であるけれども、それより広い空間が、たまたま、変数の数の次元に押し込められていると考え、カーネル関数が定義した遠近関係が成立している高次元空間があるとする。そう考えれば、そのような空間の座標を何とかして表してやり、その広い空間でbasic-PCAを実行してやればよいことになる。

広い空間への座標変換と、そこでのbasic-PCAをするには、座標変換もしないといけないし、広い空間でのbasic-PCAもしないといけないけれど、カーネル関数を使ってできる、そんなに大きくない正方行列の固有値分解をするのと結果が変わらないのならば、カーネル関数を使って、簡単にやりましょう。

それって、カーネル関数を使ったトリックだ、ということである。

また、このことは、次のようにも言える。
与えられた変数の次元の空間では、サンプルがうまく分類できないけれど、広い空間に展開する座標変換さえあれば、固有値分解のような線形の計算でうまく分類できる。そして、結果として得られる、分類の様子を限定した軸のみで眺めると、線形では分類できないような分類になっている、というのが、カーネルPCAが非線形PCAと呼ばれる理由である。

# 欠測値に対応したPCA

データが欠けているとき、行列計算はできない。
したがって、Rの関数を使おうと手作業でやろうと、basic-PCAは実行できない。

それを補うのが、欠測値を推定しながら実施する手法である。
その手法から得られるのは、推定した欠測値とPCAの結果である。

欠測値がない場合のbasic-pcaの結果と欠測値がある場合にprobabilistic-pcaと呼ばれる欠測値対応のpcaの結果を並べてプロットする。

あまり違いはない(各軸の正負逆転などがある可能性があることに留意)。
```{r}
plot(as.data.frame(scores(basic.out)), col=as.integer(type) + 1,main="pca() basic-pca")
```
```{r}
missing.out <- pca(missing.m, method="ppca", nPcs=Np)
plot(as.data.frame(scores(missing.out)), col=as.integer(type) + 1,main="pca() probabilistic-pca")
```

推定された値を、マスクした真値と比較してプロットすれば、その推定の良さについてイメージが得られる。

```{r}
val.missing.true <- data.m[missing.cell]
val.missing.estimated <- completeObs(missing.out)[missing.cell]
plot(val.missing.true,val.missing.estimated,xlab="true val",ylab="estimated val")
abline(0,1,col=2)
```

上記の例ではmethod="ppca"と指定したが、他にもいくつもの手法が知られている。

pca()関数で選択できる手法は以下のようになる。ただし、"svd"はbasic-pcaのそれである。

```{r}
listPcaMethods()
```

# latent variable-PCA

latent variable-PCAは発展の経緯などから、さまざまな名称やモデルの複合体となっているが、ここではパッケージdmtの説明の仕方を踏襲することにする。

latent variable-PCAでは、ある相互に独立なlatent variablesが、観測variablesに異なったやり方で影響することをモデルとする。
異なったやり方での影響をモデルにするときには、2つの影響の仕方を設定することが基本なので、そのようにする。

まずは、実例をlatent variable-PCAで評価し、その出力がどうなっているかを確認し、そのやり方を理解する。

dmtパッケージに付属するデータセットmodelDataを利用して、次のようなデータセットを作る。

サンプル数Nsが51、観測変数の数が1818個の観測変数が10個と8個に二分できるとして、Ns個のサンプルのそれぞれのデータセットをXm,Ymとする。

```{r}
# library(dmt)
data(modelData)
Xm <- X #10行51列の行列
Ym <- Y[1:8,] # modelDataのYは10x51行列であるが、そのうちの8行を取り出す
dim(Xm)
dim(Ym)
```


XmとYmとは行数が同じで列数が異なる行列。

## latent variable-PCAの一つprobabilistic-PCAを実行する

probabilistic-PCAと呼ばれるlatent variable-PCAの1手法で評価してみる。

```{r}
ppca.out <- ppca(Xm,Ym)
```

基本となる出力は以下の3種類、5行列。

```{r}
# W: Wx and Wy
dim(ppca.out@W$X)
dim(ppca.out@W$Y)
# Phi: Phix and Phiy
dim(ppca.out@phi$X)
dim(ppca.out@phi$Y)
# Z: latent variables
dim(ppca.out@z)
```

これらを$W_x,W_y$,$\Phi_x,\Phi_y$,$z$とすると
$$
X_m \sim W_x z + N(0,\Phi_x)\\
Y_m \sim W_y z + N(0,\Phi_y)
$$
というモデルで推定したことになる。

式を言葉で説明すれば、『$X_m,Y_m$はどちらも$z$$W_x,W_y$で一次変換したものに、多次元正規乱数の変動が付加されたものである。ただし、多次元正規乱数には、$\Phi_x,\Phi_y$で表された、サンプル間の分散共分散関係がある』となる。

上記のモデル推定の結果を実際に出力された行列で計算しておくとわかりやすいので、以下にその計算を示す。

XmとYmとに共通するlatent variableの行列zに、zとXm,Ymそれぞれとの関係を表すWx,Wyとを掛けたものが、変動のないXm,Ymの推定値となるので、それをいったん計算し、ついで、$Phi$として推定した分散共分散行列に基づいた多次元正規乱数を付加している。

```{r}
Xm. <- ppca.out@W$X %*% ppca.out@z
Xm.. <- Xm. + t(rmvnorm(ncol(Xm),sigma=ppca.out@phi$X))
Ym. <- ppca.out@W$Y %*% ppca.out@z
Ym.. <- Ym. + t(rmvnorm(ncol(Ym),sigma=ppca.out@phi$Y))
```

このようにして算出されるXm..,Ym..が観測データセットXm,Ymと近くなるように$W_x,W_y,z,\Phi_x,\Phi_y$を線形代数計算で推定しているわけである。

## probabilistic-PCAのモデル制約

推定された$Phi_x,$Phi_y$はどんな行列かを図示してみる。
この2つの正方行列を対角成分とした行列がppca.out@phi$totalとして出力されているので、以下のようになる。
```{r}
image(ppca.out@phi$total)
```

図から読み取れることは、対角行列であって、対角成分がすべて($\Phi_x,\Phi_y$ともに)等しくなっていることである。

このように厳しく制約された$\Phi_x,\Phi_y$の意味は、Xm,Ymのそれぞれの変動は、相互に独立であり、変動の大きさも同程度である、ということである。

実は、この制約を変化させることもできる。
そして変化させたlatent variable-PCAという手法もあり、それを実行することもできる。
dmtパッケージでは、それを容易に行うために、一般化した関数fit.dependency.model()を用意し、その実行オプションとしてprobabilistic-PCAを実行することもできる。
以下のようにする。

zDimensionという引数はlatent variableの数を制約しない、という条件であるし、matched=FALSEという引数条件は、Xm,Ymの行列の形が異なることを指定しており、最後のmarginalCovariatesは$\Phi_x,\Phi_y$のことで、2つの行列が、どちらも対角行列であってすべての対角成分が同じであることを"identical isotoropic"という条件で指定している。


```{r}
ppca.out.2 <- fit.dependency.model(Xm, Ym,zDimension=NULL,matched=FALSE,marginalCovariances="identical isotropic")
```

実際、ppca()の出力とfit.dependency.model(marginalCovariances="identical isotropic")とが同じであることを確認しておく。

```{r}
ppca.out@W$X - ppca.out.2@W$X
#ppca.out@W$Y - ppca.out.2@W$Y
#ppca.out@z - ppca.out.2@z
#ppca.out@phi$X - ppca.out.2@phi$X
#ppca.out@phi$Y - ppca.out.2@phi$Y
```

## probabilistic-PCA以外のlatent variable-PCA
fit.dependency.model()関数のオプション指定によって、色々な制約の仕方でlatent variable評価ができることがわかった。

どのような制約の仕方があるのかは、関数のヘルプ記事に詳しいが、要点を述べると以下のようになる。

latent variableの数をzDimensionで指定する。この意味は特に難しいことはない。

$\Phi_x,\Phi_y$の制約を"identical isotropic", "isotropic", "diagonal" , "full"の4通りで指定して実行し、その結果の$\Phi$を図示し、イメージを掴むことにする。

$\Phi_x$の場合
```{r}
m.c <- c("identical isotropic","isotropic","diagonal","full")

```{r}
par(mfrow=c(2,2))
for(i in 1:length(m.c)){
  tmp.out <- fit.dependency.model(Xm, Ym,zDimension=NULL,matched=FALSE,marginalCovariances=m.c[i])
  image(tmp.out@phi$X,main=m.c[i])
}
par(mfcol=c(1,1))
```

$\Phi_y$の場合

```{r}
par(mfrow=c(2,2))
for(i in 1:length(m.c)){
  tmp.out <- fit.dependency.model(Xm, Ym,zDimension=NULL,matched=FALSE,marginalCovariances=m.c[i])
  image(tmp.out@phi$Y,main=m.c[i])
}
par(mfcol=c(1,1))
```

$\Phi_x,\Phi_y$を連結した場合

```{r}
par(mfrow=c(2,2))
for(i in 1:length(m.c)){
  tmp.out <- fit.dependency.model(Xm, Ym,zDimension=NULL,matched=FALSE,marginalCovariances=m.c[i])
  if(i<4){
    image(tmp.out@phi$total,main=m.c[i])
  }else{
    tmp.mat <- matrix(0,nrow(tmp.out@phi$X)+nrow(tmp.out@phi$Y),ncol(tmp.out@phi$X)+ncol(tmp.out@phi$Y))
    tmp.mat[1:nrow(tmp.out@phi$X),1:ncol(tmp.out@phi$X)] <- tmp.out@phi$X
    tmp.mat[(nrow(tmp.out@phi$X)+1):nrow(tmp.mat),(ncol(tmp.out@phi$X)+1):ncol(tmp.mat)]
    image(tmp.mat,main=m.c[i])
  }
}
par(mfcol=c(1,1))
```

# 発展

## dual-PCAとGP-LVM

basic-PCAとその応用として、カーネルPCA、欠測値対応PCA、latent variable-PCAについて、それぞれ、Rでの実装を確認した。

これらの強みは、優れた計算ライブラリがある線形代数計算を使えることである。


これらを組み合わせることで、さらに複雑なモデルに基づいてデータを解釈することもできる。

たとえば、latent variableモデルを使いながらカーネルトリックを導入して非線形な評価をする手法であるGP-LVM(Gaussian Process-Latent Variable Model)などがその例である。

その手法では、latent variable-PCAでlatent variableの行列zと、それを観測変数の行列に結びつける行列Wとについて、線形代数計算をしているが、その線形代数処理の解釈を工夫することで、$K=MM^T$のような内積型の行列を取扱い対象とし、その固有値分解へと問題を変化させることができる。
そのような変形をlatent variable-PCAの一つであるprobabilistic-PCAに施したのがprobabilistic dual-PCAである。

さらに、probabilistic dual-PCAの内積行列にカーネル・関数を結びつけることで、probabilisitc dual-PCAの(そして結果的にlatent variable-PCAの一つであるprobabilistic-PCAの)非線形版を構成することが可能である。

これがGP-LVMである。

ただし、GP-LVMに至るにあたって、dual化での変形とカーネル関数の導入との2ステップを踏んだ結果、計算は単純な固有値分解に帰着されないので、固有値分解を絡めた最適化問題へと変化している点で、basic-PCAからの単純な発展とは言い切れない手法ではある。