Zenn
🌀

Rotations.jlで回転しよう！

2021/12/01に公開
Julia
math
tech

回転の導入

「回転」というのは直観的なイメージ通り、物体の形状を変化させずに向きを変える変換のことです。

線形代数では、回転を直交行列で表すことができます。
ただし正確に言えば、鏡像は回転に含めたくないので行列式が1のものに限定します。
つまり

  • R^{-1} = R^{\top}
  • \det(R) = 1

を満たす行列Rを回転行列と呼びます。
このような行列全体は群をなすので、記号SO(n)が使われます。特殊直交群と呼ばれます。
ここでのnは正方行列の一辺のサイズで、回転が行われる空間の次元に相当します。

鏡像も含むもの、つまりn次直交行列全体は記号O(n)で表されますが、以降の本記事では登場しません。

直交行列を複素数に拡張したものとして、ユニタリ行列があります。転置{R}^{\top}の代わりに随伴{R}^{*}(転置+複素共役)を考えます。回転行列と同様に行列式が1のものを考えると都合が良いので

  • R^{-1} = R^{*}
  • \det(R) = 1

を満たす行列全体を考えて、特殊ユニタリ群SU(n)と呼びます。

行列式の制限がないもの、つまりn次ユニタリ行列行列全体は記号O(n)で表されますが、こちらも以降の本記事では登場しません。

私の過去の関連発表

回転については何度かLTなどで発表しています。
本記事を読む際にも以下が部分的に役立つと思います。

  • 第1回日曜数学会
    • 日程: 2015/06/20
    • 「逆数の作図からCayley変換まで」というタイトルでLTしました。
    • SO(3)とSU(2)が準同型であることを、一次分数変換を通じて確認しました。
    • 発表資料
  • 第4回日曜数学会
    • 日程: 2016/01/30
    • 「Riemann球面に内接する直方体」というタイトルでLTしました。
    • Riemann球面での対蹠点がR\in SU(2)の不動点であることを確認しました。
    • 発表資料, 発表の様子
  • JuliaTokai #10
    • 日程: 2021/05/15
    • 「Rotations.jlで学ぶ3次元の回転」というタイトルで発表しました。
    • 本記事の内容に最も近いです。
    • 発表資料, イベントページ

次元について

回転が作用する空間の次元

  • n=1は直線なので回転行列は単位行列しかないので自明です。
  • n=2は原点中心の角度\theta回転ですべての回転行列を表せます。
  • n=3より上では回転は非可換で、紛らわしくなってきます。
  • n\ge 4は応用が少ないです。私達の住んでいる空間は3次元だからです。

そういう訳で、応用上特に重要で紛らわしいのがn=3です。

2次元と3次元の回転[1]を扱うための色々な関数を用意しているJuliaパッケージ、それがRotations.jlです！

SO(n)の次元

SO(n)は群であると書きましたが、同時に多様体で、つまりLie群です。
特殊直交群の次元は一般に\dim(SO(n)) = n(n-1)/2なので\dim(SO(3)) = 3ですね。

つまり、局所的には3つのパラメータで3次回転行列を表すことができるということになります。たとえば、軸回り回転の合成で表すEuler角は3回の回転なのでパラメータ3つ、局所的に過不足ないという訳です。
これは行列の要素数9=3\times 3に比べればかなり少ないですね。

2次元でも同様で、\dim(SO(2)) = 1なのでこちらも行列の要素数4=2\times2に比べて少ないですね。

数値的に回転を扱うときの注意点

簡単のために、2次元で確認しましょう。
適当に回転行列を作って行列式を計算すると以下のようになります。

julia> using LinearAlgebra

julia> a = 0.1;

julia> R = [cos(a) -sin(a);sin(a) cos(a)]
2×2 Matrix{Float64}:
 0.995004   -0.0998334
 0.0998334   0.995004

julia> det(R)
1.0000000000000002

行列式は1のはずなのに…厳密に1にならないですね。
これはコンピュータが浮動小数点数を使っているために起こっています。

誤差は累積していくので、べき乗を繰り返すとどんどんずれていきます。

julia> det(R)
1.0000000000000002

julia> det(R^100)
1.0000000000000133

julia> det(R^10000)
1.000000000001315

julia> det(R^1000000)
1.0000000001314495

julia> det(R^100000000)
1.000000013144949

julia> det(R^10000000000)
1.000001314495683

julia> det(R^1000000000000)
1.0001314581217502

julia> det(R^100000000000000)
1.013231722819454

このような誤差の問題を避けるためには、SO(n)の性質をうまく使い、より少ないパラメータで回転行列を表現すると良いです。

パラメータ3つで回転を表す場合は、9次元空間\mathbb{R}^{n^2}に埋め込んで得られる(冗長な)座標の代わりに、SO(n)に直接座標を入れることに相当します。
浮動小数点数を使うので、依然としてSO(n)内での誤差は発生しますが、SO(n)からはみ出すことは無くなります。

また、この方針は実行速度・メモリ使用量の観点からも効果的です。
パラメータが少なくなるので計算量が減り、同時に必要なメモリも削減できます。

では、どのようにしてパラメータを減らすのが適切でしょうか？

  • SO(2)ではパラメータとして回転角を取れば良いので、単純ですね。
  • SO(3)では少し込み入っています。色々な便利なパラメータがあり、状況に合わせて使い分ける必要があります。次節で簡単に紹介します。

忙しい人のための3次元回転パラメータ入門

詳しい話は本記事末尾の参考文献を参照してください。
概要だけ述べます。

  • 回転行列SO(3)
    • 3\times3実行列。パラメータ9つ。
  • ユニタリ行列SU(2)
    • 2\times2複素行列。パラメータ4つ。
    • 一次分数変換(Möbius変換)でSO(3)と2:1で対応。
  • 単位四元数
    • 単位長さの四元数q=w+ix+jy+kz。パラメータ4つ。
    • SU(2)と代数的に同型。
    • 3次元球面S^3と幾何的に同型。
    • 一次分数変換とSU(2)を経由せずに、直接SO(3)と対応付ける説明もある。
  • MRP (Modified Rodrigues Parameters)
    • S^3の局所座標。パラメータ3つ。
    • S^3から南極(w=-1)を除いて\mathbb{R}^3へ立体射影で対応させる。
  • Rodrigues Parameters
    • S^3の局所座標。パラメータ3つ。
    • S^3から赤道(w=0)を除いて\mathbb{R}^3へ球面原点からの射影で対応させる。
  • 軸回り回転
    • 単位長さ回転軸と回転角。パラメータ4つ。
    • 任意の3次回転行列は適当な回転軸回りの回転になるので可能。
  • 軸回り回転2
    • 回転軸×回転角。パラメータ3つ。
    • 回転軸を表す単位長さベクトルに回転角を掛ければパラメータを減らせる！
    • Lie代数\mathfrak{so}(3)と相性が良い。
  • Euler角
    • 座標軸まわりの3回の回転で回転行列を表す。パラメータ3つ。
    • 直感的だが、ジンバルロック・非可換性等の影響があり使いやすくはない。

Rotations.jlの方法

さて、ここからが本題です。
以下のコード例では、Rotations.jlがインストールされていてusing Rotationsされているとします。

回転は実行列として扱われる

julia> Rotation <: AbstractMatrix
true

これまでの説明からすれば、3次元回転を3×3行列で表すことは自然で、上のtrueは当り前のように見えます。しかし…

  • 「単位四元数(\simeq SU(2))こそ回転」と考える人もいます。
    • AbstractQuaternionなる抽象型は存在しないので、この場合はRotation <: Numberになります。
    • 単位四元数の方が都合が良いのは「単連結になっているから」等いくつか理由があります。
  • 「回転という操作は幾何的なもので、行列に変換可能であっても行列として扱うべきでない」という人もいます。
    • この場合はRotation <: Anyのみの階層関係です。
    • scipyの回転はこれに近いですね。

Rotations.jlで回転を行列として扱うのは、行列としての演算(method)をそのまま使えると便利だからです。

色々なパラメータに合わせた具象型が存在

Rotations.jlでの型の階層関係は以下のようになっています。

julia> using StaticArrays, Rotations

julia> Rotation <: StaticMatrix <: AbstractMatrix
true

julia> Rotation{2,Float64} <: StaticMatrix{2,2,Float64} <: AbstractMatrix{Float64}
true

julia> Rotation{3,Float64} <: StaticMatrix{3,3,Float64} <: AbstractMatrix{Float64}
true

Rotation{3}は3次元の回転を表すための抽象型で、以下のような具象型をsubtypeを持ちます。
ちょうど前述のパラメータのそれぞれに対応するようになっています。[2]

  • 回転行列SO(3) RotMatrix{3}
  • ユニタリ行列SU(2) (単位四元数と同一なので存在せず)
  • 単位四元数 QuatRotation
  • MRP (Modified Rodrigues Parameters) MRP
  • Rodrigues Parameters RodriguesParam
  • 軸回り回転 AngleAxis
  • 軸回り回転2 RotationVec
  • Euler角 RotXYZ, RotXYX etc.

Rotation{2}は単純で、以下の2種類があります。

  • RotMatrix{2}: 内部的に2\times 2のSMatrixを保持するもの
  • Angle2d: 内部的に回転角のみを保持するもの

これらの具象型は以下のようにインスタンスを生成して変換することができます。

julia> r = RotX(1.2)  # x軸回りで1.2(rad)回転
3×3 RotX{Float64} with indices SOneTo(3)×SOneTo(3)(1.2):
 1.0  0.0        0.0
 0.0  0.362358  -0.932039
 0.0  0.932039   0.362358

julia> AngleAxis(r)  # 軸回りの回転に変換
3×3 AngleAxis{Float64} with indices SOneTo(3)×SOneTo(3)(1.2, 1.0, 0.0, 0.0):
 1.0  0.0        0.0
 0.0  0.362358  -0.932039
 0.0  0.932039   0.362358

julia> QuatRotation(r)  # 四元数回転に変換
3×3 QuatRotation{Float64} with indices SOneTo(3)×SOneTo(3)(Quaternion{Float64}(0.825336, 0.564642, 0.0, 0.0, true)):
 1.0  0.0        0.0
 0.0  0.362358  -0.932039
 0.0  0.932039   0.362358

julia> Rotations.params(QuatRotation(r))
4-element SVector{4, Float64} with indices SOneTo(4):
 0.8253356149096782
 0.5646424733950354
 0.0
 0.0

剛体回転の例

Rotations.jlの応用として、剛体回転を実装してみましょう。
テニスラケットの定理を数値計算で確かめたかったのですが、実装が間に合いませんでした！
(あとで追記します)

今後のRotations.jl

Rotations.jlは少し歴史の長いパッケージで、Julia v0.4のころから開発が始まっていたようです。私がRotations.jlの開発に加わったのは今年に入ってからで、大きなところだと以下のPRを行いました。

  • rand関数のバグの修正(#151)
    • 以前のバージョンではランダムな回転を正しく生成できていませんでした。
  • Documenter.jlを使ったドキュメント作成(#169)
    • ドキュメントはこちらです。
    • ロゴも私が作成しました。(#172)
    • ドキュメントが埋まっていない部分が残っています。contributionチャンスです！
  • Quaternions.jlとの互換性の向上(#175)
    • 四元数を扱うためのパッケージQuaternions.jlとの変換ができるようになりました。

まだまだ開発は続いていて、将来的には以下の機能を追加する予定です。

  • RotationGeneratorの追加
    • 回転の生成子を扱うための抽象型(反対称行列)
    • v1.2.0で追加予定です。(#203)
    • 回転行列のlogを取ってRotationGeneratorを返すようになります。
    • exp(::RotationGenerator)で回転行列が得られるようになります。
  • ChainRulesCore.jlとの連携
    • v1.3.0で追加予定です。
    • Zygote.jlなどの自動微分パッケージを適用するために修正する予定です。

参考文献

脚注

  1. 将来的には4次元以上にも対応したいと考えていますが、未実装です。 ↩︎

  2. 一覧はsubtypes(Rotation{3})で取得できます。 ↩︎

GitHubで編集を提案

Discussion