facelessuser
commented
1 year ago So, once we added in all the checking it's probably comparable in speed, but we can use it for inv and it doesn't require all the complexity to ensure it can invert. It's now comparable to numpy and scipy's implementation.
With that said, both numpy and scipy, due to the perils of floating point math, sometimes try to invert singular matrices, and this is because of floating point error causing the determinants to sometimes be very close to zero, but not zero. Ours bears the same behavior.
One thing to note, it seems there must be a bug in scipy's ul. I thin they may not have partial pivoting done quite correctly, or maybe it is LAPACK, consider the following. When decomposing, PA = LU, but this doesn't always happen with Scipy. Notice that the order is all wrong.
>>> import numpy as np
>>> import scipy.linalg as la
>>> from coloraide import algebra as alg
>>> m = [[-22.966983131482493, 0.0, 0.0],
... [-21.359566335117776, 0.0, 0.0],
... [92.44620809845696, -55.206609400991134, -97.35747258042645]]
>>> p, l, u = la.lu(m)
>>> np.dot(p, m)
array([[-21.35956634, 0. , 0. ],
[ 92.4462081 , -55.2066094 , -97.35747258],
[-22.96698313, 0. , 0. ]])
>>> np.dot(l, u)
array([[ 9.24462081e+01, -5.52066094e+01, -9.73574726e+01],
[-2.29669831e+01, 0.00000000e+00, 0.00000000e+00],
[-2.13595663e+01, -1.01741598e-15, 9.14592066e-17]])Now, consider ours:
>>> import numpy as np
>>> from coloraide import algebra as alg
>>> m = [[-22.966983131482493, 0.0, 0.0],
... [-21.359566335117776, 0.0, 0.0],
... [92.44620809845696, -55.206609400991134, -97.35747258042645]]
>>> p, l, u = alg.lu(m)
>>> alg.pprint(alg.dot(p, m))
[[92.44620809845696, -55.206609400991134, -97.35747258042645],
[-22.966983131482493, 0.0, 0.0],
[-21.359566335117776, 0.0, 0.0]]
>>> alg.pprint(alg.dot(l, u))
[[92.44620809845696, -55.206609400991134, -97.35747258042645],
[-22.966983131482493, 0.0, 0.0],
[-21.359566335117776, 0.0, 0.0]]I honestly didn't expect that ours would be more reliable in this sense. I haven't looked at their code, but I know they use LAPACK, so if there is a bug, I imagine it is there, but it seems as if that library is in Fortran, and I don't know Fortran.
Now, for these kinds of singular matrices, we don't 100% match Scipy when determinants are zero, or have small error being close to zero, I imagine there are some small implementation details that affect the floating point error in our favor or in their favor, but regardless of when we fail or not fail for a singular matrix due to the determinant not being exactly zero, ours seems to always resolve the LU properly while Scipy can be hit or miss.
I have my theories as to why theirs can fail, but I'd have to dig into their code base and see exactly how they implement their pivoting.
Still needs tests