articles:matrix_vs_quaternion_1 [RoboLab.Suehiro]

この文書は読取専用です。文書のソースを閲覧することは可能ですが、変更はできません。もし変更したい場合は管理者に連絡してください。
====== 回転行列と四元数の計算速度比較 ======

===== 目的 =====

ベクトル，回転行列，座標変換行列といった三次元幾何演算をpythonのリストをベースに[[./geo.py|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|matrix_vs_quaternion.ipynb]]となっている．

===== モジュールの読み込み =====

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

<code python>
from geo import *
import time
import pandas as pd
</code>
===== データの生成 =====

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

<code python>
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])
</code>
この内容はこうなる．

<code python>
print('v1 =', v1)
print('R1 =', R1)
print('R2 =', R2)
print('q1 =', q1)
print('q2 =', q2)
print('qv1 =', qv1)
</code>
<code>
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])
</code>
===== 時間計測の関数と結果保存データ =====

<code python>
def test(n,fn):
    i=0
    start=time.time()
    while i< n :
        fn()
        i += 1
    end = time.time()
    rslt=end-start
    return rslt
</code>
<code python>
data = []
</code>
<code python>
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
</code>
===== ループ回数の決定 =====

<code python>
test(100, lambda : R1*R2)
</code>
<code>
0.0013973712921142578
</code>
<code python>
test(1000, lambda : R1*R2)
</code>
<code>
0.015026330947875977
</code>
<code python>
test(10000, lambda : R1*R2)
</code>
<code>
0.11950469017028809
</code>
<code python>
test(100000, lambda : R1*R2)
</code>
<code>
0.7195339202880859
</code>
<code python>
test(1000000, lambda : R1*R2)
</code>
<code>
9.323349952697754
</code>
<code python>
test(10000000, lambda : R1*R2)
</code>
<code>
78.96706867218018
</code>
百万回ぐらいでループ前後のオーバーヘッドの影響が少なくなってきている． まだ多少影響はあるが，一千万，一億は時間がかかるし， どうせループ内の処理の影響は消せないので百万回に決定する．

<code python>
N=1000000
</code>
===== オブジェクトの生成 =====

<code python>
m_time = test(N, lambda : MATRIX())
print(m_time)
</code>
<code>
1.5926079750061035
</code>
<code python>
q_time = test(N, lambda : QUATERNION())
print(q_time)
</code>
<code>
1.4729571342468262
</code>
<code python>
data.append(judge('オブジェクトの生成', m_time, q_time))
</code>
===== 回転の合成　QUATERNIONの勝ち =====

<code python>
R1*R2
</code>
<code>
m:[[0.8660254037844387, 0.0, 0.49999999999999994], [0.43301270189221924, 0.5000000000000001, -0.75], [-0.25, 0.8660254037844386, 0.43301270189221946]]
</code>
<code python>
(R1*R2).quaternion()
</code>
<code>
q:(0.8365163037378079, v:[0.4829629131445341, 0.22414386804201336, 0.1294095225512603])
</code>
<code python>
q1*q2
</code>
<code>
q:(0.8365163037378079, v:[0.4829629131445341, 0.2241438680420134, 0.12940952255126034])
</code>
<code python>
m_time = test(N, lambda : R1*R2)
print(m_time)
</code>
<code>
8.527657508850098
</code>
<code python>
q_time = test(N, lambda : q1*q2)
print(q_time)
</code>
<code>
7.0049989223480225
</code>
<code python>
data.append(judge('回転の合成', m_time, q_time))
</code>
===== ベクトルの回転　MATRIXの圧勝 =====

<code python>
R1*v1
</code>
<code>
v:[1.0, -1.5980762113533158, 3.2320508075688776]
</code>
<code python>
q1*qv1*q1.conjugate()
</code>
<code>
q:(0.0, v:[1.0, -1.5980762113533153, 3.2320508075688776])
</code>
<code python>
m_time = test(N, lambda : R1*v1)
print(m_time)
</code>
<code>
3.5639126300811768
</code>
<code python>
q1c=q1.conjugate()
q_time = test(N, lambda : q1*qv1*q1c)
print(q_time)
</code>
<code>
13.009745597839355
</code>
<code python>
data.append(judge('ベクトルの回転', m_time, q_time))
</code>
===== 結論 =====

まとめの表

<code python>
df = pd.DataFrame(data, columns=["項目", "MATRIX", "QUATERNION", "結果"])
df
</code>
<HTML>
<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

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

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
<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>
</table>
</div>
</HTML>
結論から言うと，ほぼ予想通りである．

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

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

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

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

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

<code python>

</code>