Pfaffian of a skew-symmetric matrix

[darijgr]

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:

• attachment: trac_15245-pfaffian-dg.patch

Depends on #14117

CC: @sagetrac-sage-combinat @sagetrac-spancratz @sagetrac-tmonteil

Component: combinatorics

Keywords: matrix, sage-combinat, pfaffian

Author: Darij Grinberg

Reviewer: Travis Scrimshaw

Merged: sage-5.13.beta1

Issue created by migration from https://trac.sagemath.org/ticket/15245

[darijgr]

Dependencies: #14117

[tscrim]

comment:3

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

[darijgr]

comment:4

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

[darijgr]

Description changed:

--- 
+++ 
@@ -3,3 +3,7 @@
 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:
+
+* [attachment: trac_15245-pfaffian-dg.patch](https://github.com/sagemath/sage-prod/files/10658462/trac_15245-pfaffian-dg.patch.gz)

[darijgr]

comment:6

for the patchbot:

apply trac_15245-pfaffian-dg.patch​

[nbruin]

comment:7

Replying to @darijgr:

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.

[darijgr]

comment:8

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?

[jdemeyer]

comment:10

darij: it looks like you reviewed your own patch. I'd say this still needs a "formal" review.

[darijgr]

comment:11

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

[tscrim]

comment:12

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

[darijgr]

comment:13

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

[nbruin]

comment:14

Replying to @darijgr:

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 we mean ".

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

[darijgr]

new version, separating skew-symmetry from alternatingness systematically

[darijgr]

comment:15

Attachment: trac_15245-pfaffian-dg.patch.gz

Done. A number of people seem to be lax about characteristic 2, among them Knuth from whom I had expected this the least.

[tscrim]

comment:16

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

[nbruin]

comment:17

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.

[tscrim]

comment:18

Then it's a positive review.

[tscrim]

Reviewer: Travis Scrimshaw

[nbruin]

comment:19

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>

[darijgr]

comment:20

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!

[nbruin]

comment:21

Replying to @darijgr:

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.

[tscrim]

comment:22

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.

[darijgr]

comment:23

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

[nbruin]

comment:24

Replying to @darijgr:

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.

[jdemeyer]

Merged: sage-5.13.beta1