SIMDライブラリの行列とクォータニオン

CPUが持っているベクトルユニットを陽に使うためのSIMDライブラリの一部が、Swift Standard Libraryの中に取り込まれていることをこないだ確認した。アセンブラを見ると作業量が少ない場合は、Standard Libraryにある関数を呼ぶのではなく、ベクトルユニットの命令が直接インライン展開されていることがわかる。ソースは確認してないので実際にどうなのかはわからないけど、@inlinable属性がついていてその通りのコードが出力されているということだろう。

またこのStandard Libraryとは別に

import simd

で、SIMD関連の追加がされることを知った。simdのインポートによって4次元までのベクトルと行列の基本的な演算がいろいろ可能になる。ちょっとだけ前回の追加．．．

追加の一つが行列演算である。行列は

public struct simd_floatNxM
public struct simd_doubleNxM

のかっこうで定義されている。Nが列（Column）の数でMが行（Row）の数でそれぞれ2、3、4の場合が網羅されている。たとえば普通に行列に書いたとき \begin{equation} \begin{pmatrix} 0 & 1 & 2\\ 3 & 4 & 5 \\ \end{pmatrix} \nonumber \end{equation} のような行列は

simd_double3x2

になる。ただし要素はColumn majorで保持されるので、上の行列を要素を列挙して初期化するとき

let m = simd_double3x2([[0, 3], [1, 4], [2, 5]])

あるいは、行方向であることを明示して

let m = simd_double3x2(rows: [[0, 1, 2], [3, 4, 5]])

になる。配列の形があっていないと実行時エラーが起きるようである。

普通のベクトル行列の和や積は演算子として定義されてるけど、Column majorなので、例えばベクトルとの積 \begin{equation} \begin{pmatrix} 1 & 2 & 3 \\ \end{pmatrix} \begin{pmatrix} 0 & 1 & 2 \\ 3 & 4 & 5 \\ 6 & 7 & 8 \\ \end{pmatrix} \end{equation}

    public static func * (lhs: double3, rhs: simd_double3x3) -> double3

と \begin{equation} \begin{pmatrix} 0 & 3 & 6 \\ 1 & 4 & 7 \\ 2 & 5 & 8 \\ \end{pmatrix} \begin{pmatrix} 1 \\ 2 \\ 3 \\ \end{pmatrix} \end{equation}

    public static func * (lhs: simd_double3x3, rhs: double2x3) -> double2x3

とでは効率が違うことに注意する必要がある（言うまでもないけど(1)の方が速い）。

それ以外にも正方行列にはプロパティとして

    public var determinant: Double { get }
    public var transpose: double4x4 { get }
    public var inverse: simd_double4x4 { get }

なんかが定義されている（前回示したように関数もあるけど、プロパティの方を使うようにというコメントがあった）。逆行列が存在しないときinverseプロパティは要素全てがNaNの行列が返されるようである。こういうのはオプショナルになってる方が使いやすかったけど、実態はC関数なのでしょうがないかもしれない。



追加のもう一つがクォータニオン（Quaternion）である。クォータニオンは3次元空間内の任意の回転をエレガントかつ簡単に表すことができるので、3Dグラフィクスではよく使われる。SIMDライブラリでは

let ax = double3(arrayLiteral: 0.0, 0.0, 1.0)
let quat = simd_quatd(angle: .pi * 0.5, axis: ax)

のように、回転させる軸（上の例では3次元ベクトルax）と回転角（同じく.pi * 0.5つまり$\pi/2$）を与えることでクォータニオン構造体を作ることができる（ほかにも色々なイニシアライザが用意されている）。そしてベクトルを実際に回転させるのは

let v0 = double3(arrayLiteral: 1.0, 0.0, 0.0)
let v1 = quat.act(v0)

とする。この例では$x$方向を向いたベクトルv0を$z$軸の周りに90°回転させたものがv1になる。

SIMDのSwiftインターフェイスは演算子が定義されていて、ふたつのクォータニオンの合成は

let ax0 = double3(arrayLiteral: 0.0, 0.0, 1.0)
let quat0 = simd_quatd(angle: .pi * 0.5, axis: ax0)
let ax1 = double3(arrayLiteral: 0.0, 1.0, 0.0)
let quat1 = simd_quatd(angle: .pi * 0.5, axis: ax1)

let quat2 = quat0 * quat1

のように積の形で書けるようになっている。ここで注意すべき点があって、クォータニオンの積は可換ではない。上の場合$z$軸周りの回転の後に$y$軸周りの回転をほどこすことになる。SIMDライブラリはColumn Major、つまり行列に限らずクォータニオンも積の右の方があとの作用になる。それとクォータニオンを回転軸と角度で作った場合、回転軸ベクトルが規格化（長さが1）されていないと正しいものが返らないようなので注意が必要である（回転軸ベクトルはsimd_normalize関数を呼んで規格化しておく）。

クォータニオンにまつわる一連の作業で、クォータニオンが要素を4つ持つベクトルだなんてことはあまり気にする必要がない。そこが便利なところである。その背景には複素数の積の演算が複素平面上の回転になるのを3次元に拡張するとクォータニオンの積が3次元での回転になる、という数学的にシンプルな対応があるけど、そんなこと知ったこっちゃない、と言う人でも使うことができるというのもその懐深さを表しているようで楽しい。

実数の世界での解析学が複素数に拡張されたとき、実数だけでは思ってもみなかった強力な定理がたくさん証明された。複素関数論の技術の領域への影響はつくづく巨大だと思い知らされることがよくある。複素数をさらに拡張したクォータニオンは解析学をさらに便利なものにするか、とクォータニオンを作った張本人のハミルトンは考えたらしい。残念ながら普通の複素数以上の解析学に追加する定理を見つけることができなかった。

ハミルトンの思惑は空振りだったけど100年経ってコンピュータグラフィクスに便利に使われることになった。よかったじゃん、ハミルトンさん。ちなみに、光学の教科書にときどき出てくるアイコナール方程式はハミルトンさんのアイデアである。ハミルトニアンと解析理学といい、ハミルトンさんは解析と代数の両方のセンスを持って間をつなぐような業績が多い。カッコいい。



脱線したのでSIMDにあるクォータニオンの話に戻す。光学でレンズ設計をするとき、光線追跡という一種のシミュレーションを行う。そこでは計算量を減らすために座標変換をたくさんする。クォータニオンは3次元空間内の任意の回転を表すことができるので光線追跡にも使える。しかし全くランダムな光線束を使うことは少なくて、対称性を利用して光線の本数を減らすことが普通である。その場合、クォータニオンよりも冗長性の多い同次座標変換を使うことの方が多い。その方が座標軸ごとに分けて考えやすいからである。

SIMDライブラリにはクォータニオンを座標変換行列に変換する関数が用意されている。

let mat3 = simd_double3x3(quat);
let mat4 = simd_double4x4(quat);

simd_double?x?は正方行列のクォータニオンを引数にしてイニシアライザを呼ぶ。その結果mat3は回転の行列、mat4は平行移動成分が0の同次座標行列が返る。

その逆も

let quat3 = simd_quatd(mat3)
let quat4 = simd_quatd(mat4)

の両方のイニシアライザが定義されている。自分でやろうとすると結構めんどくさくて間違いやすいので小さな親切である。

ただし、これで得られる行列もColumn majorである。例えば上のイニシアライザでできる同次座標変換の行列は \begin{equation} \left( \begin{array}{ccc|c} &&&0\\ &A&&0\\ &&&0\\ \hline t_x & t_y & t_z &1 \end{array} \right) \nonumber \end{equation} の形（$A$は3次元での回転の行列）で、変換したいベクトルへは右側から掛け算（ベクトルが左に）する。座標変換の合成も右側が後になるように記述する。3Dグラフィクスの世界ではこちらの方が主流かもしれない。

勉強も兼ねてと思って自前のベクトルと行列の定義を使って光線追跡エンジンを書き始めたけど、やっぱりSIMDライブラリを使う方法に書き直そう。ひっそりとではあるけどわざわざ用意してくれてるわけだし。しかしまさか、次のバージョンでobsoleteになる、なんてことはないだろうな。