Problem with computing determinant of 3x3 matrix

[enokjun]

Describe the bug
The computation of the determinant of 3x3 matrix in taichi produces incorrect results

To Reproduce

import numpy as np
import taichi as ti

ti.init(arch=ti.gpu)

# numpy method for comparison
A = np.array([[2.52825e+04, 2.38500e+04, 4.77000e+02],
       [2.38500e+04, 2.25015e+04, 4.50000e+02],
       [4.77000e+02, 4.50000e+02, 9.00000e+00]])

B = np.array([[8.10150e+03, 1.35000e+04, 2.70000e+02],
       [1.35000e+04, 2.25015e+04, 4.50000e+02],
       [2.70000e+02, 4.50000e+02, 9.00000e+00]])

# create taichi 3x3 matrix based on numpy array
temp3n3_matrix = ti.types.matrix(3, 3, ti.f32)

A_ti = temp3n3_matrix(0)
A_ti[0,0] = A[0,0]
A_ti[0,1] = A[0,1]
A_ti[0,2] = A[0,2]
A_ti[1,0] = A[1,0]
A_ti[1,1] = A[1,1]
A_ti[1,2] = A[1,2]
A_ti[2,0] = A[2,0]
A_ti[2,1] = A[2,1]
A_ti[2,2] = A[2,2]

B_ti = temp3n3_matrix(0)
B_ti[0,0] = B[0,0]
B_ti[0,1] = B[0,1]
B_ti[0,2] = B[0,2]
B_ti[1,0] = B[1,0]
B_ti[1,1] = B[1,1]
B_ti[1,2] = B[1,2]
B_ti[2,0] = B[2,0]
B_ti[2,1] = B[2,1]
B_ti[2,2] = B[2,2]

# use taichi math module
@ti.kernel
def taichi_det_v1(M: ti.types.matrix(3,3,ti.f32)) -> ti.f32:
    return ti.math.determinant(M)

# directly compute 3x3 matrix determinant from formula
@ti.kernel
def taichi_det_v2(M: ti.types.matrix(3,3,ti.f32)) -> ti.f32:
    return M[0,0]*(M[1,1]*M[2,2] - M[1,2]*M[2,1]) - M[0,1]*(M[1,0]*M[2,2] - M[1,2]*M[2,0]) + M[0,2]*(M[1,0]*M[2,1] - M[1,1]*M[2,0])

# use svd to compute determinants
@ti.kernel
def taichi_det_v3(M: ti.types.matrix(3,3,ti.f32)) -> ti.f32:
    U,s,V = ti.svd(M, ti.f32) 
    return s[0,0]*s[1,1]*s[2,2]

# fails to calculate correct for matrix A
print(f'\nnumpy - determinant of 3x3 matrix A: {np.linalg.det(A)}')
print(f'taichi.math.determinant - determinant of 3x3 matrix A: {taichi_det_v1(A_ti)}')
print(f'taichi - formula - determinant of 3x3 matrix A: {taichi_det_v2(A_ti)}')
print(f'taichi svd - determinant of 3x3 matrix A: {taichi_det_v3(A_ti)}')

# but seems to be match for matrix B
print(f'\nnumpy - determinant of 3x3 matrix B: {np.linalg.det(B)}')
print(f'taichi.math.determinant - determinant of 3x3 matrix B: {taichi_det_v1(B_ti)}')
print(f'taichi - formula - determinant of 3x3 matrix B: {taichi_det_v2(B_ti)}')
print(f'taichi svd - determinant of 3x3 matrix B: {taichi_det_v3(B_ti)}')

Log/Screenshots

[Taichi] version 1.7.0, llvm 15.0.1, commit 2fd24490, win, python 3.10.11
[W 12/01/23 11:33:38.121 21492] [cuda_driver.cpp:taichi::lang::CUDADriverBase::load_lib@36] nvcuda.dll lib not found.
[Taichi] Starting on arch=vulkan

numpy - determinant of 3x3 matrix A: 20.24999999995664
taichi.math.determinant - determinant of 3x3 matrix A: -218.25
taichi - formula - determinant of 3x3 matrix A: -218.25
taichi svd - determinant of 3x3 matrix A: 20.257919311523438

numpy - determinant of 3x3 matrix B: 20.250000000027377
taichi.math.determinant - determinant of 3x3 matrix B: 20.25
taichi - formula - determinant of 3x3 matrix B: 20.25
taichi svd - determinant of 3x3 matrix B: 20.276962280273438

Additional comments

1. depending on the matrix, the accuracy of the taichi.math.determinant would change.
2. the svd would produce a close-enough approximation to the determinant. However, this error compounds as performing the inversion of the 3x3 matrix. While solving the following least-square linear multi-regression, the error is very apparent.
   For example,

A = np.array([[29.5, 49.5,  1. ],
       [30. , 49.5,  1. ],
       [30.5, 49.5,  1. ],
       [29.5, 50. ,  1. ],
       [30. , 50. ,  1. ],
       [30.5, 50. ,  1. ],
       [29.5, 50.5,  1. ],
       [30. , 50.5,  1. ],
       [30.5, 50.5,  1. ]])
y = np.array([[81.307915],
       [81.307915],
       [81.307915],
       [81.01924 ],
       [81.01924 ],
       [81.01924 ],
       [80.73056 ],
       [80.73056 ],
       [80.73056 ]])
alpha = np.dot((np.dot(np.linalg.inv(np.dot(A.T,A)),A.T)),y)
# results = [ 0.      -0.57735   109.88696 ]

# performing equivalent matrix operation in taichi with
# local_xy_matrix is taichi matrix version of numpy.array A
# local_z_vector is taichi vector version of numpy.array y

alpha_ti = (ti.math.inverse(local_xy_matrix.transpose()@local_xy_matrix)@local_xy_matrix.transpose())@local_z_vector
# results = [ 0.86705124      -0.57683617     219.1809  ]
# error seems to arise after inversion of 3x3 matrix ( = local_xy_matrix.transpose()@local_xy_matrix )

3. When testing if the matrix inversion itself was the issue, the normalized inversion seems to be accurate

# show inversion is not the problem, the determinant computation is
@ti.kernel
def taichi_inv_norm(M: ti.types.matrix(3,3,ti.f32)) -> ti.types.matrix(3,3,ti.f32):
    return ti.math.inverse(M)*ti.math.determinant(M)

print(np.linalg.inv(A)*np.linalg.det(A))
print(taichi_inv_norm(A_ti))

print()
print(np.linalg.inv(B)*np.linalg.det(B))
print(taichi_inv_norm(B_ti))

######### results show #########
# for matrix A
numpy version
[[ 1.35000000e+01  9.79780096e-12 -7.15500000e+02]
 [-1.73176166e-11  1.35000000e+01 -6.75000000e+02]
 [-7.15500000e+02 -6.75000000e+02  7.16737500e+04]]
taichi version
[[ 1.3500001e+01  0.0000000e+00 -7.1600000e+02]
 [ 0.0000000e+00  1.3500001e+01 -6.7500006e+02]
 [-7.1600000e+02 -6.7500006e+02  7.1645750e+04]]

# for matrix B
numpy version
[[ 1.35000000e+01 -5.27235589e-13 -4.05000000e+02]
 [-8.32455030e-13  1.35000000e+01 -6.75000000e+02]
 [-4.05000000e+02 -6.75000000e+02  4.59022500e+04]]
taichi version
[[ 1.350000e+01  0.000000e+00 -4.050000e+02]
 [ 0.000000e+00  1.350000e+01 -6.750000e+02]
 [-4.050000e+02 -6.750000e+02  4.590225e+04]]

[enokjun]

The issue was due to using f32 instead of f64, which is a standard problem that causes issues in determinant computation for any framework