Trouble with 'solve' over ZZ

[lkastner]

I am using 'solve' to find the integral solution of a system of linear equations. I know that there is a unique solution. It works sometimes, other times it does not, and then I promote everything to QQ and afterwards lift back to ZZ:

im = matrix {{1, 1, 1, 1, 1, 1, 1, 1}, {1, -1, 0, 0, 0, 0, 0, 0}, {0, 0, 1, -1, 0, 0, 0, 0}, {0, 0, 0, 0, 1, -1, 0, 0}, {0, 0, 0, 0, 0, 0, 1, -1}}
v = matrix {{0, 2, 2, 0, 0, 2, 0, 2}}
solve(transpose im, transpose v)
solve(promote(transpose im, QQ), promote(transpose v, QQ))

If this is intended behaviour, maybe the error message could be changed to indicate that one should use another ring?

[DanGrayson]

The change in this commit

commit b8701f5ee83cad81534168e10f92977efd4dd994
Author: Mike Stillman <mikestillman1@gmail.com>
Date:   Wed Dec 19 12:25:10 2012 -0500

    started to connect fast linear algebra routines to top level functions.  Also: changed solve of MutableMatrices to return a MutableMatrix

diff --git a/M2/Macaulay2/m2/mutablemat.m2 b/M2/Macaulay2/m2/mutablemat.m2
index 08a6b0e..02693b7 100644
--- a/M2/Macaulay2/m2/mutablemat.m2
+++ b/M2/Macaulay2/m2/mutablemat.m2
@@ -180,9 +184,9 @@ solve(MutableMatrix,MutableMatrix) := opts -> (A,b) -> (
      if opts.ClosestFit
      then rawLeastSquares(raw A,raw b,raw x,opts.MaximalRank)
      else rawSolve(raw A,raw b,raw x);
-     matrix x)
+     x)
 solve(Matrix,Matrix) := opts -> (A,b) -> (
-     solve(mutableMatrix(A,Dense=>true),
+     matrix solve(mutableMatrix(A,Dense=>true),
                   mutableMatrix(b,Dense=>true),
          opts))

interacts with an undocumented feature of "solve": it returns "null" if there is no solution.

[lkastner]

Ok, but there is another thing here: There is a solution over ZZ, so why is null appearing?

i2 : im = matrix {{1, 1, 1, 1, 1, 1, 1, 1}, {1, -1, 0, 0, 0, 0, 0, 0}, {0, 0, 1, -1, 0, 0, 0, 0}, {0, 0, 0, 0, 1, -1, 0, 0}, {0, 0, 0, 0, 0, 0, 1, -1}}

o2 = | 1 1  1 1  1 1  1 1  |
     | 1 -1 0 0  0 0  0 0  |
     | 0 0  1 -1 0 0  0 0  |
     | 0 0  0 0  1 -1 0 0  |
     | 0 0  0 0  0 0  1 -1 |

              5        8
o2 : Matrix ZZ  <--- ZZ

i3 : v = matrix {{0, 2, 2, 0, 0, 2, 0, 2}}

o3 = | 0 2 2 0 0 2 0 2 |

              1        8
o3 : Matrix ZZ  <--- ZZ

i4 : solve(transpose im, transpose v)
stdio:4:1:(3): error: no method found for applying matrix to:
     argument   :  null (of class Nothing)

i5 : solve(promote(transpose im, QQ), promote(transpose v, QQ))

o5 = | 1  |
     | -1 |
     | 1  |
     | -1 |
     | -1 |

              5        1
o5 : Matrix QQ  <--- QQ

i6 :

[DanGrayson]

@mikestillman ?

[mikestillman]

Yes? I was emailing Lars and suggested he post this issue.

On Jun 3, 2016, at 2:20 PM, Daniel R. Grayson notifications@github.com wrote:

@mikestillman https://github.com/mikestillman ?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/Macaulay2/M2/issues/479#issuecomment-223654427, or mute the thread https://github.com/notifications/unsubscribe/AAGPR1ccg3hGQIwmD3KQyxB6woc7yTtiks5qIHBXgaJpZM4Itiq4.

[DanGrayson]

Okay.

[mikestillman]

I think I should not have switched the linear algebra for matrices over ZZ (that was just a mistake, unfortunately, there were no tests of this behavior), possibly over QQ too. I need to review that.

On Jun 3, 2016, at 2:26 PM, Daniel R. Grayson notifications@github.com wrote:

Okay.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/Macaulay2/M2/issues/479#issuecomment-223656014, or mute the thread https://github.com/notifications/unsubscribe/AAGPR_pHgaFj2-m9fJCn8l-5YzU3zCdLks5qIHHAgaJpZM4Itiq4.

[mikestillman]

The 'solve' function, as implemented, is not designed to work over ZZ. There are two workarounds:

1. promote all matrices to QQ, solve there (as done above)
2. Use //:

(transpose v) // (transpose im)

(edited: I had said (transpose im) // (transpose v), but I had the order switched around.)

I think I will change the 'solve' function for matrices over ZZ, to use this latter method. However, this could be slower (not for 'reasonable size' examples, but for larger ones). So really, we should make sure linear algebra over ZZ is implemented directly.

[mikestillman]

I am changing the code so that many linear algebra operations which involve LU decomposition, are not defined over ZZ. The code was incorrectly calling routines meant only for fields.

The problems are not over yet, as I need to revisit other rings too. For example, LUdecomposition, when called with a matrix over a polynomial ring, also gives the wrong answer.

Really, I need to revisit the templated functions in e/mat-linalg.hpp, as this is hard to understand, and the defaults are poor (i.e. incorrect, in the case when the ring is not a field).

[tom111]

I think this issue is still not solved. When faced with an inconsisten system over ZZ, solve gives non-sensical results:

i9 : A = matrix {{0_ZZ,0},{0,1}}

o9 = | 0 0 |
     | 0 1 |

              2        2
o9 : Matrix ZZ  <--- ZZ

i10 : b = matrix {{1_ZZ},{1}}

o10 = | 1 |
      | 1 |

               2        1
o10 : Matrix ZZ  <--- ZZ

i11 : solve (A,b)

o11 = | 0 |
      | 1 |

               2        1
o11 : Matrix ZZ  <--- ZZ

[jchen419]

The same error message

error: no method found for applying matrix to: argument : null (of class Nothing)

appears in using solve over RR as well: see this Google groups post. This seems to occur when M2 cannot find a solution (even if a solution exists). It would be good to have a slightly more informative warning/error message, e.g. explicitly stating that no solution exists if this is the case. Over RR, using solve(A, b, ClosestFit => true) seems to avoid triggering an error (although a non-solution may be returned).