がー。

まじかよ。回転行列の生成元って、交換子はパウリ行列の交換関係と一致するけれど、SxSy=Szみたいなのは満たさないじゃん(一方でパウリ行列は四元数の基底として正しく使用できる)。積が閉じない以上今まで書いてきた話は当然にして全部間違っていた。SU(2)とSO(3)の局所同型の話って、完全にLie群、Lie環の話であって、一方で四元数の話は、Lie群、Lie環の話とは全然向いている方向が違うんだなあ。

クオータニオン⇒回転の対応は、じゃあそんなに自然なものじゃなくて、二段構えになっているのかな？つまり、四元数が直接対応するのはIとパウリ行列を基底とする実数係数線形空間なわけで、SU(2)がその中に入っているのは、リー環がまさにパウリ行列で張られて、かつIを加えると(交換子のみならず通常の)積が閉じている、という議論から言えてこれは比較的自然なストーリーだ。一方でさらにトレースレスのエルミート行列と三次元ベクトルを
[z/2, x/2+iy/2]
[x/2-iy/2,-z/2]
みたいな感じで対応させるとこれもこの線形空間に入っていて、しかもSU(2)は``三次元ベクトルを回転''させる、という話が介在してようやく三次元ベクトルの回転に四元数が使えるわけだ。多分そういう二段構造になっている。

ラッキーな事に、対角化でなく加法定理を調べる事を選択したから、その部分はパウリ行列の線形結合の指数関数の話にそのまま使える。ただし、パウリ行列はブロッホ球上の``点''を角速度2で回転させるようなジェネレータだから回転時間が回転角の半分で済むので、θがθ/2の形で入るわけだ。

四元数がベクトルを回転させるのって、もっと自然な話だと思ってたんだけど、偶々としか言いようがないぞ、これ。高次元に多分拡張できないし。

まあいいや、なんとなくクオータニオンの使い方は分かってきたし。でも、剛体の回転シミュレーションはまだ敷居が高いかも。まあどうせ遊びだ。気長に気長に。

---

追記：パウリ行列でなく、それにiを乗じた物、つまり

J=[i, 0]
  [0,-i]

K=[ 0,1]
  [-1,0]

L=[0,i]
  [i,0]

に単位行列Iを加えた物が四元数の単位要素と対応する。考えてみればこれらをSU(2)のジェネレータとして使う場合、iを乗じてから指数の肩に乗っけるのだから、以前の加法定理の議論を流用するには、当然そうあるべきだった。これでようやく、これら単位要素は二乗すると-Iになるという虚数単位と同様の振る舞いをするようになり、I, -I, J, -J, K, -K, L, -Lの間で積が閉じる。さらに行列としての共役転置も閉じ、これは四元数の共役と対応する。後でSU(2)と四元数の対応を書くが、ユニタリ行列の逆は共役転置なので、SU(2)に対応させた四元数の逆数は共役により簡単に得られる事になる。

3次元ベクトルと行列の対応は

[ix/2    , y/2+z/2]=xJ+yK+zL
[-y/2+z/2, -ix/2  ]

となり、トレースレスの反エルミート行列となる。こうして3次元ベクトルは純四元数(I成分がない)に対応するだけでなく、見かけの上で数ベクトルとしての三次元ベクトルをそのまま四元数と同一視できる。1/2の係数は、共役転置との積が、つまり四元数としてのノルムの二乗が、ベクトルのノルムの二乗と一致するようにするためだ。

一方でSU(2)も四元数と、|(nx,ny,nz)=1|の時

exp(t(nxJ+nyK+nzL))=cos(t)I+sin(t)(nxJ+nyK+nzL)

で対応する。SU(2)が四元数に対応する、つまりI,J,K,Lが張る実係数線形空間内にある事は既に書いたI, -I, J, -J, K, -K, L, -Lの間の(交換子のみならず)積が閉じる事から直ちに言えるが、この具体的な表式からさらに次の事実を見出す事ができる。すなわち、SU(2)の要素は四元数としてノルム1であり、しかもノルム1の四元数は全てこの形に書ける。つまり、SU(2)は、四元数体の中にノルム1である物からなる部分集合として埋め込まれる。既に注意したように、SU(2)の場合、四元数としての共役と、四元数としての逆数と、行列としての共役転置と、行列としての逆行列が全て対応する。

これら二つの全く異なる意味での対応関係と、上で示した形で表したSU(2)がxJ+yK+zLの形の``三次元ベクトル''をtの2倍の角度だけ(nx,ny,nz)を軸として``回転''させるという話まで考えて、結局三つの半ば独立した話を全て合わせてようやくクオータニオンが回転変換に使用できるというわけだ。

元から全く自然なストーリではない非自明な結果であり、それ故に使い方だけを説明される場合が多いわけだが、なんとなく自然なストーリーが背後にありそうに見えてしまって、勝手な哲学を見出そうとして、それで混乱するわけだ。q v conj(q)でベクトルvが回転します、なんてあっさり書いてあると、まるで、ベクトルを四元数が``直接''回転するように見えるのでどうにか自然な説明を見出そうとしてしまいがちになるが、実際にはこのように、何段階もの対応関係を経て、ようやく四元数は回転を実行できるのだ。行列のユニタリ行列による変換だから、ユニタリ行列とその共役転置で両側から挟む形になるわけだ。

WPFのクオータニオンのコンストラクタ

リフレクタで見て見た。ちょっと修正して見やすくすると

public Quaternion(Vector3D axisOfRotation, double angleInDegrees)
{
    angleInDegrees = angleInDegrees % 360.0;
    double num2 = angleInDegrees * (2*Math.Pi/360.0);
    double length = axisOfRotation.Length;
    Vector3D vectord = (Vector3D) ((axisOfRotation / length) * Math.Sin(0.5 * num2));
    this._x = vectord.X;
    this._y = vectord.Y;
    this._z = vectord.Z;
    this._w = Math.Cos(0.5 * num2);
    this._isNotDistinguishedIdentity = true;
}

つまり、まず俺の計算はやっぱりどこか間違っていて、cos, sinの中は角度/2。それからこのコンストラクタは回転軸を正規化し、angleInDegreesは正確に回転角の意味。後、in degreesかよ！！

角速度ベクトル、回転時間が与えられた時の回転行列

対角化が面倒そうだったので、加法定理を求める方針でやってみる事にした。
つまり、

(I+δt(x Sx+y Sy+z Sz))(a0(t)I+a1(t)Sx+a2(t)Sy+a3(t)Sz)
=a0(t+δt)I+a1(t+δt)Sx+a2(t+δt)Sy+a3(t+δt)Sz ... (1)

を満たし、かつa0(0)=1, a1(0)=a2(0)=a3(0)=0であるような関数(連続性とか微分可能性とかの細かな扱いは知らない。必要なら行列の指数関数の話と併用すれば論理の隙間は埋められるだろう)を求めれば良い。それが求まるとa0(t)I+a1(t)Sx+a2(t)Sy+a3(t)Szが、角速度ベクトル(x,y,z), つまり、回転軸(x,y,z)/|(x,y,z)|, 角速度|(x,y,z)|でtだけ回転した時の回転行列、つまり、回転軸(x,y,z)/|(x,y,z)|, 回転角|(x,y,z)|tの回転行列を表す。

(x,y,z)をベクトルωで表し、(a1(t),a2(t),a3(t))をベクトルα(t)で表す事にする(ベクトルを明示するうまい方法がないのでここでだけギリシャ文字はベクトルを表すことにする)。しばしばtをオミットする。

すると、(1)式から
da0/dt=-ω・α
dα/dt=a0ω+ω×α
が得られる。
左辺にそれぞれa0, a1, a2, a3を書けて足しあげると0となる。これはa0^2+a1^2+a2^2+a3^2が不変である事を表す。t=0ではa0=1で他は0だから、連続性も考慮して、
a0=sqrt(1-|α|^2)
と書ける事が分かる。後はαの各成分だけ求めれば良い。αの微分方程式は

dα/dt=sqrt(1-|α|^2)ω+ω×α ... (2)

であるから、左からωを外積すると右辺第一項が消えて
d(ω×α)/dt=ω×(ω×α)
を得る。右辺はω×αに直交するから、|ω×α|は保存する(まじめにやるならば両辺とω×αの内積を取って右辺が消滅する事をがんばって計算すればd(|ω×α|^2)/dt=0を示せる。超めんどい。)。その値はt=0を参照すれば、0である。とすれば、ω×αは零ベクトルである。(2)式はいまや

dα/dt=sqrt(1-|α|^2)ω ... (3)

となる。αの各成分は初期値が等しく、微分係数が定数にωの各成分を乗じた形になっているので、αは全ての時刻でωに比例する。α=f(t)ω, とおくとf(t)が満たす微分方程式は
f(0)=0,
df/dt=sqrt(1-f^2|ω|^2)
である。必要ならばさらにg(t)=f(t)|ω|と変数変換しても良いが、ともかく、解が
f(t)=|ω|^(-1)sin|ω|t
と求まる。結局、

exp(t(x Sx+y Sy+z Sz))=cos(sqrt(x^2+y^2+z^2)t)
   +sin(sqrt(x^2+y^2+z^2)t)(x Sx+y Sy+z Sz)/sqrt(x^2+y^2+z^2)

と求まる。特に、sqrt(x^2+y^2+z^2)=1ならばexp(t(x Sx+y Sy+z Sz))は軸(x,y,z)周りの角度tだけの回転行列を表すが、それは、

exp(t(x Sx+y Sy+z Sz))=cos(t)I
   +sin(t)(x Sx+y Sy+z Sz)

である。

で、記憶によればこの結果は間違っていて、cos(t/2), sin(t/2)みたいになるはずなんだけど・・・結局、そのうち対角化の方もやってみなきゃならなそう。あと、そろそろWPFのQuaternionクラスのコンストラクタに色々なパラメタを入れたり、それからRotateTransform3Dを作って行列を調べたりして、ここまでの話と照らし合わせて見よう。あっちは、右手系だからちゃんと整合せいると思うんだけど・・・。Quaternionのコンストラクタには軸と角度を与えるコンストラクタがあるのだが、正規化されていないベクトルを与えた時の挙動が少し気になる。場合によっては、これらの引数は角速度ベクトルと回転時間としての役割を持つかもしれないし、内部で強制的に正規化されるかもしれない。この辺り全くドキュメント化されていない。

あ、足が、目が・・・

調子に乗ってたら著しく体調が悪くなってきた。最近、朝起きて数十分とその後との落差が酷い。

せっかくだから剛体回転でも、と思って準備

MSDNのクオータニオン関係の説明が詳細を欠いていて、係数の役割や回転変換との関係付けなど細かな点で不明な所があるので、剛体の回転も含めて整理しようとしたのだが、やっぱり集中力が阻害されて本をまともに読めない。まあ復習だしと思って自分ではじめからやってみる。

まず、三次元回転行列とは三次元の基底ベクトルを並べて作る。つまり、その行または列は正規直交する。従って回転行列RはRR^T=Iで特徴付けられる。

単位行列Iの近くの回転行列だけを考える。つまり、式(I+εM)(I+εM)^T=Iが、εを小さくしていった時にどんな条件式と等価になるかを考える。それは、
M=-M^T
である。この条件式はMの定数倍に対して保存される。つまり、Mに対して成り立てばεMに対しても成り立つ。つまり、単位行列Iの``近く''の回転行列の集合は
{I+M, where M=-M^T, Mは``小さい''}
と書ける。そこで、線形空間
L={M, where M=-M^T}
を考える事が重要になる。

反対称という事は、対角成分は0で上三角部分を決めると下三角部分が決まるから自由に決められるパラメタは三つしかない。これは、この線形空間が三次元空間である事を意味する。なんて勿体つけなくても、基底は直ぐに見つかって、

[0, 0,0]
[0, 0,1]
[0,-1,0]

[ 0,0,1]
[ 0,0,0]
[-1,0,0]

[ 0,1,0]
[-1,0,0]
[ 0,0,0]

がLを張る。

x,y,zが対称的に扱われるように、以下、右手系における右ネジ回転を基準とする事にする。これでこのエントリにおける不明瞭な点はなくなるだろう。後でWPFのクオータニオンと対応付ける必要はあるが。また、点の座標の変換と座標系の変換を明確に区別する。

(x,y,z)にある点をx軸でθ回転した時移る先の座標(x',y',z')は

[x'] [1,0  　,0   　]
[y']=[0,cosθ,-sinθ]
[z'] [0,sinθ,cosθ ]

と表される。この事から、x軸中心の微小なθ回転は

  [0,0,0]
I+[0,0,-θ]
  [0,θ,0]

と表される。他の軸中心の微小回転も同様に求められる。そこでLの基底を次のように取り直すと便利だ。

   [0,0, 0]
Sx=[0,0,-1]
   [0,1, 0]

   [ 0,0,1]
Sy=[ 0,0,0]
   [-1,0,0]

   [0,-1,0]
Sz=[1, 0,0]
   [0, 0,0]

この場合、積はIを加えれば閉じていて、反対称である。つまり
SxSx=-I, SxSy=Sz, SySx=-Sz, ... .

恒等変換の``近く''の回転変換を次々と作用させていく事で、恒等変換から``遠く''の回転変換も作る事ができる。剛体の回転運動はこのようにイメージできる。対応する概念は単位行列の``近く''の回転行列の複数回の積だ。つまり、x(t), y(t), z(t)を与えて、
Π(I+δt(x(t_n)Sx+y(t_n)Sy+z(t_n)S_z)), where δt=t_n-t_(n-1) ...(1)
によって単位行列から``連続に''別の回転行列に移る事ができる。ここでSx, S_y, S_zが単位ラジアンあたりの各軸周りの回転に対応していた事を思い出せば、x(t), y(t), z(t)は角速度ベクトルの各成分を与えると分かる。

この計算を実行するのは難しいが、ひとつだけ分かる事は、I, Sx, Sy, Szの積が閉じているために、この結果は必ず、I, Sx, Sy, Szの線形結合になるという事だ(Lの基底と違い、Iが加わる事に注意！)。これは重要な事だ。というのも、回転行列の積を繰り返すと計算誤差により徐々に回転行列の条件＝正規直交性を満たさなくなってしまうからだ。しかし、これら四つの基底に関する係数成分を時々刻々と計算するならば、計算誤差がどれだけあろうとも、それは回転変換に厳密な意味で対応し続ける！こうして四元数が登場する。

少なくとも剛体のシミュレーションに関する限り式(1)の計算結果を具体的に表す方法がない事はほとんど問題にならない。なぜならば、実際に、その時点での回転行列に恒等変換の近くの回転行列を乗じて、回転変換を更新する処理を繰り返せば良いからだ。変数に更新量を加算したり乗じたりしてループを回す類の処理は、シミュレーションにおける最も基本的なアルゴリズムで誰もが良く慣れているはずだ。また、すでに注意したように、回転行列自身の代わりに基底行列I, Sx, Sy, Szの係数を保存し、回転行列の積を計算する時には、基底行列の間の積の関係を使って係数の組から結果の係数を求めるようにすれば、係数に対応する回転行列は厳密に何らかの回転変換に対応し続ける。

存在する問題は二つだ。一つは、(1)式の積に現れる各項は近似的にしか回転行列に対応しないという事。もう一つは初期姿勢を決めるために(1)式の極限の具体的な結果が一つ欲しいという事。二つ目が解決されれば、その結果は、一つ目の問題の解決にも貢献するだろう。つまり、本物の微小な回転変換を繰り返し合成して一般の回転変換を作る事ができるようになる。

さて、(1)式の計算結果を具体的に求められる簡単なケースが一つある。それは、x(t), y(t), z(t)が一定である場合で、結果は行列の指数関数になる。その意味も明確だ。これは、角速度ベクトル[x, y, z]による、軸・角速度一定の回転を表す。特に、|[x,y,z]|=1ならば、時刻t=θにおいて角度θだけ回転した事になる。これがいわゆるクオータニオンと、軸／回転角で表した回転の関係だ。

ふう。次回は、せっかくだからSU(2)でやらずにこのまま気合で対角化して指数関数を求めるか。