so3.log gives an incorrect result for some matrix

The log of a rotation matrix obtained from so3.log for the following matrix is incorrect:

import torch

import ptlk

mat = torch.Tensor([
    [-1.00000396e+00, -9.55433245e-07,  1.04267154e-06],
    [ 1.04267254e-06, -9.99052394e-01,  4.36201482e-02],
    [ 9.55432245e-07,  4.36191482e-02,  9.99051394e-01]])

print(ptlk.so3.log(mat))
# prints:
# tensor([0., 0., 0.])

print(ptlk.so3.exp(ptlk.so3.log(mat)))
# prints:
# tensor([[1., 0., 0.],
#         [0., 1., 0.],
#         [0., 0., 1.]])

print(torch.allclose(mat, ptlk.so3.exp(ptlk.so3.log(mat)), atol=1e-2))
# prints:
# False

I tried to change eps in so3.log function, but it does not seem to help.

Ok, I think that I understand what happened. Essentially, this is a numerical precision error.
I will close the issue because I did not find any better way to calculate the coefficients in this case.

Here is a detailed explanation, in case one might find it useful in the future.

The generators' coefficients are calculated using the following:

Usually, the second term (red) in (15) is nonzero or far enough from zero. In such a case, we can use:

(Where R=exp(\omega_x) ).
However, when sin(\theta)/(\theta) is zero or very small the above does not work properly. It happens when approximately \theta=(+/-)\pi, as in the case above. In this case, we have to use the quadratic order in \theta (blue in (15), (14) is still valid):

The numerical problem occurs because \omega_1^2=0 but \omega_1\omega_2 is not zero. Unfortunately, it looks like this method is also not stable in this case.

The part in the code of so3 which performs the calculation:

PointNetLK/ptlk/so3.py

Lines 99 to 116 in 8b851c1

    
           if idx0.any(): 
        
               # t[idx0] == math.pi 
        
               t2 = t[idx0] ** 2 
        
               A = (R[idx0] + torch.eye(3).type_as(R).unsqueeze(0)) * t2.view(-1, 1, 1) / 2 
        
               aw1 = torch.sqrt(A[:, 0, 0]) 
        
               aw2 = torch.sqrt(A[:, 1, 1]) 
        
               aw3 = torch.sqrt(A[:, 2, 2]) 
        
               sgn_3 = torch.sign(A[:, 0, 2]) 
        
               sgn_3[sgn_3 == 0] = 1 
        
               sgn_23 = torch.sign(A[:, 1, 2]) 
        
               sgn_23[sgn_23 == 0] = 1 
        
               sgn_2 = sgn_23 * sgn_3 
        
               w1 = aw1 
        
               w2 = aw2 * sgn_2 
        
               w3 = aw3 * sgn_3 
        
               w = torch.stack((w1, w2, w3), dim=-1) 
        
               W = mat(w) 
        
               X[idx0] = W

(Especially L102.)

The above is based on:
http://www.ethaneade.org/lie.pdf

	if idx0.any():
	# t[idx0] == math.pi
	t2 = t[idx0] ** 2
	A = (R[idx0] + torch.eye(3).type_as(R).unsqueeze(0)) * t2.view(-1, 1, 1) / 2
	aw1 = torch.sqrt(A[:, 0, 0])
	aw2 = torch.sqrt(A[:, 1, 1])
	aw3 = torch.sqrt(A[:, 2, 2])
	sgn_3 = torch.sign(A[:, 0, 2])
	sgn_3[sgn_3 == 0] = 1
	sgn_23 = torch.sign(A[:, 1, 2])
	sgn_23[sgn_23 == 0] = 1
	sgn_2 = sgn_23 * sgn_3
	w1 = aw1
	w2 = aw2 * sgn_2
	w3 = aw3 * sgn_3
	w = torch.stack((w1, w2, w3), dim=-1)
	W = mat(w)
	X[idx0] = W