Opened 4 years ago
Closed 3 years ago
#26237 closed defect (fixed)
Speed up IntegerMatrix
Reported by:  tscrim  Owned by:  

Priority:  major  Milestone:  sage8.8 
Component:  matroid theory  Keywords:  
Cc:  Stefan, Rudi, yomcat, msaaltink  Merged in:  
Authors:  Travis Scrimshaw  Reviewers:  Daniel Krenn, Jeroen Demeyer, Stefan Van Zwam, Rudi Pendavingh 
Report Upstream:  N/A  Work issues:  
Branch:  3c0195b (Commits, GitHub, GitLab)  Commit:  3c0195bc7302b36431d1c4715582fd33ed5cc0c0 
Dependencies:  Stopgaps: 
Description (last modified by )
We explicitly declare things to be int
type so Cython can generate better C code.
Change History (49)
comment:1 Changed 4 years ago by
comment:2 Changed 4 years ago by
Hmm...strange. Let me try from a fresh session to see if I somehow unintentionally corrupted things.
comment:3 Changed 4 years ago by
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.
comment:4 followup: ↓ 6 Changed 4 years ago by
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.
comment:5 followup: ↓ 9 Changed 4 years ago by
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})]
comment:6 in reply to: ↑ 4 Changed 4 years ago by
Replying to tscrim:
So what was going wrong is that there is a
bint is_nonzero
method which simply wasreturn 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"?
comment:7 Changed 4 years ago by
Yes, that is correct (well, an int
to a bint
to be specific).
comment:8 followup: ↓ 10 Changed 4 years ago by
At least, once I changed that to an explicit check != 0, it fixed the problem.
comment:9 in reply to: ↑ 5 ; followup: ↓ 13 Changed 4 years ago by
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?
comment:10 in reply to: ↑ 8 ; followup: ↓ 12 Changed 4 years ago by
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...
comment:11 followup: ↓ 14 Changed 4 years ago by
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().
comment:12 in reply to: ↑ 10 ; followup: ↓ 18 Changed 4 years ago by
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.
comment:13 in reply to: ↑ 9 Changed 4 years ago by
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.
comment:14 in reply to: ↑ 11 ; followup: ↓ 15 Changed 4 years ago by
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
diff git a/src/sage/matroids/lean_matrix.pxd b/src/sage/matroids/lean_matrix.pxd index c284f6a..122dab1 100644
a b cdef class QuaternaryMatrix(LeanMatrix): 102 102 cdef class IntegerMatrix(LeanMatrix): 103 103 cdef int* _entries 104 104 105 cdef inline get(self, long r, long c) # Not a Sage matrix operation105 cdef inline int get(self, long r, long c) # Not a Sage matrix operation 106 106 cdef inline void set(self, long r, long c, int x) # Not a Sage matrix operation 107 107 108 108 cdef inline long row_len(self, long i) except 1 # Not a Sage matrix operation 
src/sage/matroids/lean_matrix.pyx
diff git a/src/sage/matroids/lean_matrix.pyx b/src/sage/matroids/lean_matrix.pyx index 93a32bc..d6022c2 100644
a b cdef class IntegerMatrix(LeanMatrix): 2844 2844 """ 2845 2845 return "IntegerMatrix instance with " + str(self._nrows) + " rows and " + str(self._ncols) + " columns" 2846 2846 2847 cdef inline get(self, long r, long c): # Not a Sage matrix operation2847 cdef inline int get(self, long r, long c): # Not a Sage matrix operation 2848 2848 return self._entries[r * self._ncols + c] 2849 2849 2850 2850 cdef inline void set(self, long r, long c, int x): # Not a Sage matrix operation … … cdef class IntegerMatrix(LeanMatrix): 2877 2877 return 0 2878 2878 2879 2879 cdef bint is_nonzero(self, long r, long c) except 2: # Not a Sage matrix operation 2880 return self.get(r, c) 2880 return self.get(r, c) != 0 2881 2881 2882 2882 cdef LeanMatrix copy(self): # Deprecated Sage matrix operation 2883 2883 cdef IntegerMatrix M = IntegerMatrix(self._nrows, self._ncols) … … cdef class IntegerMatrix(LeanMatrix): 2982 2982 ignored. 2983 2983 """ 2984 2984 cdef long i 2985 cdef int sval 2985 2986 if s is None: 2986 2987 for i from 0 <= i < self._ncols: 2987 2988 self.set(x, i, self.get(x, i) + self.get(y, i)) 2988 2989 else: 2990 sval = int(s) 2989 2991 for i from 0 <= i < self._ncols: 2990 self.set(x, i, self.get(x, i) + s * self.get(y, i))2992 self.set(x, i, self.get(x, i) + sval * self.get(y, i)) 2991 2993 return 0 2992 2994 2993 2995 cdef int swap_rows_c(self, long x, long y) except 1: … … cdef class IntegerMatrix(LeanMatrix): 3010 3012 compatibility, and is ignored. 3011 3013 """ 3012 3014 cdef long i 3015 cdef int sval = int(s) 3013 3016 for i from 0 <= i < self._ncols: 3014 self.set(x, i, s * self.get(x, i))3017 self.set(x, i, sval * self.get(x, i)) 3015 3018 return 0 3016 3019 3017 3020 cdef int rescale_column_c(self, long y, s, bint start_row) except 1: … … cdef class IntegerMatrix(LeanMatrix): 3020 3023 compatibility, and is ignored. 3021 3024 """ 3022 3025 cdef long j 3026 cdef int sval = int(s) 3023 3027 for j from 0 <= j < self._nrows: 3024 self.set(j, y, self.get(j, y) * s )3028 self.set(j, y, self.get(j, y) * sval) 3025 3029 return 0 3026 3030 3027 3031 cdef int pivot(self, long x, long y) except 1: # Not a Sage matrix operation
comment:15 in reply to: ↑ 14 ; followup: ↓ 16 Changed 4 years ago by
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.
comment:16 in reply to: ↑ 15 Changed 4 years ago by
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.
comment:17 Changed 4 years ago by
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)
).
comment:18 in reply to: ↑ 12 Changed 4 years ago by
Replying to tscrim:
So when I explicitly tell Cython that
get
should return anint
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) != 0otherwise 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
.
comment:19 Changed 4 years ago by
So the solution for really converting an int
to a bint
in Cython is indeed adding an explicit != 0
.
comment:20 Changed 4 years ago by
Should I close this ticket or do you want to recycle it?
comment:21 followups: ↓ 23 ↓ 24 Changed 4 years ago by
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
?
comment:22 Changed 4 years ago by
 Description modified (diff)
 Summary changed from SystemError with LinearMatroid with an IntegerMatrix to Speed up IntegerMatrix
comment:23 in reply to: ↑ 21 ; followup: ↓ 26 Changed 4 years ago by
Replying to tscrim:
However, the fact that
bint
andint
are the same, then why do I get the error I am getting when I do not explicitly make it abint
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 isbint
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
.
comment:24 in reply to: ↑ 21 Changed 4 years ago by
Replying to tscrim:
I am guessing instead of
int
, you would say we should have them all bePy_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!
comment:25 followup: ↓ 27 Changed 4 years ago by
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.
comment:26 in reply to: ↑ 23 Changed 4 years ago by
Replying to jdemeyer:
Replying to tscrim:
However, the fact that
bint
andint
are the same, then why do I get the error I am getting when I do not explicitly make it abint
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 isbint
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 isint(bool(2)) = 1
. With your changes, the return value for the sameis_nonzero()
call is just2
. The Cython wrapper interprets this2
as an exception return value (that is whatexcept 2
does). But there hasn't been an exception raised, which gives theSystemError
.
Ah, I see. Thank you for the explanation. I fully agree that it is not a bug.
comment:27 in reply to: ↑ 25 ; followup: ↓ 37 Changed 4 years ago by
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.
comment:28 Changed 4 years ago by
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.
comment:29 Changed 4 years ago by
comment:30 Changed 4 years ago by
 Branch set to public/matroids/speedup_integer_matrix26237
 Commit set to 31a87a3194b9c77febb961c4821d1873a4b01b2a
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 passbyreference or passbycopy (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.

comment:31 Changed 4 years ago by
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
.
comment:32 Changed 4 years ago by
 Commit changed from 31a87a3194b9c77febb961c4821d1873a4b01b2a to 7effef06c3c1f821c1bfa3316d55e88c0026333a
Branch pushed to git repo; I updated commit sha1. New commits:
7effef0  Removing superfluous int's.

comment:33 Changed 4 years ago by
Done. I didn't know that. Thanks.
comment:34 Changed 4 years ago by
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)
comment:35 Changed 4 years ago by
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.
comment:36 followup: ↓ 38 Changed 4 years ago by
Can you move "Implement RationalMatrix
" to a new ticket instead?
The commits that you added here to improve IntegerMatrix
make sense independently.
comment:37 in reply to: ↑ 27 Changed 4 years ago by
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 perhapsUnitIntegerMatrix
orPlusMinusOneMatrix
. Thoughts?
+1 to PlusMinusOneMatrix
.
comment:38 in reply to: ↑ 36 Changed 4 years ago by
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.
comment:39 Changed 4 years ago by
 Description modified (diff)
#26269 for the RationalMatrix
implementation.
comment:40 Changed 4 years ago by
 Commit changed from 7effef06c3c1f821c1bfa3316d55e88c0026333a to fb64e43cd4e5eaa28b3156b39c7cb42909ac0075
Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:
fb64e43  Removing superfluous int's.

comment:41 Changed 4 years ago by
 Status changed from new to needs_review
I've split the ticket into the two parts. This part is now ready for review.
comment:42 Changed 4 years ago by
(Essentially) Green patchbot.
comment:43 Changed 4 years ago by
So how about renaming to PlusMinusOneMatrix
? I'd like the opinion of Stefan on that.
comment:44 Changed 4 years ago by
 Commit changed from fb64e43cd4e5eaa28b3156b39c7cb42909ac0075 to 3c0195bc7302b36431d1c4715582fd33ed5cc0c0
Branch pushed to git repo; I updated commit sha1. New commits:
3c0195b  Renaming IntegerMatrix to PlusMinusOneMatrix.

comment:45 Changed 4 years ago by
I did the refactoring. Stefan, Rudi, any objections?
comment:46 Changed 4 years ago by
I am happy with the new name. Don’t have time for a long review.
comment:47 Changed 3 years ago by
 Reviewers set to Daniel Krenn
 Status changed from needs_review to positive_review
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.
comment:48 Changed 3 years ago by
 Milestone changed from sage8.4 to sage8.8
 Reviewers changed from Daniel Krenn to Daniel Krenn, Jeroen Demeyer, Stefan Van Zwam, Rudi Pendavingh
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).
comment:49 Changed 3 years ago by
 Branch changed from public/matroids/speedup_integer_matrix26237 to 3c0195bc7302b36431d1c4715582fd33ed5cc0c0
 Resolution set to fixed
 Status changed from positive_review to closed
Works for me: