高次元での回転と疎行列

R 形・動き Rpackage Rを使いこなす 数学

MikuHatsune2016-07-11

3次元物体を回転させたい。
2次元であれば回転行列は回転角\thetaを用いてR=\begin{pmatrix}\cos\theta&\sin\theta\\-\sin\theta&\cos\theta\end{pmatrix} と書ける。
高次元での回転ならば、正方対角行列(対角成分は1)Rについて、(i, i), (i, j), (j, i), (j, j) 成分をひらすら上のRで置換していけば、ある軸での回転を表す。
特に3次元では、x軸まわりに\alpha、y軸まわりに\beta、z軸まわりに\gamma 回転することを考えると、v=(x,y,z)を回転させて得られる座標V=(X,Y,Z)は
^{t}V=\begin{pmatrix}\cos\gamma&\sin\gamma&0\\-\sin\gamma&\cos\gamma&0\\0&0&1\end{pmatrix}\begin{pmatrix}\cos\beta&0&-\sin\beta\\0&1&0\\\sin\beta&0&\cos\beta\end{pmatrix}\begin{pmatrix}1&0&0\\0&\cos\alpha&\sin\alpha\\0&-\sin\alpha&\cos\alpha\end{pmatrix}{^{t}v}
と表される。


　
これを疎行列でやる。疎行列自体の説明は飛ばして、実際に使ってみる。
基本的な疎行列の形は、例えば対角成分だけ何か値があって他は0というような、
S=\begin{pmatrix}s_1&0&0&0\\0&s_2&0&0\\0&0&s_3&0\\0&0&0&s_4\end{pmatrix}
というものがまず思いつくやつである。ここで、s_iも行列だと思えば、例えばこれがxyz 座標に相当するとして、s_i=\begin{pmatrix}x_i&0&0\\0&y_i&0\\0&0&z_i\end{pmatrix} と考えれば、
S=\begin{pmatrix}x_1&0&0&&&&\\0&y_1&0&\dots&&0&\\0&0&z_1&&&&\\&\vdots&&\ddots&&\vdots&\\&&&&x_n&0&0\\&0&&\dots&0&y_n&0\\&&&&0&0&z_n\end{pmatrix}
というような対角ブロック行列ができる。
ブロックに相当する回転行列を持った、これまた疎行列を用意すれば、すべてが行列演算でできる。
しかし、xyz座標をベタに呼び出して計算したほうが速いっぽい…メモリサイズ的には疎行列のほうがよさそうなんだが。

rotate_polar <- function(xyz, alpha=pi/3, beta=pi/3, gamma=pi/3, sparse=FALSE){
  # angle alpha around x axis, beta y axis, and gamma z axis
  X <- matrix(c(cos(gamma), -sin(gamma), 0, sin(gamma), cos(gamma), 0, 0, 0, 1), nr=3)
  Y <- matrix(c(cos(beta), 0, sin(beta), 0, 1, 0, -sin(beta), 0, cos(beta)),     nr=3)
  Z <- matrix(c(1, 0, 0, 0, cos(alpha), -sin(alpha), 0, sin(alpha), cos(alpha)), nr=3)
  if (!sparse){
    res <- t(mapply(function(w) X %*% Y %*% Z %*% xyz[w,], seq(nrow(xyz))))
  } else {
    ij <- seq(nrow(g1)*3)
    M <- sparseMatrix(i=ij, j=ij, x=c(t(g1)))
    zeroX <- sparseMatrix(i=rep(1:(nrow(g1)*3), each=3), j=rep(1:(nrow(g1)*3), 3), x=as.numeric(X))
    zeroY <- sparseMatrix(i=rep(1:(nrow(g1)*3), each=3), j=rep(1:(nrow(g1)*3), 3), x=as.numeric(Y))
    zeroZ <- sparseMatrix(i=rep(1:(nrow(g1)*3), each=3), j=rep(1:(nrow(g1)*3), 3), x=as.numeric(Z))
    res <- matrix(rowSums(zeroX %*% zeroY %*% zeroZ %*% M), nc=3, byrow=TRUE)
  }
  return(res)
}

xyz <- matrix(rnorm(1000*3), nc=3)
# 疎行列にしてみる
system.time(rotate_polar(xyz, sparse=TRUE))

   ユーザ   システム       経過  
     0.257      0.000      0.258 

# 疎行列を使わずベタにfor loop で呼び出す
system.time(rotate_polar(xyz, sparse=FALSE))

   ユーザ   システム       経過  
     0.001      0.001      0.003