回転行列と四元数の計算速度比較

目的

ベクトル，回転行列，座標変換行列といった三次元幾何演算をpythonのリストをベースにgeo.pyというモジュールを自作している．

最近，回転の表現に四元数(QUATERNION)がよく用いられているが，回転行列(MATRIX)と比較して本当にメリットがあるのかを，演算速度の点で比較してみることにした．

ちなみに，積・和の数の比較で言うと，

	積	和
回転行列の積	27	18
四元数の積	16	12
回転行列による<html><br></html>ベクトルの回転	9	6
四元数による<html><br></html>ベクトルの回転<html><br></html>下は実部計算を除く	32 <html><br></html> (28)	24 <html><br></html> (20)

であり，四元数に計算上のメリットはほとんど見られない．

このテストは以下の条件で行った．

ProBook 474s
メモリ：8 GB
CPU：Core™ i5-3230M
OS: Ubuntu 20.04
jupyter notebook，python3

このipynb自身は，matrix_vs_quaternion.ipynbとなっている．

モジュールの読み込み

geo.py は自作モジュール. 座標の姿勢（回転）に関して直交行列(MATRIX)と回転四元数(QUATERNION)の両方をサポートしている．

from geo import *
import time
import pandas as pd

データの生成

geo.pyには三次元ベクトル：VECTOR，三次元回転行列：MATRIX，四元数：QUATERNIONのクラスが定義されている． MATRIX,QUATERNIONのa, bはそれぞれx軸，y軸周りの回転を指定する．

v1=VECTOR(1,2,3)
R1=MATRIX(a=pi/3)
R2=MATRIX(b=pi/6)
q1=QUATERNION(a=pi/3)
q2=QUATERNION(b=pi/6)
qv1=QUATERNION([0,v1])

この内容はこうなる．

print('v1 =', v1)
print('R1 =', R1)
print('R2 =', R2)
print('q1 =', q1)
print('q2 =', q2)
print('qv1 =', qv1)

v1 = v:[1.0, 2.0, 3.0]
R1 = m:[[1.0, 0.0, 0.0], [0.0, 0.5000000000000001, -0.8660254037844386], [0.0, 0.8660254037844386, 0.5000000000000001]]
R2 = m:[[0.8660254037844387, 0.0, 0.49999999999999994], [0.0, 1.0, 0.0], [-0.49999999999999994, 0.0, 0.8660254037844387]]
q1 = q:(0.8660254037844387, v:[0.49999999999999994, 0.0, 0.0])
q2 = q:(0.9659258262890683, v:[0.0, 0.25881904510252074, 0.0])
qv1 = q:(0, v:[1.0, 2.0, 3.0])

時間計測の関数と結果保存データ

def test(n,fn):
    i=0
    start=time.time()
    while i< n :
        fn()
        i += 1
    end = time.time()
    rslt=end-start
    return rslt

data = []

def judge(test_name, m_time, q_time) :
    if m_time < q_time :
        judgment = "MATRIXの勝ち"
    elif m_time > q_time :
        judgment = "QUATERNIONの勝ち"
    else :
        judgment = "引き分け"
    return test_name, m_time, q_time, judgment

ループ回数の決定

test(100, lambda : R1*R2)

0.0013973712921142578

test(1000, lambda : R1*R2)

0.015026330947875977

test(10000, lambda : R1*R2)

0.11950469017028809

test(100000, lambda : R1*R2)

0.7195339202880859

test(1000000, lambda : R1*R2)

9.323349952697754

test(10000000, lambda : R1*R2)

78.96706867218018

百万回ぐらいでループ前後のオーバーヘッドの影響が少なくなってきている．まだ多少影響はあるが，一千万，一億は時間がかかるし，どうせループ内の処理の影響は消せないので百万回に決定する．

N=1000000

オブジェクトの生成

m_time = test(N, lambda : MATRIX())
print(m_time)

1.5926079750061035

q_time = test(N, lambda : QUATERNION())
print(q_time)

1.4729571342468262

data.append(judge('オブジェクトの生成', m_time, q_time))

回転の合成　QUATERNIONの勝ち

R1*R2

m:[[0.8660254037844387, 0.0, 0.49999999999999994], [0.43301270189221924, 0.5000000000000001, -0.75], [-0.25, 0.8660254037844386, 0.43301270189221946]]

(R1*R2).quaternion()

q:(0.8365163037378079, v:[0.4829629131445341, 0.22414386804201336, 0.1294095225512603])

q1*q2

q:(0.8365163037378079, v:[0.4829629131445341, 0.2241438680420134, 0.12940952255126034])

m_time = test(N, lambda : R1*R2)
print(m_time)

8.527657508850098

q_time = test(N, lambda : q1*q2)
print(q_time)

7.0049989223480225

data.append(judge('回転の合成', m_time, q_time))

ベクトルの回転　MATRIXの圧勝

R1*v1

v:[1.0, -1.5980762113533158, 3.2320508075688776]

q1*qv1*q1.conjugate()

q:(0.0, v:[1.0, -1.5980762113533153, 3.2320508075688776])

m_time = test(N, lambda : R1*v1)
print(m_time)

3.5639126300811768

q1c=q1.conjugate()
q_time = test(N, lambda : q1*qv1*q1c)
print(q_time)

13.009745597839355

data.append(judge('ベクトルの回転', m_time, q_time))

結論

まとめの表

df = pd.DataFrame(data, columns=["項目", "MATRIX", "QUATERNION", "結果"])
df

  .dataframe tbody tr th:only-of-type {
      vertical-align: middle;
  }

  .dataframe tbody tr th {
      vertical-align: top;
  }

  .dataframe thead th {
      text-align: right;
  }

</style>

<thead> <tr style=“text-align: right;”> <th> </th> <th> </HTML> 項目 <HTML> </th> <th> </HTML> MATRIX <HTML> </th> <th> </HTML> QUATERNION <HTML> </th> <th> </HTML> 結果 <HTML> </th> </tr> </thead> <tbody> <tr> <th> </HTML> 0 <HTML> </th> <td> </HTML> オブジェクトの生成 <HTML> </td> <td> </HTML> 1.592608 <HTML> </td> <td> </HTML> 1.472957 <HTML> </td> <td> </HTML> QUATERNIONの勝ち <HTML> </td> </tr> <tr> <th> </HTML> 1 <HTML> </th> <td> </HTML> 回転の合成 <HTML> </td> <td> </HTML> 8.527658 <HTML> </td> <td> </HTML> 7.004999 <HTML> </td> <td> </HTML> QUATERNIONの勝ち <HTML> </td> </tr> <tr> <th> </HTML> 2 <HTML> </th> <td> </HTML> ベクトルの回転 <HTML> </td> <td> </HTML> 3.563913 <HTML> </td> <td> </HTML> 13.009746 <HTML> </td> <td> </HTML> MATRIXの勝ち <HTML> </td> </tr> </tbody>

</div> </HTML> 結論から言うと，ほぼ予想通りである．

今回の実装では四元数の虚部をVETCORクラスとして定義しているため，そのオーバヘッドが少しあるように見える．

ベクトルの回転では，共役四元数を都度生成するのは控えた．多数のベクトルに同じ回転を施すイメージである． 3DCGなどの回転ではこのような使われ方が多いだろう．ロボットの場合は必ずしもそれは当てはまらないかもしれない．

この手のプログラムは実装の影響が大きいのは致し方ない．まだそれぞれ工夫の余地があると思う．

当初の実装では四元数の乗算を見た目のわかりやすさからベクトルのスカラー倍，ベクトルの外積を用いていたが，それだとVECTORオブジェクトの生成が3回余分に入るため，回転行列の乗算より遅くなってしまっていた．

ここでの比較ではそこはVECTORの生成を行わない形に書き下している．

目次