R de Isomap再掲

研究

前回のRで書いたIsomapはk近傍グラフの作り方が決定的に間違っていました。
なんかブクマつけてくれている人もいるので修正版を・・・

# データ取得
swiss <- read.table('swiss.data')

# 近傍の数
k <- 8

# 各点同士のユークリッド距離を計算
# データ点×データ点の行列を作成
distmat <- diag(0,length(swiss[,1]))

# 各点同士のユークリッド距離を求める関数
distance <- function(x,y) {
  tmp <- 0
  for (i in 1:length(x)) {
    tmp <- tmp + (x[i] - y[i])^2
  }
  sqrt(tmp)
}


# 関数の適用
for (i in 1:length(swiss[,1])){
  print(i)
  for (j in i:length(swiss[,1])){
    distmat[i,j] <- distance(swiss[i,],swiss[j,])
  }
}
for (i in 1:length(swiss[,1])){
  for (j in i:length(swiss[,1])) {
    distmat[j,i] <- distmat[i,j]
  }
}

print(distmat)

# k近傍を求める
ddistmat <- diag(0,length(swiss[,1]))
for (i in 1:length(swiss[,1])){
  print(i)
  odistmat <- order(distmat[i,])
  for (j in 1:length(swiss[,1])){
    for (l in 1:k  ) {
      if (distmat[i,j] == distmat[i,odistmat[l]]) {
        ddistmat[i,j] <- distmat[i,odistmat[l]]
        break
      }
      else {
        ddistmat[i,j] <- Inf
      }
    }
  }
}
for (i in 1:length(distmat[,1])) {
  for (j in i:length(distmat[,1])) {
    ddistmat[j,i] <- ddistmat[i,j]
  }
}

# Floyd Warshall法で近傍グラフ上の最短経路探索
for (k in 1:length(ddistmat[1,])){
  print(k)
  for (i in 1:length(ddistmat[1,])){
    for (j in i:length(ddistmat[1,])){
      if (ddistmat[i,j] > ddistmat[i,k] + ddistmat[k,j]) {
        ddistmat[i,j] <- ddistmat[i,k] + ddistmat[k,j]
      }
    }
  }
}
for (i in 1:length(ddistmat[,1])){
  for (j in i:length(ddistmat[,1])){
    ddistmat[j,i] <- ddistmat[i,j]
  }
}

for (i in 1:length(ddistmat[1,])){
  for (j in 1:length(ddistmat[1,])){
    if (ddistmat[i,j] == Inf) ddistmat[i,j] <- 0
  }
}

D <- ddistmat^2

# Centering Matrix作成
H <- diag(1, length(ddistmat[1,])) - (1/length(ddistmat[1,]))

# 測地線距離行列にdouble centeringを行う
gmat <- -0.5 * (H %*% D %*% H)

# 測地線距離行列の固有値分解
gmat.pc <- prcomp(gmat)

# 固有値の平方根
eigenvalues <- gmat.pc$sdev

# 各次元の値の計算
dimension <- 2
projection <- matrix(0,dimension,length(swiss[,1]))
for (i in 1:dimension) {
  projection[i,] <- eigenvalues[i]*gmat.pc$rotation[,i]
}

# データの出力
write(projection, 'iso_swiss.data', ncolumns=2)

これを使うと、まきまきのデータ（通称スイスロール）が


こんな感じで二次元に写像されます。


色を変えて見やすくするのが面倒くさかったので、こちらを見るとわかりやすいですよ。
Isomap Homepage
http://d.hatena.ne.jp/audioswitch/searchdiary?word=Isomap&.submit=%B8%A1%BA%F7&type=detail


まきまきしているので、線形写像である主成分分析で写像を行ってもうまくいかないことが予想されます。
カーネルを使うと・・・どうなんでしょうね。少なくともRBFカーネルなんかではうまくいきそうにないですし（同心円なので）、他のカーネルでもうまくパラメータを決めてやらないと難しいのでは？


というわけで、多様体学習は結構おもしろくて、実用的だと思われます。
自分のまわりではあまり研究している人はいないのですが、けっこう盛んに研究されているようなことを聞いたような気がします。
非線形なデータをどう扱うかというのは情報科学の一つのトピックではありますからね。