R de Isomap

RでIsomapを書いてみた。
ただそれだけ。

まだあんまりRのことは分かってないんだけど、for文を使うと明らかに実効速度的に不利であることは判明した。
applyとかでうまく回避するんだろうけど、C言語育ちの私にとっては「行列の全ての要素に何らかの処理を行う」ってなるとすぐにfor文が頭に浮かんでしまう。
というわけで、僕の書いたIsomapには二重ループがやたらと登場してきて実行速度的に速度的に非常にだめだめです。
どうしたものか。

まともに固有値・固有ベクトルを求めてソートをかけるのがめんどくさかったので、主成分分析の関数を代用してみたんだけどこれでいいのだろうか？
まあ、前にPythonで書いたやつと結果が大きく違わないからいいんだろうけど...
あと、eigen(A)とprincomp(A)とprcomp(A)で固有値が違う気がするのは俺だけ？

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

# 近傍の数
k <- 5

# 各点同士のユークリッド距離を計算
# データ点×データ点の行列を作成
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近傍のみ値を持たせる(k近傍グラフ作成)
for (i in 1:length(swiss[,1])){
  print(i)
  odistmat <- order(distmat[i,])
  for (j in 1:length(swiss[,1])){
    if (odistmat[j] > k) distmat[i,j] <- Inf
  }
}
for (i in 1:length(distmat[,1])){
  for (j in i:length(distmat[,1])){
    distmat[j,i] <- distmat[i,j]
  }
}

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

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

D <- distmat^2

# Centering Matrix作成
H <- diag(1, length(distmat[1,])) - (1/length(distmat[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)