sagemath / sage

Main repository of SageMath
https://www.sagemath.org
Other
1.28k stars 438 forks source link

Incorrect matrix inverse over SR #38336

Open galashin opened 1 month ago

galashin commented 1 month ago

Steps To Reproduce

M = Matrix(SR,[[-0.193469699433335, -0.183934430242923, -9.23026395957241, -234.680519425251, -4.56220425542616, -1.00000000000000, 0],
[-4.41006516415701, -4.04624755225825, -34.3402432065932, -712.818780249295, -9.03026688522083, 0, 1.00000000000000],
[-2.65271365792984, -2.42946892291712, -15.3627946531470, -289.367653048543, -2.62301553224534, 0, 0],
[-2.42930274338041, -2.22351202430932, -12.6512341012439, -230.654739849570, -1.95615689827990, 0, 0],
[-5.02850517118872, -0.940972164125432, -19.9103986565546, -236.148509713831, -8.90077623301954, 7.46410161513775, -6.46410161513775],
[-0.0799721747292845, -8.59828369225104, 3.94544109797720, -103.092869883371, 13.5928135282255, -12.9282032302755, 10.1961524227066],
[-5.35162679183391, 2.77378312755603, -12.0987003087950, -73.6532885076776, -7.76222579473772, 6.46410161513775, -4.73205080756887]])
print(M.inverse() * M)

Expected Behavior

Identity matrix expected

Actual Behavior

Does not output anything close to the identity matrix. Removing SR, from the matrix definition produces the identity matrix.

Additional Information

No response

Environment

- **OS**: Mac OS Ventura 13.5.1
- **Sage Version**: 10.0

Checklist

dvadym commented 1 month ago

Note: that’s my first encounter with Sage codebase, maybe I’ll miss something trivial or use incorrect terminology

I’ve debugged why it happens:

  1. Inverse is computed with attaching identity matrix and using echolize method ([code]())
  2. For matrix on doubles the scaled_partial_pivoting method is used (code).
  3. For matrix on SR classical method is used (code)
  4. The reason, by classical is used for SR is implementation of echelonize for Matrix_symbolic_dense (code). Which is probably correct, since there is no absolute value in SR.
  5. Classical method fails here, on 4th row the pivot value is too small, as a result there are too large numerical errors.

So in general it seems that everything worked as expected, it seems there is no way in general to use scaled_partial_pivoting for SR (or am I missing something?). It’s just that classical method can’t work here. Maybe some other method?

One of the possible approaches here is to use Maxima for inverting, AFAIU Matrix_symbolic_dense heavily uses Maxima anyway. And Maxima handles this inversion.

I’ve checked adding to Matrix_symbolic_dense class method:

  def __invert__(self):
        return self._maxima_().invert()._sage_()

fixes this problem. Though, this probably requires thorough testing (including performance) for comparing Maxima inversion vs the existing algorithm.