Opened 8 years ago

Closed 8 years ago

#15245 closed enhancement (fixed)

Pfaffian of a skew-symmetric matrix

Reported by: darij Owned by:
Priority: major Milestone: sage-5.13
Component: combinatorics Keywords: matrix, sage-combinat, pfaffian
Cc: sage-combinat, spancratz, tmonteil Merged in: sage-5.13.beta1
Authors: Darij Grinberg Reviewers: Travis Scrimshaw
Report Upstream: N/A Work issues:
Branch: Commit:
Dependencies: #14117 Stopgaps:

Status badges

Description (last modified by darij)

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)

trac_15245-pfaffian-dg.patch (25.7 KB) - added by darij 8 years ago.
new version, separating skew-symmetry from alternatingness systematically

Download all attachments as: .zip

Change History (26)

comment:1 Changed 8 years ago by darij

  • Status changed from new to needs_review

comment:2 Changed 8 years ago by darij

  • Dependencies set to #14117

comment:3 Changed 8 years ago by tscrim

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 8 years ago by tscrim (previous) (diff)

comment:4 follow-up: Changed 8 years ago by darij

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 8 years ago by darij

  • Description modified (diff)
  • Status changed from needs_review to positive_review

comment:6 Changed 8 years ago by darij

for the patchbot:

apply trac_15245-pfaffian-dg.patch​

comment:7 in reply to: ↑ 4 Changed 8 years ago by nbruin

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 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 8 years ago by darij

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 8 years ago by jdemeyer

  • Status changed from positive_review to needs_work

comment:10 Changed 8 years ago by jdemeyer

  • Status changed from needs_work to needs_review

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

comment:11 Changed 8 years ago by darij

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

comment:12 Changed 8 years ago by tscrim

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: Changed 8 years ago by darij

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 8 years ago by nbruin

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 8 years ago by darij

new version, separating skew-symmetry from alternatingness systematically

comment:15 Changed 8 years ago by darij

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: Changed 8 years ago by tscrim

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

comment:17 in reply to: ↑ 16 Changed 8 years ago by nbruin

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 8 years ago by tscrim

  • Reviewers set to Travis Scrimshaw
  • Status changed from needs_review to positive_review

Then it's a positive review.

comment:19 Changed 8 years ago by nbruin

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: Changed 8 years ago by darij

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!

Last edited 8 years ago by darij (previous) (diff)

comment:21 in reply to: ↑ 20 Changed 8 years ago by nbruin

Replying to darij:

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 8 years ago by tscrim

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: Changed 8 years ago by 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).

comment:24 in reply to: ↑ 23 Changed 8 years ago by nbruin

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 8 years ago by jdemeyer

  • Merged in set to sage-5.13.beta1
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.