数学

[数学]N次元回転行列の成分表示

TensorFlow実験をしようと思ってN次元回転行列を返す関数をNumpyで探していたのですがどうも実装されていないので自分で実装することにしました.N次元回転行列を検索するとどうしても成分表示している文献を見つけられなかったので導出しました.

前提

「2つのN次元ベクトルが張る平面内の回転」とは,N次元ベクトルをその平面に射影して二次元の回転をしたあと,射影された時に潰した成分を加えて,再びN次元ベクトルを作る操作をいうこととします.

表記

  • {\bf u,v}…回転させたい平面を与えるベクトル
  • {\bf x}…一般のN次元ベクトル
  • R…2つのN次元ベクトルが張る平面内の回転を与える行列
  • r_{ij}…Rの(i,j)成分.成分表示の基底は{\bf x}と同じとする.
  • \theta…回転角
  • {\bf e_u}=\frac{{\bf u}}{|{\bf u}|},{\bf e_v}=\frac{{\bf v-(v\cdot e_u)e_u}}{|{\bf v-(v\cdot e_u)e_u}|}
  • \{{\bf e_i}\}…N次元の正規直交基底
  • u_i, v_i{\bf e_u,e_v}の第i成分(i=1~N)

結論

Rの成分の表式は以下のようになります.

r_{ij}=u_i [(\cos \theta -1)u_j-\sin \theta v_j]+v_i[\sin \theta u_j+(\cos\theta-1)v_j]+\delta_{ij}

導出

無題 1

図が汚くてすみません笑.これはベクトルの配置を3次元でイメージするための図です.まずu,vが張る平面へのxの正射影を計算しましょう.

{\bf p}={\bf (x\cdot e_u)e_u+(x\cdot e_v)e_v}

平面内の成分表示されたベクトルの回転行列はよく知られているように\begin{pmatrix}\cos\theta & -\sin\theta \\\ \sin\theta & \cos\theta \end{pmatrix}と表せます.したがって{\bf p}を平面内で回転させると

{\bf [\cos\theta(x\cdot e_u)-\sin\theta(x\cdot e_v)]e_u+[\sin\theta(x\cdot e_u)+\cos\theta(x\cdot e_v)]e_v}

に移ります.これに垂線ベクトル{\bf h}={\bf x-p}を加えると回転後の{\bf x}が求まるので

R{\bf x}={\bf [(\cos\theta-\textrm{1})(x\cdot e_u)-\sin\theta(x\cdot e_v)]e_u+[\sin\theta(x\cdot e_u)+(\cos\theta-\textrm{1})(x\cdot e_v)]e_v+x}

なります.さて,{\bf x,e_u,e_v}をN次元の正規直交基底\{{\bf e_i}\}を用いて

{\bf x}=\sum_i x_i{\bf e_i}, {\bf e_u}=\sum_i u_i{\bf e_i}, {\bf e_v}=\sum_i v_i{\bf e_i}

と表しましょう(u_i,v_iは紛らわしいですが{\bf e_u,e_v}の成分です).先ほどの式の右辺にこれらを代入し,えいやっと変形して{\bf e_i},{\bf x_i}についてまとめると次の式が得られます.

\sum_i\left[\sum_j\left\{u_i[(\cos\theta-1)u_j-\sin\theta v_j]+v_i[\sin\theta v_j+(\cos\theta-1)v_j]+\delta_{ij}\right\}x_j\right]{\bf e_i}

一方左辺は次のようになります.

\sum_i(\sum_j r_{ij}x_j){\bf e_i}

左辺と右辺を比較することで

r_{ij}=u_i [(\cos \theta -1)u_j-\sin \theta v_j]+v_i[\sin \theta u_j+(\cos\theta-1)v_j]+\delta_{ij}

を得ます.

追記:

{\bf u,v}を与えて{\bf u}{\bf v}に重ねる関数の実装です.回転角を適当に変えたい場合は関数中のthetaを関数の引数にするなどすると良いです.


#fromVecをrotatedVecに重ねる回転行列を返す
def RotationMatrix(fromVec,rotatedVec):
    if len(fromVec)!=len(rotatedVec):
        raise ArithmeticError('Incompatible vector length')
    eu=fromVec/np.linalg.norm(fromVec)
    ev=rotatedVec-np.dot(rotatedVec,eu)*eu
    ev/=np.linalg.norm(ev)
    theta=np.arccos(np.dot(eu,rotatedVec)/np.linalg.norm(rotatedVec))
    ret=np.zeros([len(fromVec),len(fromVec)])
    for i in range(len(fromVec)):
        for j in range(len(rotatedVec)):
            ret[i,j]=eu[i]*((np.cos(theta)-1)*eu[j]-np.sin(theta)*ev[j])+ev[i]*(np.sin(theta)*eu[j]+(np.cos(theta)-1)*ev[j])
    ret[i,i]+=1
    return ret