Speed up IntegerMatrix

[tscrim]

We explicitly declare things to be int type so Cython can generate better C code.

[jdemeyer]

Works for me:

sage: from sage.matroids.linear_matroid import LinearMatroid
sage: from sage.matroids.lean_matrix import IntegerMatrix
sage: M = LinearMatroid(matrix=IntegerMatrix(4,3,matrix([[0,0,0],[0,0,0],[0,1,1],[1,0,1]])))
sage: M.circuits()
Iterator over a system of subsets

[tscrim]

Hmm...strange. Let me try from a fresh session to see if I somehow unintentionally corrupted things.

[tscrim]

Okay, I had one little seemingly innocuous change: I added an explicit return type int on the corresponding get method of IntegerMatrix. Once I removed that, it worked.

[tscrim]

So what was going wrong is that there is a bint is_nonzero method which simply was return self.get(r, c), and Cython could somehow not automatically convert that to a boolean check. Would that count as a Cython bug?

If it is not, I might also recycle this ticket with some Cython tweaks to IntegerMatrix to try and squeeze a bit more speed out of it.

[tscrim]

Well, there is a bug somewhere (on vanilla 8.4.beta4):

sage: M = LinearMatroid(matrix=IntegerMatrix(4,4,matrix([[0,0,0,0],[1,1,1,1],[0,1,1,2],[0,0,1,0]])))
sage: list(M.bases())
[frozenset({0, 1, 2}), frozenset({0, 2, 3})]

sage: Mp = LinearMatroid(matrix=GenericMatrix(4,4,matrix([[0,0,0,0],[1,1,1,1],[0,1,1,2],[0,0,1,0]])))
sage: list(Mp.bases())
[frozenset({0, 1, 2}), frozenset({0, 2, 3}), frozenset({1, 2, 3})]

[jdemeyer]

Replying to tscrim:

So what was going wrong is that there is a bint is_nonzero method which simply was return self.get(r, c), and Cython could somehow not automatically convert that to a boolean check.

I'm not following here... what do you mean with "Cython could somehow not automatically convert that to a boolean check"?

[tscrim]

Yes, that is correct (well, an int to a bint to be specific).

[tscrim]

At least, once I changed that to an explicit check != 0, it fixed the problem.

[Stefan]

Replying to tscrim:

Well, there is a bug somewhere (on vanilla 8.4.beta4):

sage: M = LinearMatroid(matrix=IntegerMatrix(4,4,matrix([[0,0,0,0],[1,1,1,1],[0,1,1,2],[0,0,1,0]])))
sage: list(M.bases())
[frozenset({0, 1, 2}), frozenset({0, 2, 3})]

sage: Mp = LinearMatroid(matrix=GenericMatrix(4,4,matrix([[0,0,0,0],[1,1,1,1],[0,1,1,2],[0,0,1,0]])))
sage: list(Mp.bases())
[frozenset({0, 1, 2}), frozenset({0, 2, 3}), frozenset({1, 2, 3})]

Most (all?) LinearMatroid? methods assume that we are allowed to pivot in the matrix, and do so without leaving the ring. Maybe that's what went wrong?

[jdemeyer]

Replying to tscrim:

At least, once I changed that to an explicit check != 0, it fixed the problem.

It's still not clear which check you mean and what the problem really is...

[Rudi]

The method GenericMatrix?.is_nonzero seems to be called only by LeanMatrix?.gauss_jordan_reduce() and LeanMatrix?.nonzero_positions_in row().

LeanMatrix?.pivot() does not even use is_nonzero(), but ‘s=self.get_unsafe(i,y)’ followed by ‘if s and ..’ (line 303 of lean_matrix.pyx). Apparently the conversion to a bool implicit in ‘if s’ does work in that context. It may be more effcient too, since depending on the base ring evaluating s!=0 may involve casting 0 as a ring element. I’m quite sure that for finite fields ‘if s’ is about 5 times faster than ‘if (s!=0)’.

So it may also be a solution to rewrite LeanMatrix?.gauss_jordan_reduce() and LeanMatrix?.nonzero_positions_in row() in the more efficient way used in pivot(), completely avoiding the use of is_nonzero().

[tscrim]

Replying to jdemeyer:

Replying to tscrim:

At least, once I changed that to an explicit check != 0, it fixed the problem.

It's still not clear which check you mean and what the problem really is...

Sorry, I was trying to answer quickly while I was out.

So when I explicitly tell Cython that get should return an int, then I have to do this change

     cdef bint is_nonzero(self, long r, long c) except -2:   # Not a Sage matrix operation
-        return self.get(r, c)
+        return self.get(r, c) != 0

otherwise I get the original error message in the ticket description.

[tscrim]

Replying to Stefan:

Replying to tscrim:

Well, there is a bug somewhere (on vanilla 8.4.beta4):

sage: M = LinearMatroid(matrix=IntegerMatrix(4,4,matrix([[0,0,0,0],[1,1,1,1],[0,1,1,2],[0,0,1,0]])))
sage: list(M.bases())
[frozenset({0, 1, 2}), frozenset({0, 2, 3})]

sage: Mp = LinearMatroid(matrix=GenericMatrix(4,4,matrix([[0,0,0,0],[1,1,1,1],[0,1,1,2],[0,0,1,0]])))
sage: list(Mp.bases())
[frozenset({0, 1, 2}), frozenset({0, 2, 3}), frozenset({1, 2, 3})]

Most (all?) LinearMatroid? methods assume that we are allowed to pivot in the matrix, and do so without leaving the ring. Maybe that's what went wrong?

But in the latter case, the matrix is over ZZ, so I would not think that is the problem.

[tscrim]

Replying to Rudi:

The method GenericMatrix?.is_nonzero seems to be called only by LeanMatrix?.gauss_jordan_reduce() and LeanMatrix?.nonzero_positions_in row().

LeanMatrix?.pivot() does not even use is_nonzero(), but ‘s=self.get_unsafe(i,y)’ followed by ‘if s and ..’ (line 303 of lean_matrix.pyx). Apparently the conversion to a bool implicit in ‘if s’ does work in that context. It may be more effcient too, since depending on the base ring evaluating s!=0 may involve casting 0 as a ring element. I’m quite sure that for finite fields ‘if s’ is about 5 times faster than ‘if (s!=0)’.

Yes, as you surmised, the s != 0 is slower because invokes the coercion framework. So if s is the best thing to do for speed in general. It is just in this case that Cython does not seem to be converting an int (from an inline function) into a bint.

So it may also be a solution to rewrite LeanMatrix?.gauss_jordan_reduce() and LeanMatrix?.nonzero_positions_in row() in the more efficient way used in pivot(), completely avoiding the use of is_nonzero().

What I was hoping to do was take advantage of the known return types in the integer matrix to improve the Cython code. With the explicit casting, my computation goes from 13.6s to 3s.

src/sage/matroids/lean_matrix.pxd

@@ -102 +102 @@
 cdef class IntegerMatrix(LeanMatrix):
     cdef int* _entries
 
-    cdef inline get(self, long r, long c)   # Not a Sage matrix operation
+    cdef inline int get(self, long r, long c)   # Not a Sage matrix operation
     cdef inline void set(self, long r, long c, int x)   # Not a Sage matrix operation
 
     cdef inline long row_len(self, long i) except -1   # Not a Sage matrix operation

src/sage/matroids/lean_matrix.pyx

@@ -2844 +2844 @@
         """
         return "IntegerMatrix instance with " + str(self._nrows) + " rows and " + str(self._ncols) + " columns"
 
-    cdef inline get(self, long r, long c):   # Not a Sage matrix operation
+    cdef inline int get(self, long r, long c):   # Not a Sage matrix operation
         return self._entries[r * self._ncols + c]
 
     cdef inline void set(self, long r, long c, int x):   # Not a Sage matrix operation
@@ -2877 +2877 @@
         return 0
 
     cdef bint is_nonzero(self, long r, long c) except -2:   # Not a Sage matrix operation
-        return self.get(r, c)
+        return self.get(r, c) != 0
 
     cdef LeanMatrix copy(self):   # Deprecated Sage matrix operation
         cdef IntegerMatrix M = IntegerMatrix(self._nrows, self._ncols)
@@ -2982 +2982 @@
         ignored.
         """
         cdef long i
+        cdef int sval
         if s is None:
             for i from 0 <= i < self._ncols:
                 self.set(x, i, self.get(x, i) + self.get(y, i))
         else:
+            sval = int(s)
             for i from 0 <= i < self._ncols:
-                self.set(x, i, self.get(x, i) + s * self.get(y, i))
+                self.set(x, i, self.get(x, i) + sval * self.get(y, i))
         return 0
 
     cdef int swap_rows_c(self, long x, long y) except -1:
@@ -3010 +3012 @@
         compatibility, and is ignored.
         """
         cdef long i
+        cdef int sval = int(s)
         for i from 0 <= i < self._ncols:
-            self.set(x, i, s * self.get(x, i))
+            self.set(x, i, sval * self.get(x, i))
         return 0
 
     cdef int rescale_column_c(self, long y, s, bint start_row) except -1:
@@ -3020 +3023 @@
         compatibility, and is ignored.
         """
         cdef long j
+        cdef int sval = int(s)
         for j from 0 <= j < self._nrows:
-            self.set(j, y, self.get(j, y) * s)
+            self.set(j, y, self.get(j, y) * sval)
         return 0
 
     cdef int pivot(self, long x, long y) except -1:   # Not a Sage matrix operation

[Rudi]

Replying to tscrim:

Replying to Rudi: What I was hoping to do was take advantage of the known return types in the integer matrix to improve the Cython code. With the explicit casting, my computation goes from 13.6s to 3s.

That looks like a good solution. Great!

Could GenericMatrix perhaps also benefit from such explicit casting? Anything that improves the efficiency of row operations should have an immediate effect on overall efficiency of linear matroid code.

[tscrim]

Replying to Rudi:

Replying to tscrim:

Replying to Rudi: What I was hoping to do was take advantage of the known return types in the integer matrix to improve the Cython code. With the explicit casting, my computation goes from 13.6s to 3s.

That looks like a good solution. Great!

Could GenericMatrix perhaps also benefit from such explicit casting? Anything that improves the efficiency of row operations should have an immediate effect on overall efficiency of linear matroid code.

You cannot explicitly cast to something you don't know. Plus the arithmetic operations being done for IntegerMatrix are explicitly as C integers and not going through Python (much less Sage) at all. So I don't see how.

[tscrim]

I guess something you could do is also special case QQ and go directly with the C data using the mpq commands. It looks like the GF(2) and GF(3) matrices are already quite specialized. At least I would think that is the next most common field to work over (other than maybe GF(4)).

[jdemeyer]

Replying to tscrim:

So when I explicitly tell Cython that get should return an int

First of all, that's a bad idea: a C int has limited range (typically 32 bits) and you cannot guarantee that all entries fit in that.

then I have to do this change

     cdef bint is_nonzero(self, long r, long c) except -2:   # Not a Sage matrix operation
-        return self.get(r, c)
+        return self.get(r, c) != 0

otherwise I get the original error message in the ticket description.

I don't consider that a bug. The reason is that the conversion of a Python object to a bint essentially translates to int(bool(obj)) while the conversion to an int translates to int(obj). These are obviously different. But once on the C level, the types int and bint are exactly the same, so no conversion is done when casting int to bint.

[jdemeyer]

So the solution for really converting an int to a bint in Cython is indeed adding an explicit != 0.

[jdemeyer]

Should I close this ticket or do you want to recycle it?

[tscrim]

I will recycle this, at least for the speedups to IntegerMatrix, but maybe with some more tweaks or improvements or an implementation of a RationalMatrix.

However, the fact that bint and int are the same, then why do I get the error I am getting when I do not explicitly make it a bint with the != 0? From what you're saying, it should not result in an error. Unless there is some code elsewhere secretly expecting the result to be 0 or 1, which seems unlikely (or is bint allowed to take on different values)?

I am guessing instead of int, you would say we should have them all be Py_ssize_t?

[jdemeyer]

Replying to tscrim:

However, the fact that bint and int are the same, then why do I get the error I am getting when I do not explicitly make it a bint with the != 0? From what you're saying, it should not result in an error. Unless there is some code elsewhere secretly expecting the result to be 0 or 1, which seems unlikely (or is bint allowed to take on different values)?

Suppose that an entry in the matrix equals the Python integer -2. Currently (without applying any changes): when calling is_nonzero() on that -2, the return value is int(bool(-2)) = 1. With your changes, the return value for the same is_nonzero() call is just -2. The Cython wrapper interprets this -2 as an exception return value (that is what except -2 does). But there hasn't been an exception raised, which gives the SystemError.

[jdemeyer]

Replying to tscrim:

I am guessing instead of int, you would say we should have them all be Py_ssize_t?

I just looked at the Cython sources for lean_matrix.pyx and IntegerMatrix does indeed use int internally. Therefore, returning int in get() is correct in that sense.

However... this just shows that using IntegerMatrix in general looks dangerous: operations can overflow. This is actually documented in the docstring of IntegerMatrix:

        This class is mainly intended for use with the RegularMatroid class,
        so entries are assumed to be small integers. No
        overflow checking takes place!

[Stefan]

Thing is, the LeanMatrix? classes are internal datatypes. Regular matroids have matrices with entries equal to -1, 0, 1, and are such that pivoting preserves that property. For a potentially regular matroid on which you want to run is_valid(), we check this condition through pivoting, which means we get to a "bad subdeterminant" of size at most 2x2, hardly enough to cause overflows if the entries are 0,1, or -1.

We took some effort to keep the LeanMatrix? away from the end user (LinearMatroid?.Representation returns a Sage matrix, for instance, and the datatypes don't get imported into the Sage namespace). Speed is of the essence here, since we create and destroy tons of these (Sage matrices have a lot of overhead at creation time), and we do tons of row operations on them. If you dig deep enough to use IntegerMatrix? at this spot in the code, you'll have some idea of what you're getting yourself into.

Again, since you want to use IntegerMatrix? in the LinearMatroid? subclasses, and since you need a matrix where pivoting is safe, really the only use case is Regular Matroids where overflows are never an issue.

[tscrim]

Replying to jdemeyer:

Replying to tscrim:

However, the fact that bint and int are the same, then why do I get the error I am getting when I do not explicitly make it a bint with the != 0? From what you're saying, it should not result in an error. Unless there is some code elsewhere secretly expecting the result to be 0 or 1, which seems unlikely (or is bint allowed to take on different values)?

Suppose that an entry in the matrix equals the Python integer -2. Currently (without applying any changes): when calling is_nonzero() on that -2, the return value is int(bool(-2)) = 1. With your changes, the return value for the same is_nonzero() call is just -2. The Cython wrapper interprets this -2 as an exception return value (that is what except -2 does). But there hasn't been an exception raised, which gives the SystemError.

Ah, I see. Thank you for the explanation. I fully agree that it is not a bug.

[tscrim]

Replying to Stefan:

Thing is, the LeanMatrix? classes are internal datatypes. Regular matroids have matrices with entries equal to -1, 0, 1, and are such that pivoting preserves that property. For a potentially regular matroid on which you want to run is_valid(), we check this condition through pivoting, which means we get to a "bad subdeterminant" of size at most 2x2, hardly enough to cause overflows if the entries are 0,1, or -1.

I see. IntegerMatrix is a misnomer because it requires entries to be +-1 or 0, which is also undocumented. In fact, in the .. NOTE::, it says this works for LinearMatroid, but that is not quite the case. Thank you for the explanation.

Now I understand why it is working for the generic matrix over ZZ: it is secretly converting things to QQ when doing the matrix operations. My computation was not working using IntegerMatrix because of the assumption of +-1 for the nonzero entries.

I think this needs to be documented at the very least (which I will do here). I might also like to change the name of IntegerMatrix to perhaps UnitIntegerMatrix or PlusMinusOneMatrix. Thoughts?

Again, since you want to use IntegerMatrix? in the LinearMatroid? subclasses, and since you need a matrix where pivoting is safe, really the only use case is Regular Matroids where overflows are never an issue.

So it seems like I should write a version of LeanMatrix that directly uses mpq. I will start to work on that.

[jdemeyer]

Also: if the entries are really only 0, 1 or -1 then using an int to store the entries wastes memory. You could use int8_t (8 bits) instead of int (typically 32 bits). I'll leave it to you to decide whether it would be worth to change that.

[tscrim]

I will leave it to Rudi or Stefan to decide for that.

So I have gotten a somewhat working RationalMatrix version working, but I am having a bit of an issue with getting things to work like I expected. In Cython, are things pass-by-reference or pass-by-copy (like C)? In particular, I wanted to have an inlined get that returned the correct mpq_t entry, and I want to pass some of those around (and sometimes assign them to local variables). However, this does not seem to work. Does that mean I have to handle them always as pointers when I want to pass them around?

---

New commits:

​0d3105e	Explicit casts to int for IntegerMatrix.
​101663d	Documenting IntegerMatrix only should have +-1 or 0 entries.
​44538fa	Adding RationalMatrix that directly manipulates an mpq_t array.
​31a87a3	Removing get in place of index and some bug fixes.

[jdemeyer]

Replace

cdef int sval = int(s)

by

cdef int sval = s

The int() is superfluous and will only slow things down. Cython already implicitly converts when you assign to a C int.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​7effef0

Removing superfluous int's.

[tscrim]

Done. I didn't know that. Thanks.

[jdemeyer]

You will see that no function ever returns an mpq_t (or mpz_t or whatever). The typical calling convention is passing a return argument, for example the declaration of mpq_add:

void mpq_add(mpq_t sum, mpq_t addend1, mpq_t addend2)

[tscrim]

Yea, I noticed that and so I scrapped the get in place of the index for better maintainability should the internal ordering change. However, I do have some functions where I am passing an mpq_t, e.g., in pivot, and this is not working as I would expect. Also, I had tried to do something like

cdef mpq_t s = self._entries[self.index(i,j)]
mpq_add(s, s, t)

where t is some other mpq_t, but that wasn't working like I expected either.

[jdemeyer]

Can you move "Implement RationalMatrix" to a new ticket instead?

The commits that you added here to improve IntegerMatrix make sense independently.

[jdemeyer]

Replying to tscrim:

I think this needs to be documented at the very least (which I will do here). I might also like to change the name of IntegerMatrix to perhaps UnitIntegerMatrix or PlusMinusOneMatrix. Thoughts?

+1 to PlusMinusOneMatrix.

[tscrim]

Replying to jdemeyer:

Can you move "Implement RationalMatrix" to a new ticket instead?

The commits that you added here to improve IntegerMatrix make sense independently.

Yep, I can do that. I will do that when I get into my office tomorrow.

[tscrim]

#26269 for the RationalMatrix implementation.

[git]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

​fb64e43

Removing superfluous int's.

[tscrim]

I've split the ticket into the two parts. This part is now ready for review.

[tscrim]

(Essentially) Green patchbot.

[jdemeyer]

So how about renaming to PlusMinusOneMatrix? I'd like the opinion of Stefan on that.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​3c0195b

Renaming IntegerMatrix to PlusMinusOneMatrix.

[tscrim]

I did the refactoring. Stefan, Rudi, any objections?

[Stefan]

I am happy with the new name. Don’t have time for a long review.

[dkrenn]

I've looked through the code and I am happy with it, so LGTM.

However contributed to this ticket (mainly review), please insert your name in the reviewer field.

[tscrim]

Stefan, Rudi I added you as reviewers based on our conversions above. Jeroen is obviously a reviewer for looking at the Cython code (as well as his additional useful insights).