Problem

I want to compare two rotation matrices $R_A$ and $R_B$ both representing the orientation of the same point cloud in space, but computed from different methods. The idea is to have an estimation of the error between those two matrices.

Method

My idea was to do it as follows:

1. Compute the rotation $R_{AB}$ between $R_A$ and $R_B$ as $R_{AB} = R_A^TR_B$
2. Compute the axis-angle ($\omega$, $\theta$) representation of $R_{AB}$ using the following formula: $$Tr(R_A) = 1 + 2cos(\theta)$$
3. Use the angle $\theta$ as the rotation error.

In Python, I do:

r_oa = import(R_A) // See R_A below, in "Data"
r_ob = import(R_B) // See R_B below, in "Data"

r_oa_t = np.transpose(r_oa)
r_ab = r_oa_t * r_ob



That seems quite straightforward to me, but I get $\theta = 23.86\unicode{xb0}$, which seems really unrealistic. This intuition is confirmed by the use of a software comparing the matrices, as described below in "Verification".

Am I doing something wrong there? Is it just not the right way to compare my matrices?

Verification

In a processing software, I have the possibility to "compare" two matrices (each being a 4x4 matrix defining a position and an orientation: [R | t]).

The documentation does not explain exactly how it works, but the output of this comparison is:

• 4 values for the rotation [deg]: Omega, Phi, Kappa, Total
• 4 values for the translation: $\Delta$x, $\Delta$y, $\Delta$z, Total.

I assumed that the "Total" angle for the rotation is my $\theta$ (as described above). For the same matrices $R_A$ and $R_B$ as above, the software outputs: Total = 0.036477551. That seems much more realistic, but then I don't know why it is not what I get.

Just to check my computation of $\theta$, I tried to compute it on $R_A$ and $R_B$ (instead of $R_{AB}$) and to compare that to the output of the software when running its comparison between $R_A$ and $I$, and also between $R_B$ and $I$. Even though I don't get the exact same $\theta$ as the software, my result is really close (to the 4th decimal).

Assuming that the software loses precision somewhere, I believe that my computation of $\theta$ is correct. And my computation of $R_{AB}$ is quite straightforward, so I don't understand why I don't find the same $\theta$ for $R_{AB}$.

Data

My matrices are:

$R_A = \begin{bmatrix} -0.956395958000000 & 0.292073230000000 & 0.000014880000000 \\ -0.292073218000000 & -0.956395931000000 & 0.000242173000000 \\ 0.000084963000000 & 0.000227268000000 & 0.999999971000000\end{bmatrix}$

$R_B = \begin{bmatrix} -0.956227882000000 & 0.292623030000000 & -0.000013768000000 \\ -0.292623029000000 & -0.956227882000000 & -0.000029806000000 \\ -0.000021887000000 & -0.000024473000000 & 0.999999999000000 \end{bmatrix}$

The software outputs, for the comparison between $R_A$ and $R_B$: "Total = 0.036477551". Following my method, I get $\theta = 23.86$, which is completely different.

Related questions

• Don’t you have some documentation for this software that explains the values it produces? – amd Jan 25 '17 at 19:29
• There is some documentation, but not for that. At least not that I found =(. But one point is that I don't think it is realistic to have an angle of 23 degrees between both matrices. 0.03 seems much more reasonable! – JonasVautherin Jan 25 '17 at 23:50
• Please show us how you get to this angle of $23.86°$. I don't think this is right. – Yves Daoust Jan 26 '17 at 11:32
• In python, I do: np.rad2deg(np.arccos((np.trace(r_ab) - 1) / 2)), with r_ab = r_oa_t * r_ob, r_oa_t = np.transpose(r_oa), r_oa = $R_A$ and r_ob = $R_B$ :-/. – JonasVautherin Jan 26 '17 at 11:47
• But the way, it seems that * means array multiplication, not matrix multiplication in NumPy. – user1551 Jan 26 '17 at 12:38

Apparently you have mistaken array multiplication for matrix multiplication. In Octave, acos((trace(A'*B)-1)/2)*360/2/pi gives 0.0365 degree, but acos((trace(A'.*B)-1)/2)*360/2/pi (note it's .* here instead of *) gives 23.86 degrees.

In NumPy, * means array multiplication, not matrix multiplication. Your Python code seems to be incorrect.

• And I was looking for an error in the maths =/. – JonasVautherin Jan 26 '17 at 12:54

A direct way to measure the angle between matrices is to view them as vectors in $\mathbb{R}^{n^2}$ and compute the cosine between these vectors as usual. Notation:
$\qquad x, y$ : 1d vectors
$\qquad A, B$ : $n \times n$ matrices
$\qquad A_{flat}, B_{flat}$ : the corresponding vectors in $\mathbb{R}^{n^2}$

$\qquad dot( x, y ) = \sum x_i y_i$
$\qquad \|x\| = \sqrt \, dot( x, x )$
$\qquad cos( x, y ) = {dot( x, y ) \over {\|x\| \, \|y\|}}$

$\qquad Fdot( x, y ) = dot( \, A_{flat}, \, B_{flat} \, ) = \sum \sum A_{ij} B_{ij}$
$\qquad \|A\|_F = \sqrt \, Fdot( A, A ) \ \$ -- Frobenius norm
$\qquad Fcos( A, B ) = {Fdot( \, A, \, B\, ) \over {\|A\|_F \, \|B\|_F}}$

Rotation matrices have $\|R\|_F = \sqrt n$ (check $\|I\|_F$), so $Fcos( R, S ) = {1 \over n} Fdot( R, S )$. ($trace( A^T B )$ is also $\sum \sum A_{ij} B_{ij}$, but your formula is missing the factor ${1 \over n}$.)

In python numpy, the same block of numbers can be viewed as a matrix (actually ndarray), or as a 1d vector:

def cos( A, B ):
""" comment cos between vectors or matrices """
Aflat = A.reshape(-1)  # views
Bflat = B.reshape(-1)
return (np.dot( Aflat, Bflat )
/ max( norm(Aflat) * norm(Bflat), 1e-10 ))