tSNE

  • こちらで多次元のcloud dataを可視化するcytoSPADEのことを書いた
  • 今日は多様体学習で2次元埋め込みをするtSNEを
  • 資料はこちら
  • 基本処理
    • データ点が多くない場合は全データ点を用いる
    • 高次元空間でのペアワイズな点間距離を、観測データから定める
    • 2次元(もしくはその他の下げた次元)の対応点の座標について、ペアワイズな点間距離を同様に、低次元空間での観測点から定める(ここでKL divergenceを使う:この点が別のこの点として観察される確率は?のようにして、ある点がどの点と遠いか近いかを相対的に(尤度的に)評価する)
    • 高次元でのペアワイズ点間距離分布と低次元でのそれとの異同をペアワイズなKL divergenceをまとめたもの(和だったり、joint distributionのそれだったりとする)で測り、それが最小になるように低次元座標を求める(勾配法)
    • 高次元・低次元の2点間距離は、それ以外の点との距離との相対的なものとして定める。そのとき、正規分布(距離の2乗の負の指数に比例)を考慮するのが、基本形
    • 本例では、KLdivergenceの非対称性を解消するように、相対的な距離の定義を変更し、また、正規分布ではなくてt分布を用いることで当てはまりが変わらず(もしくはよりよく)、かつ、計算が単純になるようにした
  • データ点が多い場合の工夫
    • 点の数が多いと大変
    • 一部の点をランダムに抜き出す(ダウンサンプリングとかはしない)
    • その上で、選ばれた点のペアワイズ距離は、全ての点のneighbor graph上のランダムウォーク(エッジの選び方はエッジの距離の2乗の負数の指数で採択確率の重みづけをする)によって、到達する割合によって定める
    • そのランダムウォークを実際にモンテカルロでやってもよいし、疎な行列に関して解析的に解いてもよいらしい
    • 点を選んで処理した後、残った点はどうするの?
    • ランダムウォークが重くて話にならないソースだけれど、意味あいとしては以下のような感じ
library(vegan)
library(GPArotation)
library(rgl)
library(igraph)
library(MCMCpack)
# データの座標によらず、視野内に収めるための補正
vision.scale <- function(x,scl){
	(x-min(x))/(max(x)-min(x)) * scl + (1-scl)/2
}
# 観察対象をきちんと視野内に収めるために、周辺部にはシグナルがないように調整する
vision.standard <- function(x,scl = 0.9){
	rg <- apply(x,2,range)
	tmp <- apply(x,2,vision.scale, scl)
	return(list(st.x = tmp,rg=rg))
}
# 視野全体のピクセル数を指定して、該当する超立方体の番地にする

round.coords <- function(x,n.cell=50){
	d <- length(x[1,])
	matrix(floor(x * n.cell),ncol=d)
}


# 適当に線状の点の集まりを作る
d <- 3

n.steps <- 10
n.pt.range <- 50:100
NN <- 100
X <- matrix(0,0,d)

for(i in 1:n.steps){
	tmp.n <- sample(n.pt.range,1)
	tmp.X <- matrix(0,tmp.n,d)
	tmp.L <- runif(1,5,10)
	tmp.V <- seq(from=0,to=1,length=tmp.n*NN)
	tmp.prob <- (tmp.V-mean(tmp.V))^2
	tmp.prob <- tmp.prob/sum(tmp.prob)
	tmp.S <- tmp.V[sort(sample(1:length(tmp.V),tmp.n,replace=TRUE,prob=tmp.prob))]*tmp.L
	tmp.X[,1] <- tmp.S
	tmp.R <- Random.Start(d)
	tmp.X <- tmp.X %*% tmp.R

	if(i > 1){
		tmp.st <- sample(1:length(X[,1]),1)
		#tmp.st.2 <- sample(1:tmp.n,1)
		tmp.st.2 <- 1
		st.pt <- X[tmp.st,]
		connect.pt <- tmp.X[tmp.st.2,]
		tmp.X <- t(t(tmp.X) + st.pt - connect.pt)
	}
	
	X <- rbind(X,tmp.X)
}


X <- X + rnorm(length(X),0,0.3)
plot3d(X)
X.curve <- X
for(i in 1:d){
	X.curve[,i] <- sign(X.curve[,i]) * abs(X.curve[,i]^2)
}
X <- X.curve
plot3d(X)

# Distance matrix
D <- as.matrix(dist(X,method="manhattan"))
diag(D) <- max(D)+1
neighb <- 5

D.neighbor <- matrix(0,length(D[,1]),neighb)
D.neighbor.dist <- matrix(0,length(D[,1]),neighb)
for(i in 1:length(D.neighbor[,1])){
	tmp <- D[i,]
	tmp2 <- sort(tmp)[neighb]
	tmp.neighb <- which(tmp<=tmp2)
	D.neighbor[i,] <- tmp.neighb
	D.neighbor.dist[i,] <- tmp[tmp.neighb]
}
D.neighbor.prob <- exp(-D.neighbor.dist^2)
D.neighbor.prob.sum <- apply(D.neighbor.prob,1,sum)
D.neighbor.prob <- D.neighbor.prob/D.neighbor.prob.sum

# No. points
n <- length(X[,1])
print(n) # あまり多くなくやっておこう
# n.pick個の点について最近点距離を求める
n.pick <- 100 # 2000ではなく少なめでやってみる
picked <- sample(1:n,n.pick)

n.iter.walk <- 10
end.walk <- matrix(0,n.pick,n.pick)
for(i in 1:n.pick){
	for(j in 1:n.iter.walk){
		loop <- TRUE
		wlk <- c(picked[i])
		while(loop){
			tail <- wlk[length(wlk)]
			tmp.neighbor <- D.neighbor[tail,]
			tmp.prob <- D.neighbor.prob[tail,]
			non.start <- which(tmp.neighbor!=picked[i])
			tmp.neighbor <- tmp.neighbor[non.start]
			tmp.prob <- tmp.prob[non.start]
			tmp.prob <- tmp.prob/sum(tmp.prob)
			tmp <- sample(tmp.neighbor,1,prob=tmp.prob)
			#if(tmp==wlk[1]){
				#wlk <- c(picked[i])
			#}else{
				wlk <- c(wlk,tmp)
				if(prod(picked-wlk[length(wlk)])==0){
					loop <- FALSE
				}
			#}
			
		}
		tmp <- which(picked==wlk[length(wlk)])
		end.walk[i,tmp] <- end.walk[i,tmp] + 1
	}
}