Opened 9 years ago
Closed 9 years ago
#15245 closed enhancement (fixed)
Pfaffian of a skew-symmetric matrix
Reported by: | Darij Grinberg | Owned by: | |
---|---|---|---|
Priority: | major | Milestone: | sage-5.13 |
Component: | combinatorics | Keywords: | matrix, sage-combinat, pfaffian |
Cc: | Sage Combinat CC user, spancratz, Thierry Monteil | Merged in: | sage-5.13.beta1 |
Authors: | Darij Grinberg | Reviewers: | Travis Scrimshaw |
Report Upstream: | N/A | Work issues: | |
Branch: | Commit: | ||
Dependencies: | #14117 | Stopgaps: |
Description (last modified by )
I couldn't believe my eyes when I saw we don't have the Pfaffian implemented in Sage.
Attached is an implementation that computes it over any ring in the combinatorial way using perfect matchings. This is probably not an optimal algorithm (and I feel like a lot could already be gained by improving up the PerfectMatchings?(n) iterator without even changing the algorithm) but it's enough for my combinatorial needs.
Other than this, the patch makes perfect matchings iterable (no, they weren't) and improves some docstrings related to rook polynomials. The #14117 dependency is because both patches edit matrix2.py and would probably cause some fuzz.
Apply:
Attachments (1)
Change History (26)
comment:1 Changed 9 years ago by
Status: | new → needs_review |
---|
comment:2 Changed 9 years ago by
Dependencies: | → #14117 |
---|
comment:4 follow-up: 7 Changed 9 years ago by
Hi Travis,
thanks once again for the reviews! I fixed all of your issues apart from not using is_skew_symmetric()
because that method doesn't check diagonal entries to be 0 (it only checks them to satisfy x = -x, which is not the same in characteristic 2):
sage: M = Matrix(Zmod(2), [[0,1],[1,1]]) sage: M.is_skew_symmetric() True
If you ask me, this is a bug, but I don't want to grasp into yet another hornet's nest. But thanks for having me look at my check code again; it contained an ugly indentation error...
I was pretty sure that matrix algorithms are the least busy part of the code, given how long #14117 took to be reviewed. In hindsight I shouldn't have done unrelated docstring fixes, but the rook stuff was just pretty close and caught my eyes.
Best regards,
Darij
comment:5 Changed 9 years ago by
Description: | modified (diff) |
---|---|
Status: | needs_review → positive_review |
comment:7 Changed 9 years ago by
Replying to darij:
thanks once again for the reviews! I fixed all of your issues apart from not using
is_skew_symmetric()
because that method doesn't check diagonal entries to be 0 (it only checks them to satisfy x = -x, which is not the same in characteristic 2):
I think that's the definition of skew symmetric/antisymmetric: that A^{T}=-A. It just happens to be the case that the concepts "symmetric" and "antisymmetric" coincide in characteristic 2.
In other words: skew symmetric matrices don't have to have 0 on their diagonal. Note that the terminology comes from bilinear forms, where alternating means (v,v)=0, antisymmetric means (v,w)=-(w,v) and symmetric means (v,w)=(w,v). "alternating" is not "antisymmetric" in characteristic 2.
comment:8 Changed 9 years ago by
That's the hornet's nest I was talking about. With the definition you give, the Pfaffian lacks many of its nice properties like squaring to the determinant. It is the definition used on http://en.wikipedia.org/wiki/Skew-symmetric_matrix but not the definition used on http://en.wikipedia.org/wiki/Pfaffian . Since mathematicians can't agree, I figured it is easier to check for the kind of skew-symmetry needed in the definition of the Pfaffian rather than push that definition into the rest of the code. Is there a better solution?
comment:9 Changed 9 years ago by
Status: | positive_review → needs_work |
---|
comment:10 Changed 9 years ago by
Status: | needs_work → needs_review |
---|
darij: it looks like you reviewed your own patch. I'd say this still needs a "formal" review.
comment:11 Changed 9 years ago by
Oops -- I just realized what Travis gave were comments, not a review. Sorry!
comment:12 Changed 9 years ago by
Hey Darij,
Sorry for letting this slip away.
For the skew-symmetric, how about instead calling is_skew_symmetric()
and then checking that the diagonal entries are 0 with
all(d == 0 for d in self.diagonal())
One more minor thing: the AUTHORS:
block shouldn't be indented.
Best,
Travis
comment:13 follow-up: 14 Changed 9 years ago by
I've given the is_skew_symmetric
method an additional keyword variable now. All the rest is OK?
comment:14 Changed 9 years ago by
Replying to darij:
I've given the
is_skew_symmetric
method an additional keyword variable now. All the rest is OK?
That keyword should really be is_alternating
, so the more appropriate thing would be to supply a method is_alternating
.
If pfaffians need alternating rather than skew-symmetric (I haven't checked) then the confusion in terminology just comes from the fact that people looking at pfaffians haven't considered characteristic 2. That kind of thing happens all the time, and it's the kind of thing that computer algebra systems need to be a little more pedantic about than math literature, since you don't get to say "in this paper, with <THIS> we mean <SOMETHING ELSE>".
It may well be that pfaffians are simply not all that useful in characteristic 2, so that people didn't bother with them (Cayley certainly wouldn't have).
Changed 9 years ago by
Attachment: | trac_15245-pfaffian-dg.patch added |
---|
new version, separating skew-symmetry from alternatingness systematically
comment:15 Changed 9 years ago by
Done. A number of people seem to be lax about characteristic 2, among them Knuth from whom I had expected this the least.
comment:16 follow-up: 17 Changed 9 years ago by
I'm happy with it (having an is_alternating()
method). Nils?
comment:17 Changed 9 years ago by
Replying to tscrim:
I'm happy with it (having an
is_alternating()
method). Nils?
Yep. It's a pedantic difference, but since this gets exposed in our official matrix API I think it's worth being precise.
comment:18 Changed 9 years ago by
Reviewers: | → Travis Scrimshaw |
---|---|
Status: | needs_review → positive_review |
Then it's a positive review.
comment:19 Changed 9 years ago by
Oh, a comment that may be worthwhile for future work:
The generic implementations of is_skewsymmetric and is_alternating are probably horribly slow compared to what one can do on specific classes (see, e.g. #15104). Since the difference between is_alternating and is_skewsymmetric only is apparent in characteristic 2, we'd probably get better performance if one of the two calls the other if required. If the specific implementation only applies to characteristic not equal to 2, it only needs to implement one fast method and the generic other one will call it.
Probably is_skewsymmetric
is the better choice for being the "main" method, because it's the more widely used term. We'd get something like:
def is_alternating(self): if self.base_ring().characteristic() !=2: return self.is_skewsymmetric() <rest of code>
comment:20 follow-up: 21 Changed 9 years ago by
Not sure about it. is_skewsymmetric
should be a tad slower than is_alternating
(not seriously so -- it just calls diagonal elements twice rather than once). And I never understood what a characteristic of a ring is; chances are this is another thing not consistently understood in Sage. What we probably cannot do is ask whether the characteristic is not 2; if anything, we should ask for 2 to be invertible. But even then, I fear that excepting the NotImplementedError
errors will ruin the speed benefits we get from unifying the code.
Thanks, Travis and Nils, for the review and the helpful comments!
comment:21 Changed 9 years ago by
Replying to darij:
Not sure about it.
is_skewsymmetric
should be a tad slower thanis_alternating
(not seriously so -- it just calls diagonal elements twice rather than once).
If that matters people can always implement both.
What we probably cannot do is ask whether the characteristic is not 2;
uh ...
sage: GF(2).characteristic() !=2 False sage: ZZ.characteristic() != 2 True
if anything, we should ask for 2 to be invertible.
No, that's not the correct check. In ZZ and (ZZ/6ZZ), 2 is not invertible and yet skew symmetric is the same as alternating.
You could check whether 1+1==0, but that's more expensive.
comment:22 Changed 9 years ago by
I think what we really should be checking is if it has positive even characteristic (see the ZZ/4 example in the patch). For example doing something like
def is_alternating(self): if not self.is_skew_symmetric(): return False c = self.base_ring().characteristic() return c != 0 or c % 2 != 0 \ # If past here, we have pos. even char. or all(self.get_unsafe(i,i) == 0 for i in range(self._ncols))
I do agree that is_skew_symmetric()
should be the "main" method, also because it is a larger class of matrices in positive even characteristic. However, I don't see a way to speed up is_skew_symmetric()
for the mod case, so I'm happy with the two ([more] optimized) implementations.
comment:23 follow-up: 24 Changed 9 years ago by
Skew symmetric is NOT the same as alternating in ZZ/(6 ZZ); think of a 3 on the main diagonal. Just having 1 + 1 != 0 is not enough. Infinite characteristic in the sense of "ZZ embeds into the ring" is not enough since we can have things like ZZ[X] / (2X).
comment:24 Changed 9 years ago by
Replying to darij:
Skew symmetric is NOT the same as alternating in ZZ/(6 ZZ); think of a 3 on the main diagonal. Just having 1 + 1 != 0 is not enough. Infinite characteristic in the sense of "ZZ embeds into the ring" is not enough since we can have things like ZZ[X] / (2X).
Oops, sorry. I should use non-integral domains a little more. I withdraw my proposal then. Deciding whether alternating is the same as skew symmetric is not so easy after all, so we should just have two separate methods and if people want the tests to be quicker they will have to implement both themselves.
comment:25 Changed 9 years ago by
Merged in: | → sage-5.13.beta1 |
---|---|
Resolution: | → fixed |
Status: | positive_review → closed |
Hey Darij,
A few things I'd like to see changes:
self.is_skew_symmetric()
instead of implementing your own test?3x3
in the examples block, could you put it in full latex`3 \times 3`
?algorithm
input, I'd rather see something like or a variant with a better name for the algorithm.Also since this file is so big and often subject to changes, it would be better to not include as many whitespace changes.
Thanks,
Travis