Hi,
Recently I have migrated a simple parallel orthogonalization function from numpy to pytorch. However, the torch.eig() function seems very unstable (or the precision is very low).
The core of the code is the following:
# Step1: Initialization:
P_cuda = torch.randn([300,40000]).type(dtype)

# Step2: orthogonalization:
L = torch.mm(P_cuda,P_cuda.t())
w,v = torch.eig(L,eigenvectors=True)
U = torch.mm(torch.mm(v,torch.diag(torch.sqrt(1./w[:,0]))),v.t())
P_cuda = torch.mm(U,P_cuda)

# Step3: Now the rows of P_cuda should be orthogonal and normalized that L2 norm is 1. So the output of the following code should be very close to 300.
normSquareSum = torch.pow(P_cuda,2).sum()
print(normSquareSum)

The first iteration of step2 and step seems successful. However, when I repeat them, I get a bunch of different numbers range from 300 to 370. Is this due a low precision of torch.eig()? If this is the case, how can I improve the precision? This  result makes the algorithm almost unusable since if I loop these 2 steps a little while, it will give me a nan.
The above code has a numpy implementation, which is a direct translation. The numpy version on the other hand is very stable and gives constantly 300. It is attached here too for comparison:
def orthogonalization_parallel(P):
    #This method orthogonalizes and normalizes the rows of P and the original matrix will not be changed.
    #P -> P_ such that P_*P_^T = I
    w, v = la.eig(np.dot(P,P.T))
    U = np.dot(np.dot(v,np.diag(np.sqrt(1./w))),v.T)
    return np.dot(U,P)

My pytorch version is ‘0.2.0_4’.
Thanks a lot for your help in advance!