Opened 16 months ago
Closed 15 months ago
#29839 closed enhancement (fixed)
Implement `zero_pattern_matrix`
Reported by:  ghkliem  Owned by:  

Priority:  major  Milestone:  sage9.2 
Component:  basic arithmetic  Keywords:  incidence matrix, zero pattern 
Cc:  Merged in:  
Authors:  Jonathan Kliem, Travis Scrimshaw  Reviewers:  Travis Scrimshaw, Jonathan Kliem 
Report Upstream:  N/A  Work issues:  
Branch:  82c6b4a (Commits, GitHub, GitLab)  Commit:  82c6b4ad9afe2ae7be5e5c2fb2e47147c2b2bb89 
Dependencies:  Stopgaps: 
Description (last modified by )
We implement a method on matrices that gives a new matrix with specified base ring and entries 1 where and only where the initial matrix had zero entries, and 0 elsewhere.
Motivation: This method is useful to obtain the incidence matrix of a polyhedron from the slack matrix in #29837.
Along the way this ticket also fixes a bug:
M = Matrix(ZZ, [[0,1],[0,1]]) sage: M._mod_int(2).transpose()  TypeError Traceback (most recent call last) <ipythoninput48587f9dc80b2> in <module>() > 1 M._mod_int(Integer(2)).transpose() /srv/public/kliem/sage/local/lib/python3.7/sitepackages/sage/matrix/matrix_mod2_dense.pyx in sage.matrix.matrix_mod2_dense.Matrix_mod2_dense.transpose (build/cythonized/sage/matrix/matrix_mod2_dense.c:10873)() 1400 1 1401 """ > 1402 cdef Matrix_mod2_dense A = self.new_matrix(ncols = self._nrows, nrows = self._ncols) 1403 if self._nrows == 0 or self._ncols == 0: 1404 return A TypeError: Cannot convert sage.matrix.matrix_modn_dense_float.Matrix_modn_dense_float to sage.matrix.matrix_mod2_dense.Matrix_mod2_dense
To optimize this method, a matrix class can use the optimized get_is_zero_unsafe
.
In this ticket this is done for most matrix classes. We also use this new optimized zero check in places where it makes sense.
Change History (49)
comment:1 Changed 16 months ago by
comment:2 Changed 16 months ago by
 Branch set to public/29839
 Commit set to 43e6dfbd64c1372d0829ff18df7a9c7aa9491738
 Status changed from new to needs_review
New commits:
43e6dfb  implement _zero_matrix

comment:3 Changed 16 months ago by
zero_pattern_matrix?
comment:4 followup: ↓ 6 Changed 16 months ago by
Better anyway. With public or semiprivate? (leading underscore).
Btw, there is a bug here. I copied this from _mod_two
which apparently sets up the parent incorrect:
M = Matrix(ZZ, [[0,1],[0,1]]) sage: M._mod_int(2).transpose()  TypeError Traceback (most recent call last) <ipythoninput48587f9dc80b2> in <module>() > 1 M._mod_int(Integer(2)).transpose() /srv/public/kliem/sage/local/lib/python3.7/sitepackages/sage/matrix/matrix_mod2_dense.pyx in sage.matrix.matrix_mod2_dense.Matrix_mod2_dense.transpose (build/cythonized/sage/matrix/matrix_mod2_dense.c:10873)() 1400 1 1401 """ > 1402 cdef Matrix_mod2_dense A = self.new_matrix(ncols = self._nrows, nrows = self._ncols) 1403 if self._nrows == 0 or self._ncols == 0: 1404 return A TypeError: Cannot convert sage.matrix.matrix_modn_dense_float.Matrix_modn_dense_float to sage.matrix.matrix_mod2_dense.Matrix_mod2_dense sage: MS = matrix_space.MatrixSpace(IntegerModRing(2), self._nrows, self._ncols)
comment:5 Changed 16 months ago by
shouldn't def _zero_matrix(self):
(or rather named as in comment:3) be cpdef
, to facilitate calling from Cython?
comment:6 in reply to: ↑ 4 Changed 16 months ago by
Replying to ghkliem:
Better anyway. With public or semiprivate? (leading underscore).
Btw, there is a bug here. I copied this from
_mod_two
which apparently sets up the parent incorrect:M = Matrix(ZZ, [[0,1],[0,1]]) sage: M._mod_int(2).transpose()  TypeError Traceback (most recent call last) <ipythoninput48587f9dc80b2> in <module>() > 1 M._mod_int(Integer(2)).transpose() /srv/public/kliem/sage/local/lib/python3.7/sitepackages/sage/matrix/matrix_mod2_dense.pyx in sage.matrix.matrix_mod2_dense.Matrix_mod2_dense.transpose (build/cythonized/sage/matrix/matrix_mod2_dense.c:10873)() 1400 1 1401 """ > 1402 cdef Matrix_mod2_dense A = self.new_matrix(ncols = self._nrows, nrows = self._ncols) 1403 if self._nrows == 0 or self._ncols == 0: 1404 return A TypeError: Cannot convert sage.matrix.matrix_modn_dense_float.Matrix_modn_dense_float to sage.matrix.matrix_mod2_dense.Matrix_mod2_dense sage: MS = matrix_space.MatrixSpace(IntegerModRing(2), self._nrows, self._ncols)
why can't you use mod
?
sage: M = Matrix(ZZ, [[0,1],[0,1]]) sage: M.mod(2) [0 1] [0 1]
comment:7 Changed 16 months ago by
Because it is slow. I just copied the construction of the parent from _mod_two
.
I really don't know why, but change_ring
is really slow (e.g. for ZZ
> GF(2)
). It should be improved, but maybe not in this ticket.
_mod_two
uses the ring with two elements to obtain the parent. However, get_matrix_class
gives you modn
in this case. It gives you mod2
if you pass the field with two elements.
comment:8 Changed 16 months ago by
And M.mod
is really not what I want. I just illustrated this bug.
comment:9 Changed 16 months ago by
 Description modified (diff)
 Keywords pattern added; entries removed
comment:10 Changed 16 months ago by
 Commit changed from 43e6dfbd64c1372d0829ff18df7a9c7aa9491738 to 7b708307f2aa4860299a9b7470e1344b88b59c12
comment:11 Changed 16 months ago by
 Summary changed from Implement `_zero_matrix` for dense integer and rational matrices to Implement `zero_pattern_matrix` for dense integer and rational matrices
comment:12 Changed 16 months ago by
 Status changed from needs_review to needs_work
I really want a new integer matrix.
comment:13 Changed 16 months ago by
 Commit changed from 7b708307f2aa4860299a9b7470e1344b88b59c12 to 62f1ff9bcbf0bc221d209ef92e893c375b76bb02
Branch pushed to git repo; I updated commit sha1. New commits:
62f1ff9  use integer matrix for zero_pattern_matrix

comment:14 Changed 16 months ago by
 Description modified (diff)
 Status changed from needs_work to needs_review
comment:15 followup: ↓ 16 Changed 16 months ago by
Have you tried doing either x = 1  x
or x = not x
? These should be fast as they avoid intermediate function calls and the first should work for any ring (maybe replace 1
with one
where one = self.base_ring().one()
).
comment:16 in reply to: ↑ 15 Changed 16 months ago by
Replying to tscrim:
Have you tried doing either
x = 1  x
orx = not x
? These should be fast as they avoid intermediate function calls and the first should work for any ring (maybe replace1
withone
whereone = self.base_ring().one()
).
When? As in slack_matrix.apply_map
?
This is slow, whatever I do.
If slack_matrix
is the slack_matrix
of the polytopes.associahedron(['A', 9])
obtained e.g. with
sage: def get_slack(P): ....: self = P ....: Vrep_matrix = matrix(self.base_ring(), self.Vrepresentation()) ....: Hrep_matrix = matrix(self.base_ring(), self.Hrepresentation()) ....: hom_helper = matrix(self.base_ring(), [1 if v.is_vertex() else 0 for v ....: in self.Vrepresentation()]) ....: hom_Vrep = hom_helper.stack(Vrep_matrix.transpose()) ....: slack_matrix = (Hrep_matrix * hom_Vrep).transpose() ....: return slack_matrix sage: slack_matrix = get_slack(polytopes.associahedron(['A', 9])
than the solution of this ticket needs about 15 ms. In comparison to
sage: %time incidence_matrix = M.apply_map(lambda y: not y, R=M.base_ring()) CPU times: user 281 ms, sys: 15.8 ms, total: 296 ms Wall time: 296 ms
(which is 100 ms slower than the trivial lambda y: y
).
comment:17 Changed 16 months ago by
For the x = 1  x
, I was thinking the matrix was already a 01 matrix (and perhaps you were trying to do it inplace). This is also slower than simply doing this?
cdef Py_ssize_t i, j for i from 0 <= i < self._nrows: for j from 0 <= j < self._ncols: if not self.get_unsafe(i, j): M.set_unsafe_si(i, j, 1)
This is more generic, but perhaps it is slower for special cases (e.g., rationals) because it goes partially through Python.
I have some questions:
Should this be lifted up to a generic matrix method with a generic implementation of get_is_zero_unsafe
?
Do you want to be an integer matrix or a mod2 matrix? The latter seems like it would be more memory efficient and possibly faster.
I am not sure it is safe to assume the initialized matrix is 0. Is this documented somewhere?
comment:18 Changed 16 months ago by
Good point. It is documented that the entries are not initialized.
So the following would be better:
cdef Matrix_integer_dense M = self._new(self._nrows, self._ncols) cdef Py_ssize_t i, j for i from 0 <= i < self._nrows: for j from 0 <= j < self._ncols: if self.get_is_zero_unsafe(i, j): M.set_unsafe_si(i, j, 1) else: M.set_unsafe_si(i, j, 0)
For integers the difference between if not self.get_unsafe(i, j)
and if get_is_zero_unsafe(i, j)
is a factor of three. For rationals the difference is a factor of about 20.
comment:19 Changed 16 months ago by
M.set_unsafe_si(i, j, self.get_is_zero_unsafe(i, j))
is a bit faster yet and almost the speed of not correctly initializing (26 ms vs 23 ms).
comment:20 Changed 16 months ago by
 Commit changed from 62f1ff9bcbf0bc221d209ef92e893c375b76bb02 to b5837452ce2e3584c1ec2ca6f98ab4a2dd17de1c
Branch pushed to git repo; I updated commit sha1. New commits:
b583745  make methods more generic

comment:21 Changed 16 months ago by
 Description modified (diff)
 Summary changed from Implement `zero_pattern_matrix` for dense integer and rational matrices to Implement `zero_pattern_matrix`
the
comment:22 Changed 16 months ago by
 Commit changed from b5837452ce2e3584c1ec2ca6f98ab4a2dd17de1c to 535dc85e9d8a04baf28a646e11f80651b9c16e8a
Branch pushed to git repo; I updated commit sha1. New commits:
535dc85  note for developers

comment:23 Changed 16 months ago by
If you feel like get_is_zero_unsafe
should be improved in other classes as well, please tell me. At the moment I would just leave it being optimized for those two classes that I need.
comment:24 followup: ↓ 26 Changed 16 months ago by
Thank you. It is looking a lot better. One other place where it might make sense to do additional specialized get_is_zero_unsafe
would be cyclo, gf2e, gfpn, and modn matrices, but these can wait if they turn out to be harder to do (I can also do some of them tomorrow). A few more little things:
I don't understand this line:
cdef object one = M._coerce_element(1)
Why not use R.one()
? Also, similar to one of my previous questions, do we want this to be over the general ring? In principle, the matrix could be over a rng, and so the result might not be valid over R
but the operation/output still makes sense. I think whatever output type we want (either as a ZZ
matrix or a Zmod(2)
matrix) is fine as long as we document it.
 NOTE:: + .. NOTE::  This method can be optimized by improving :meth:`get_is_zero_unsafe`` for derived matrix classes. + This method can be optimized by improving :meth:`get_is_zero_unsafe` for derived matrix classes.
.. warning:: +.. WARNING::
Do we want get_is_zero_unsafe
to return a bint
instead of an int
? I don't hold a strong opinion, but I generally prefer the most restrictive/informative return type possible.
comment:25 Changed 16 months ago by
I prefer int
for the following reason: It can be useful, if we guarantee that the value is either 1
or 0
. I'm not using it like this anymore, but I wouldn't burn bridges. This way we really guarantee that the method returns 1
or 0
. If we return a bint
this would leave the impression that the exact value doesn't matter.
sage: cython(''' ....: def foo(int a): ....: cdef bint b = a ....: return 2  b ....: ''')
This is an example that bint
really can take any integer value. If you run foo(10)
, it will return 8
.
comment:26 in reply to: ↑ 24 Changed 16 months ago by
Replying to tscrim:
Thank you. It is looking a lot better. One other place where it might make sense to do additional specialized
get_is_zero_unsafe
would be cyclo, gf2e, gfpn, and modn matrices, but these can wait if they turn out to be harder to do (I can also do some of them tomorrow). A few more little things:I don't understand this line:
cdef object one = M._coerce_element(1)Why not use
R.one()
? Also, similar to one of my previous questions, do we want this to be over the general ring? In principle, the matrix could be over a rng, and so the result might not be valid overR
but the operation/output still makes sense.
I don't find an example of a matrix with base ring not a ring. It is mentioned in the constructor that such things exists, but I can't find it. Can you help me out?
We could have another optional argument, where you specify what you want to fill the matrix with, if not ones.
I think whatever output type we want (either as a
ZZ
matrix or aZmod(2)
matrix) is fine as long as we document it. NOTE:: + .. NOTE::  This method can be optimized by improving :meth:`get_is_zero_unsafe`` for derived matrix classes. + This method can be optimized by improving :meth:`get_is_zero_unsafe` for derived matrix classes... warning:: +.. WARNING::Do we want
get_is_zero_unsafe
to return abint
instead of anint
? I don't hold a strong opinion, but I generally prefer the most restrictive/informative return type possible.
comment:27 Changed 16 months ago by
 Commit changed from 535dc85e9d8a04baf28a646e11f80651b9c16e8a to 00aed0243abfd6a37b3ec9cc7c7e8c775abfab9d
comment:28 Changed 16 months ago by
 Status changed from needs_review to needs_work
 Work issues set to documentation, more classes, more general baserings
comment:29 Changed 16 months ago by
 Commit changed from 00aed0243abfd6a37b3ec9cc7c7e8c775abfab9d to fb49574339d302455e9e919d04fa6546f61c92ec
Branch pushed to git repo; I updated commit sha1. New commits:
fb49574  typo

comment:30 Changed 16 months ago by
Sage explicitly only allows rings to be passed to matrices, despite what the documentation currently says:
sage: S = Semigroups().example() sage: A = S.algebra(QQ) sage: M = matrix(A, [[A.zero()]])  TypeError Traceback (most recent call last) <ipythoninput45bed397dad10> in <module>() > 1 M = matrix(A, [[A.zero()]]) /home/uqtscrim/sage/local/lib/python3.7/sitepackages/sage/matrix/constructor.pyx in sage.matrix.constructor.matrix (build/cythonized/sage/matrix/constructor.c:2479)() 633 """ 634 immutable = kwds.pop('immutable', False) > 635 M = MatrixArgs(*args, **kwds).matrix() 636 if immutable: 637 M.set_immutable() /home/uqtscrim/sage/local/lib/python3.7/sitepackages/sage/matrix/args.pyx in sage.matrix.args.MatrixArgs.matrix (build/cythonized/sage/matrix/args.c:7668)() 650 True 651 """ > 652 self.finalize() 653 654 cdef Matrix M /home/uqtscrim/sage/local/lib/python3.7/sitepackages/sage/matrix/args.pyx in sage.matrix.args.MatrixArgs.finalize (build/cythonized/sage/matrix/args.c:10264)() 932 933 if self.space is None: > 934 self.space = MatrixSpace(self.base, self.nrows, self.ncols, 935 sparse=self.sparse, **self.kwds) 936 /home/uqtscrim/sage/local/lib/python3.7/sitepackages/sage/misc/classcall_metaclass.pyx in sage.misc.classcall_metaclass.ClasscallMetaclass.__call__ (build/cythonized/sage/misc/classcall_metaclass.c:1741)() 320 """ 321 if cls.classcall is not None: > 322 return cls.classcall(cls, *args, **kwds) 323 else: 324 # Fast version of type.__call__(cls, *args, **kwds) /home/uqtscrim/sage/local/lib/python3.7/sitepackages/sage/matrix/matrix_space.py in __classcall__(cls, base_ring, nrows, ncols, sparse, implementation, **kwds) 489 """ 490 if base_ring not in _Rings: > 491 raise TypeError("base_ring (=%s) must be a ring"%base_ring) 492 nrows = int(nrows) 493 if ncols is None: TypeError: base_ring (=Algebra of An example of a semigroup: the left zero semigroup over Rational Field) must be a ring
So that might not something we need to worry about right now. Irregardless, I am still not convinced that we should be returning a matrix over the base ring and not either Z or Z_{2}.
I see your argument for not using bint
; it is actually a little disheartening to see that code block; but I see that as garbagein, garbageout. My counter point would be the is_zero_unsage
is suppose to be recording a true/false value, which is what the bint
is suppose to be for (equivalent to a bool
in C). Plus it is standard that 1
is True
. So even if it is not enforacble, it still carries contextual significance.
comment:31 Changed 16 months ago by
As for other rings I think I have questions with all of them:
gfpn
isget_unsafe_int
guaranteed to be0
if the entry is zero, or might it be0
mod cardinality of field, if not, where do I find the cardinality? similar for
gf2e
andmodn
cyclo
(I just realized this is probably easy)
EDIT:
For gf2e
this is indirectly confirmed in rings/finite_rings/element_givaro.pyx
in fetch_int
.
comment:32 Changed 16 months ago by
 Commit changed from fb49574339d302455e9e919d04fa6546f61c92ec to 8482c5ae0107d1715c56c4bc29fbbc907b1ef29c
comment:33 Changed 16 months ago by
 Work issues changed from documentation, more classes, more general baserings to more classes
comment:34 Changed 16 months ago by
 Commit changed from 8482c5ae0107d1715c56c4bc29fbbc907b1ef29c to 169fa6c85de38be107d75481372f71cbf8b11c98
Branch pushed to git repo; I updated commit sha1. New commits:
13ff08b  Adding faster get_is_zero_unsafe for modn matrices.

20bd4e0  Adding get_is_zero_unsafe() and __bool__ to gf2e matrices.

cbe266a  Adding get_is_zero_unsafe to cyclotomic matrices.

e506816  Using get_is_unsafe_zero() in more places.

169fa6c  Adding get_is_zero_unsafe() to gfpn matrices.

comment:35 Changed 16 months ago by
 Reviewers set to Travis Scrimshaw
 Status changed from needs_work to needs_review
 Work issues more classes deleted
Okay, this is an implementation for all of the specialized dense matrices I mentioned before.
I really don't like how in the gf2e it redirects through python for such a key feature with get_unsafe
. This needs to be removed and is now #29853. Likely the way to do it will be based on what is done for gfpn.
Some additions of better __bool__
methods:
sage: K.<z> = CyclotomicField(3) sage: A = matrix(K, 4, 3, [0, z, 2, 2*z + 2, 2*z, z, z, 1z, 2+3*z, z, 1+z, 0]) sage: %timeit bool(A) 10000000 loops, best of 5: 86.2 ns per loop sage: %timeit bool(Z) 10000000 loops, best of 5: 152 ns per loop
versus before:
sage: %timeit bool(A) The slowest run took 70.37 times longer than the fastest. This could mean that an intermediate result is being cached. 100000 loops, best of 5: 8.49 µs per loop sage: Z = matrix(K, 4, 3, 0) sage: %timeit bool(Z) The slowest run took 7.43 times longer than the fastest. This could mean that an intermediate result is being cached. 10000 loops, best of 5: 50 µs per loop
I also made changes in other places where get_is_unsafe_zero()
could be used. Like the above and in your use case, this likely has (potentially significantly) sped things up.
comment:36 Changed 16 months ago by
 Commit changed from 169fa6c85de38be107d75481372f71cbf8b11c98 to 946237697b17bf1a2024326f854b1ef6e27dec04
comment:37 Changed 16 months ago by
 Commit changed from 946237697b17bf1a2024326f854b1ef6e27dec04 to 94febdeb4ae5782a97793d9446fe8951f0c78bf9
Branch pushed to git repo; I updated commit sha1. New commits:
94febde  Custom get_is_zero_unsafe() for modn sparse matrices.

comment:38 Changed 16 months ago by
 Commit changed from 94febdeb4ae5782a97793d9446fe8951f0c78bf9 to bf0352fdad5db674c3c2aa36f188f996600f7ee3
Branch pushed to git repo; I updated commit sha1. New commits:
bf0352f  See if we can check trivial zeros first for symbolic matrices.

comment:39 Changed 16 months ago by
 Description modified (diff)
Okay, that takes care of the sparse matrices and symbolic matrices.
Slightly surprisingly, the gap matrices don't benefit from not going through Sage (well, at least for integer matrices).
I am done now.
comment:40 Changed 16 months ago by
 Commit changed from bf0352fdad5db674c3c2aa36f188f996600f7ee3 to 2a08abe750c622e191a5849979bde1a041a31d52
Branch pushed to git repo; I updated commit sha1. New commits:
2a08abe  small mistakes

comment:41 Changed 16 months ago by
 Reviewers changed from Travis Scrimshaw to Travis Scrimshaw, Jonathan Kliem
IMO this is good now. The doctests should also pass.
comment:42 Changed 16 months ago by
 Commit changed from 2a08abe750c622e191a5849979bde1a041a31d52 to 333c4baf305e19b3a12e99a48e3140459521f339
Branch pushed to git repo; I updated commit sha1. New commits:
333c4ba  wrong doctest

comment:43 followup: ↓ 45 Changed 16 months ago by
Thank you. Two more little things:
INPUT:   `ring`  (optional); base ring of the output; default is `ZZ` +  `ring`  (optional); base ring of the output; default is ``ZZ``
and we should add a test for the empty row/column and empty matrix for Matrix_gf2e_dense.__bool__
(as that was a good catch). I can do this tomorrow.
comment:44 Changed 16 months ago by
 Commit changed from 333c4baf305e19b3a12e99a48e3140459521f339 to 82c6b4ad9afe2ae7be5e5c2fb2e47147c2b2bb89
Branch pushed to git repo; I updated commit sha1. New commits:
82c6b4a  docstring

comment:45 in reply to: ↑ 43 Changed 16 months ago by
Replying to tscrim:
and we should add a test for the empty row/column and empty matrix for
Matrix_gf2e_dense.__bool__
(as that was a good catch).
Well the bot was segfaulting (and me to, when running the test).
comment:46 Changed 16 months ago by
 Status changed from needs_review to positive_review
Okay, then positive review. Thank you.
comment:47 Changed 16 months ago by
I like that you made a polyhedron only ticket, into a more general improvement for other parts as well.
comment:48 Changed 15 months ago by
 Description modified (diff)
(Slightly amending ticket description, please check that works for you.).
comment:49 Changed 15 months ago by
 Branch changed from public/29839 to 82c6b4ad9afe2ae7be5e5c2fb2e47147c2b2bb89
 Resolution set to fixed
 Status changed from positive_review to closed
Better suggestions for a name would be welcome.
_zero_matrix
is somewhat strange.