2017-01-08

回転群 SO3のランダムな要素生成

回転群回転球面

こちらにある通り、４次元空間にあるS3球面上の点に相当する四元数 $(\sqrt{1-u_1}\cos{2\pi u_2},\sqrt{1-u_1}\sin{2\pi u_2},\sqrt{u_1}\cos{2\pi u_3},\sqrt{u_1}\sin{2\pi u_3})$ が表す回転がそれである
3つの一様乱数で定まっている。 $u_2,u_3$ はそれぞれS1 (円)を一様に定める変数。 $u_1$ は $(\sqrt{1-u_1})^2 + (\sqrt{u_1})^2=1$ であるから、 $\cos{\theta}$ (と $\sin{\theta}$ と)の分布

my.random.so3 <-function(n){
	# require(onion)
	u <- matrix(runif(3*n),ncol=3)
	h <- sqrt(1-u[,1]) * sin(2*pi*u[,2]) + Hi * sqrt(1-u[,1]) * cos(2*pi*u[,2]) + Hj * sqrt(u[,1]) * sin(2*pi*u[,3]) + Hk * sqrt(u[,1])*cos(2*pi*u[,3])
	return(h)
}

n <- 10^4
hs <- my.random.so3(n)
# 四元数の回転方向は３次元空間で一様
plot3d(i(hs),j(hs),k(hs))

x <- rnorm(3)
x.h <- Hi*x[1]+Hj*x[2]+Hk*x[3]

rot.x <- Conj(hs) * x.h * hs
rot.x. <- cbind(i(rot.x),j(rot.x),k(rot.x))
# ある点が、この一様要素で回転すると、
# 回転して移動した先は一様
plot3d(rot.x.)

これをよく見ると、回転軸は３次元単位ベクトルをランダムに発生させて決まるので、それは、３次元標準正規乱数で発生させることにすれば、回転角だけを $\sqrt{1-n} * \sin{2\pi n}$ で発生させればよいことが解る

my.random.so3.2 <- function(n){
	R <- matrix(rnorm(n*3),ncol=3)
	R <- R/sqrt(apply(R^2,1,sum))
	thetas <- acos(sqrt(1-runif(n)) * sin(2*pi*runif(n))) * 2
	h <- cos(thetas/2) + sin(thetas/2) * (Hi*R[,1]+Hj*R[,2]+Hk*R[,3])
	return(h)
}
-行列表現は次のように作れる
>|r|
my.random.so3.matrix <- function(n){
	h <- my.random.so3(n)
	ret <- list()
	e1 <- 1*Hi
	e2 <- 1*Hj
	e3 <- 1*Hk
	for(i in 1:n){
		this.h <- h[i]
		tmp1 <- Conj(this.h) * e1 * this.h
		tmp2 <- Conj(this.h) * e2 * this.h
		tmp3 <- Conj(this.h) * e3 * this.h
		ret[[i]] <- cbind(c(i(tmp1),j(tmp1),k(tmp1)),c(i(tmp2),j(tmp2),k(tmp2)),c(i(tmp3),j(tmp3),k(tmp3)))
	}
	return(ret)
}

QR分解を使って、ランダムな回転行列を作ることもある

library(GPArotation)
Random.Start(3)

そのソースは、以下のように、要素が正規乱数である正方行列のQR分解によっている

> Random.Start
function (k) 
{
    qr.Q(qr(matrix(rnorm(k * k), k)))
}
<environment: namespace:GPArotation>

ただし、これによって回すと…

# qr分解を使ったランダムな回転行列 RGProtation::Random.Start()だと
# ある点の移動先が一様にならない
rot.x.. <- rot.x.
for(i in 1:n){
	RS <- Random.Start(3)
	rot.x..[i,] <- RS %*% x
}
plot3d(rot.x..)

移動先は、球面に一様分布しない
四元数の方をもう少し調べる

plot(as.data.frame(cbind(i(hs),j(hs),k(hs))))
# 四元数の四成分は、４次元空間の球面(原点からの距離が１の点の集合)になっている
# しかも一様の模様
plot(as.data.frame(cbind(Re(hs),i(hs),j(hs),k(hs))))
Mod(hs)

では、４次元空間の単位球面上の一様乱点を発生させて、それに対応する四元数を使って回転すればよいだけでは…

n <- 10^4
RR <- matrix(rnorm(n*4),ncol=4)
RR <- RR/sqrt(apply(RR^2,1,sum))
hs <- RR[,1] + Hi*RR[,2] + Hj*RR[,3] + Hk*RR[,4]
Mod(hh)
x <- rnorm(3)
x.h <- Hi*x[1]+Hj*x[2]+Hk*x[3]
rot.x <- Conj(hs) * x.h * hs
rot.x. <- cbind(i(rot.x),j(rot.x),k(rot.x))
plot3d(rot.x.)