Opened 9 years ago

Closed 9 years ago

# Pfaffian of a skew-symmetric matrix

Reported by: Owned by: Darij Grinberg major sage-5.13 combinatorics matrix, sage-combinat, pfaffian Sage Combinat CC user, spancratz, Thierry Monteil sage-5.13.beta1 Darij Grinberg Travis Scrimshaw N/A #14117

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:

### comment:1 Changed 9 years ago by Darij Grinberg

Status: new → needs_review

### comment:2 Changed 9 years ago by Darij Grinberg

Dependencies: → #14117

### comment:3 Changed 9 years ago by Travis Scrimshaw

Hey Darij,

A few things I'd like to see changed:

• Could you call self.is_skew_symmetric() instead of implementing your own test?
• Instead of doing things such as 3x3 in the examples block, could you put it in full latex 3 \times 3?
• For the algorithm input, I'd rather see something like
- algorithm -- string, the algorithm to use; currently the following
algorithms have been implemented:

* 'definition' - using the definition given by perfect matchings

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

Last edited 9 years ago by Travis Scrimshaw (previous) (diff)

### comment:4 follow-up:  7 Changed 9 years ago by Darij Grinberg

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 Darij Grinberg

Description: modified (diff) needs_review → positive_review

### comment:6 Changed 9 years ago by Darij Grinberg

for the patchbot:

apply trac_15245-pfaffian-dg.patch​

### comment:7 in reply to:  4 Changed 9 years ago by Nils Bruin

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 AT=-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 Darij Grinberg

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 Jeroen Demeyer

Status: positive_review → needs_work

### comment:10 Changed 9 years ago by Jeroen Demeyer

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 Darij Grinberg

Oops -- I just realized what Travis gave were comments, not a review. Sorry!

### comment:12 Changed 9 years ago by Travis Scrimshaw

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 Darij Grinberg

I've given the is_skew_symmetric method an additional keyword variable now. All the rest is OK?

### comment:14 in reply to:  13 Changed 9 years ago by Nils Bruin

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 Darij Grinberg

new version, separating skew-symmetry from alternatingness systematically

### comment:15 Changed 9 years ago by Darij Grinberg

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 Travis Scrimshaw

I'm happy with it (having an is_alternating() method). Nils?

### comment:17 in reply to:  16 Changed 9 years ago by Nils Bruin

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 Travis Scrimshaw

Reviewers: → Travis Scrimshaw needs_review → positive_review

Then it's a positive review.

### comment:19 Changed 9 years ago by Nils Bruin

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 Darij Grinberg

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.

Last edited 9 years ago by Darij Grinberg (previous) (diff)

### comment:21 in reply to:  20 Changed 9 years ago by Nils Bruin

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).

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 Travis Scrimshaw

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 Darij Grinberg

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).