M2 icon indicating copy to clipboard operation
M2 copied to clipboard

Mutable Matrix Null Space Concern

Open davikrehalt opened this issue 2 years ago • 5 comments

If we calculate the nullspace of the matrix a=mutableMatrix{{0.0,1.0,2.0,3.0,0.0,0.0,0.0,0.0},{4.0,3.0,2.0,1.0,0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0,1.0,2.0,3.0,4.0},{0.0,0.0,0.0,0.0,3.0,2.0,1.0,0.0}} using 'nullSpace(a)' it outputs 'i2 : nullSpace(a)

o2 = | -1 -2 0 0 | | 2 3 0 0 | | -1 0 0 0 | | 0 -1 0 0 | | 0 0 2 4 | | 0 0 .5 0 | | 0 0 -1 0 | | 0 0 0 -1 |'

which has the following property i3 : a*nullSpace(a)

o3 = | 0 0 0 0 | | 0 0 0 0 | | 0 0 0 0 | | 0 0 6 12 |

o3 : MutableMatrix

which does not seem correct.

davikrehalt avatar Oct 01 '23 06:10 davikrehalt

Very strange! In the meantime, you can use ker(Matrix), which seems to work much better:

i2 : mutableMatrix gens ker matrix a
-- warning: experimental computation over inexact field begun
--          results not reliable (one warning given per session)

o2 = | 1  2  0  0  |
     | -2 -3 0  0  |
     | 1  0  0  0  |
     | 0  1  0  0  |
     | 0  0  1  2  |
     | 0  0  -2 -3 |
     | 0  0  1  0  |
     | 0  0  0  1  |

o2 : MutableMatrix

i3 : a * oo

o3 = | 0 0 0           0            |
     | 0 0 0           0            |
     | 0 0 4.44089e-16 4.44089e-16  |
     | 0 0 0           -8.88178e-16 |

o3 : MutableMatrix

d-torrance avatar Oct 01 '23 11:10 d-torrance

Yes that works! However, I think this should still be labeled as a bug!

davikrehalt avatar Oct 05 '23 02:10 davikrehalt

if you let me know where the bug it, I can try to take a look and try to fix it!

davikrehalt avatar Oct 08 '23 18:10 davikrehalt

I think it's somewhere in the engine:

i1 : code (nullSpace, MutableMatrix)

o1 = -- code for method: nullSpace(MutableMatrix)
     /usr/share/Macaulay2/Core/mutablemat.m2:323:32-323:70: --source code:
     nullSpace(MutableMatrix) := (M) -> map(ring M, rawLinAlgNullSpace raw M)

Any time something starts with raw, it means that it's likely written in C++ in the e directory of the repository. In this case:

https://github.com/Macaulay2/M2/blob/379c9634000dd87221dc353e1aaaca413f75c81c/M2/Macaulay2/e/interface/mutable-matrix.cpp#L824-L834

d-torrance avatar Oct 11 '23 01:10 d-torrance

Hmm, is https://github.com/Macaulay2/M2/blob/379c9634000dd87221dc353e1aaaca413f75c81c/M2/Macaulay2/e/dmat-ffpack.cpp#L85C4-L126C6 getting called by the code? Edit: Guess not bc it's for finite fields? I still don't get what's the underlying code that's getting called lol. like where is the nullSpace() in your code actually implemented?

davikrehalt avatar Oct 14 '23 19:10 davikrehalt