本記事はJuliaLang Advent Calendar 2022の12日目の記事です!
2022年にメンテナンスしていたQuaternions.jlについて書きます！

はじめに

四元数(quaternion)とは $q = w+ix+jy+kz$ のように表される数のことで、複素数を拡張したものと考えることができます。四元数全体を記号 $\mathbb{H}$ で表します。
ここで $(i,j,k)$ は虚数単位に相当するもので、以下の演算規則を充たします。

\begin{gathered} ii = jj = kk = -1 \\ ij = k, \quad jk = i, \quad ki = j \\ \end{gathered}

あとは通常の分配法則などを使えば和や積が計算できます。^[1]

\begin{aligned} q_1 + q_2 &= (w_1+ix_1+jy_1+kz_1) + (w_2+ix_2+jy_2+kz_2) \\ &= (w_1+w_2) + i(x_1 + x_2) + j(y_1 + y_2) + k(z_1 + z_2) \\ \\ q_1q_2 &= (w_1+ix_1+jy_1+kz_1) (w_2+ix_2+jy_2+kz_2) \\ &= (w_1w_2 - x_1x_2 - y_1y_2 - z_1z_2) +i(w_1x_2 + x_1w_2 + y_1z_2 - z_1y_2) +j(w_1y_2 + y_1w_2 + z_1x_2 - x_1z_2) +k(w_1z_2 + z_1w_2 + x_1y_2 - y_1x_2) \end{aligned}

これを計算できるJuliaのパッケージがQuaternions.jlです！
インストールするにはJuliaのREPLで

]add Quaternions

を実行します。
最初の計算例として以下を実行しましょう。

julia> using Quaternions

julia> Quaternion(1,2,3,4) + Quaternion(5,6,7,8)  # 足し算
Quaternion{Int64}(6, 8, 10, 12)

julia> Quaternion(1,2,3,4) * Quaternion(5,6,7,8)  # 掛け算
Quaternion{Int64}(-60, 12, 30, 24)

julia> sin(Quaternion(1,2,3,4))  # 解析関数もOK
QuaternionF64(91.78371578403465, 21.886486853029176, 32.82973027954377, 43.77297370605835)

より詳細な例などはQuaternions.jlのドキュメントも参照してください！

Julia言語コニュニティにおける四元数の利用状況

Julia言語のパッケージの仕組みによってコードの再利用が推奨されていますが、多くのパッケージではQuaternions.jlに依存せずに独自に四元数の型を定義しています。その例：

四元数の実装を目的にしたパッケージもいくつか存在します。

このような車輪の再実装を招いてしまった状況の原因は、Quaternions.jlが適切にメンテナンスされていなかったことにあります。Quaternions.jlは

パッケージ名が適切すぎること
JuliaGeometryのOrganizationでメンテナンスされていること

の性質を充たしており、他のパッケージで決定打になるものが出にくかったこともあります。

2022年2月くらいからSethさんと私がメンテナ権限を貰ってQuaternions.jlのメンテナンスを続けてきました。
v0.7.0のリリース辺りからQuaternions.jlはようやく使える形に整ってきた感じがあります。
本記事ではQuaternions.jlのメンテナンスでの変更点や教訓などを記述しようと思います。

メンテナンスでの変更点

私が開発に参加したのはv0.4.6辺りからです。最新リリースがv0.7.1なので、破壊的変更を伴ったリリースを3回行ったことになります。^[2]
以降ではこれらの変更の概要を説明します。

`norm`フラグの削除

四元数の一番の応用は3次元回転で、その用途に使える四元数は絶対値が1のものに限ります。
数学的に言えば、次の群 $U(1,\mathbb{H})$ が $SO(3)$ に準同型になっていることが役立っています。

\left\{w + ix + jy + kz \in \mathbb{H} \ | \ x, y, z \in \mathbb{R} \right\} = U(1,\mathbb{H}) \simeq S^3

従来は四元数の型Quaternionが以下のように実装されていました。

struct Quaternion{T<:Real} <: Number
    s::T
    v1::T
    v2::T
    v3::T
    norm::Bool
end

ここで、フィールドnorm^[4]は「四元数の絶対値が1であるか」を表すフラグで、現在の実装ではこのnormは削除されています。
削除された理由はいくつかありますが、代表的なものは以下です。

長さ1でない四元数を使う場合に不要なフィールドを避けたい
単位四元数を扱う際に必須という訳ではない
余分なフィールドによってパフォーマンスが低下する
実装次第でQuaternion(2,0,0,0,true)のような無意味なフラグを生成してしまう

詳細な議論はissue#60をご覧ください。

教訓

必須ではないフィールドを構造体に入れるのはやめましょう。
Baseの実装に類似するものがあれば合わせましょう。(今回はBase.Complex)

`Octonion`の削除

実数→複素数→四元数の拡張の延長に八元数があります。
Octonionはこのに八元数を表す型で、従来はQuaternions.jlで提供されていました。

これを削除した理由は以下です。

1つのパッケージはなるべく軽量であるべき。特にQuaternions.jlのようなプリミティブなものは。
Octonionを提供したいなら別のOctonions.jlパッケージを用意するべき。
四元数と八元数を同時に使いたいことは無いはず

そういう訳で、Octonions.jlが用意された後にOctonionsは削除になりました。

教訓

パッケージはなるべく小さく作りましょう。
「数学的な類似性があるから」というのは同一のパッケージで提供する理由には弱く、実用上の側面から考える方が良いでしょう。

`DualQuaternion`の削除

二重四元数は、二重数と四元数を組み合わせたもので、剛体変換の記述に便利です。
数学的な背景についてはTokoroさんのブログ記事やQuaternions.jlのドキュメントを参照ください。

このような二重四元数を扱う型として、従来はDualQuaternionが提供されていました。
これが削除された理由は、単純にQuaternion{ForwardDiff.Dual}が代わりになったからです。

教訓

Juliaの型パラメータは便利で、それで代用できるものがあれば使いましょう。

`Complex`との互換性の削除

従来はComplexとQuaternionの変換として $\mathbb{C}\to\mathbb{H} ; x+iy \mapsto x+iy+0j+0k$ が以下のように（中途半端に）実装されていました。

julia> using Quaternions  # v0.6.0

julia> quat(1,2,0,0) + Complex(1,2)  # 足し算はtype promotionされる
Quaternion{Int64}(2, 4, 0, 0, false)

julia> convert(Complex, quat(1,2,0,0))  # しかし明示的な型変換はエラーになる
ERROR: MethodError: no method matching Complex(::Quaternion{Int64})
Closest candidates are:
  (::Type{T})(::T) where T<:Number at boot.jl:772
  (::Type{T})(::AbstractChar) where T<:Union{AbstractChar, Number} at char.jl:50
  (::Type{T})(::Base.TwicePrecision) where T<:Number at twiceprecision.jl:266
  ...

現在はこのような互換性が削除されています。その理由：

写像 $\mathbb{C}\to\mathbb{H}$ として自然なものを1つ選ぶことはできない
ComplexとQuaternionsの互換性は実用的に便利という訳ではない

教訓

中途半端な実装はやめましょう。
テストはちゃんと書いておきましょう。
数学的な妥当性の低い実装は避けましょう。

他の変更事項

LinearAlgebra.normalize(::Quaternion)の削除
- normalize(::Number)は実装されておらず、sign(::Number)を使うのが一般的
- 教訓: 既存の関数にメソッドを追加する場合は元々のメソッドと類似性があるようにしましょう。^[5]
linpolの削除
- 同一機能のslerpが実装されていた
- Slerp(spherical linear interpolation)の方がよく使わる名称
- 教訓：同じ機能を別の関数で提供しないようにしましょう。
rotationmatrixなどの削除
- 従来は四元数から回転行列に変換( $U(1,\mathbb{H}) \to SO(3)$ )するための関数rotationmatrixが用意されていた
- この関数は戻り値がMatrix{Float64}であってStaticMatrixでないので非効率。
- 「Quaternions.jlは四元数に関連するプリミティブな演算のみ提供すべし」と考えて削除
- 教訓: Quaternions.jlのような基礎的なパッケージは軽い方が良く、複雑な機能は別のパッケージに切り分けた方が他の人からも使いやすいです。^[6]
テストの追加
- Coverage: 43% → 99% 🎉🎉
ドキュメントの作成
- Documenter.jlでドキュメントを書きました。
- 便利で良いですね。
- ロゴも作りました。

おわりに

本記事で書いた内容の要約はJulia Discourseでのアナウンスにもあります。

今後の開発予定としては以下があります。

ChainRules.jlへの対応
v1.0.0のリリース
Quaternions.jlに依存してもらえるよう他のパッケージにPRを作成

開発に協力して下さる方を引き続き募集しておりますので、興味がある人はQuaternions.jl/issuesを覗いてみて下さい！

脚注

積が可換でないことに注意して下さい。 ↩︎
Semantic versioningです。v0.Y.ZではYのincrementで破壊的変更、Zのincrementで機能追加・バグ修正を表します。(x-ref: https://twitter.com/Hyrodium/status/1536769231463022594?s=20&t=1oKHaAmXu2_nP6yRIsJpXQ) ↩︎
自分の書いた実装の問題点は、間違いにどの程度一般性があるのか不明瞭で記事として書きにくいですが、他の人の書いたコードであれば書きやすいというのもあります…。 ↩︎
ノルムの値(実数)ではなくis_normalizedの略と考えた方が良いです。trueで単位四元数を表します。 ↩︎
元々の関数と異なる目的のメソッドを追加する行為はType III piracyと呼ばれます。詳細はHands-On Design Patterns and Best Practices with Juliaもご覧ください。 ↩︎
StaticMatrixとして回転行列を表すパッケージとしてRotations.jlがあります。こちらも私がメンテナンスしています。 ↩︎

Discussion