2013-11-16

移動とマッチングに関するメモ

R 移動マッチング

空間に点集合が２つ( $V_1,V_2$ )ある。一つ目の集合( $V_1$ )から二つ目の集合( $V_2$ )への移動であると考えて、どの点がどの点に移動したとみなすことにするかを考える
同数の点に関するマッチングの決め方であるとみなせる
マッチングにはハンガリアン・アルゴリズムが使える場合がある
- その条件は、n点からn点へのマッチングであるとして、nxn行列の非負正方行列( $K$ )が定義できて、第１集合( $V_1$ )の第i点( $x_{1,i}$ を第２集合( $V_2$ )の第j点( $x_{2,j}$ )にマッチングさせたときにスコア $k_{i,j}$ が与えられて、 $(1,2,...,n) \to (s_1,s_2,...,s_n)$ なるマッチングにおいて $\sum_{i=1}^n k_{i,s_i}$ の最大化もしくは最小化問題にできる場合である
の決め方は一通りではない
- $k_{i,j}=\sqrt{\sum_{u=1}(x_{1,i}(u)-x_{2,j}(u))^2}$ のように、２点のユークリッド距離にしてもよいし
- $k_{i,j}=\sum_{u=1}|x_{1,i}(u)-x_{2,j}(u)|$ のように、２点のマンハッタン距離にしてもよい
- (内積)のようなものを考えてもよい。この場合は、式変形に注意すれば[tex:k_{i,j} = \sum_{u=1} *1 - (\frac{\sum_{i=1}^n x_{1,i}(u)}{n} - \frac{\sum_{i=1}^n x_{2,i}(u)}{n}))^2]と大小関係が保たれていることがわかる
  - この例は、すべての２点ペアを結ぶベクトルの多次元分散の大小に相当している
- さらに、２点ペアのベクトルが似通っていることを求めるのであれば、すべてのペアのベクトルが似ていることを求めるかわりに、近い関係にあるそれにのみ求めるという方法も考えられ(それは非線形なアプローチなのだが)るが、それをnxn行列に $i,j$ のペアとして実現することはできない(もしくは困難)なので、それはハンガリアン・アルゴリズムでの処理には合わないもののようだ

# ノード対応づけをするために試行錯誤した関数
# マッチングでできたペアのユークリッド距離
library(clue)
my.move <- function(x,y,w.x=NULL,w.y=NULL){
	if(!is.matrix(x)){
		x <- matrix(x,ncol=1)
	}
	if(!is.matrix(y)){
		y < matrix(y,ncol=1)
	}
	if(is.null(w.x)){
		w.x <- rep(1,length(x[,1]))
	}
	if(is.null(w.y)){
		w.y <- rep(1,length(y[,1]))
	}
	D <- x %*% t(y)
	D <- D-min(D)
	xy <- rbind(x,y)
	m.xy <- apply(xy,2,mean)
	x.d <- t(t(x)-m.xy)
	y.d <- t(t(y)-m.xy)
	W <- outer(w.x,w.y,"*")
	D <- D * W
	
	as.vector(solve_LSAP(D,maximum=TRUE))
}
# ペアのベクトルの分散
my.move.D <- function(x,y,w.x=NULL,w.y=NULL){
	if(!is.matrix(x)){
		x <- matrix(x,ncol=1)
	}
	if(!is.matrix(y)){
		y < matrix(y,ncol=1)
	}
	if(is.null(w.x)){
		w.x <- rep(1,length(x[,1]))
	}
	if(is.null(w.y)){
		w.y <- rep(1,length(y[,1]))
	}
	D <- (as.matrix(dist(rbind(x,y))))[1:length(x[,1]),(1+length(x[,1])):(2*length(x[,1]))]
	#D <- D-min(D)
	xy <- rbind(x,y)
	m.xy <- apply(xy,2,mean)
	W <- outer(w.x,w.y,"*")
	D <- D * W
	
	as.vector(solve_LSAP(D,maximum=FALSE))
}
# 相互のベクトルの近さで工夫する？？
my.move.L <- function(x,y,w.x=NULL,w.y=NULL){
	if(!is.matrix(x)){
		x <- matrix(x,ncol=1)
	}
	if(!is.matrix(y)){
		y < matrix(y,ncol=1)
	}
	if(is.null(w.x)){
		w.x <- rep(1,length(x[,1]))
	}
	if(is.null(w.y)){
		w.y <- rep(1,length(y[,1]))
	}
	D <- x %*% t(y)
	D <- D-min(D)
	xy <- rbind(x,y)
	m.xy <- apply(xy,2,mean)
	x.d <- t(t(x)-m.xy)
	y.d <- t(t(y)-m.xy)
	L.x <- sqrt(apply(x.d^2,1,sum))
	L.y <- sqrt(apply(y.d^2,1,sum))
	LL <- outer(L.x,L.x,"*")
	L.q <- quantile(LL,0.05)
	LL.q <- sign(LL < L.q)
	W <- outer(w.x,w.y,"*")
	D <- D * W * LL.q
	
	as.vector(solve_LSAP(D,maximum=TRUE))
}
v1 <- runif(2)
v2 <- runif(2)
t <- seq(from=0,to=1,length=100)
x <- matrix(0,length(t),2)
for(i in 1:length(t)){
	x[i,] <- t[i] * v1 + (1-t[i]) * v2
}
x[,1] <- x[,1]-mean(x[,1])
x[,2] <- x[,2]-mean(x[,2])
theta <- pi/3
M <- matrix(c(cos(theta),sin(theta),-sin(theta),cos(theta)),2,2)
x <- rbind(x,x%*%M)
theta <- pi/6
M <- matrix(c(cos(theta),sin(theta),-sin(theta),cos(theta)),2,2)

y <- x %*% M
y <- y + 0.3
y <- (abs(y))^1.3
out <- my.move(x,y)
out.D <- my.move.D(x,y)
out.L <- my.move.L(x,y)
par(mfcol=c(2,3))
plot(rbind(x,y),pch=20,col=rep(1:2,each=length(x[,1])),main="内積")
for(i in 1:length(out)){
	segments(x[i,1],x[i,2],y[out[i],1],y[out[i],2],col=3)
}
plot(out,main="内積")
plot(rbind(x,y),pch=20,col=rep(1:2,each=length(x[,1])),main="ユークリッド距離")
for(i in 1:length(out)){
	segments(x[i,1],x[i,2],y[out.D[i],1],y[out.D[i],2],col=3)
}
plot(out.D,main="ユークリッド距離")
plot(rbind(x,y),pch=20,col=rep(1:2,each=length(x[,1])))
for(i in 1:length(out)){
	segments(x[i,1],x[i,2],y[out.L[i],1],y[out.L[i],2],col=3,main="近傍重視・内積")
}
plot(out.L,main="近傍重視・内積")
par(mfcol=c(1,1))

*1:x_{1,i}(u)-x_{2,j}(u