#30681 closed enhancement (fixed)

Fast Pfaffian using Faddeev-LeVerrier

Reported by: gh-mjungmath Owned by:
Priority: major Milestone: sage-9.3
Component: linear algebra Keywords: pfaffian
Cc: tscrim, mkoeppe, chapoton, vdelecroix Merged in:
Authors: Michael Jung Reviewers: Travis Scrimshaw
Report Upstream: N/A Work issues:
Branch: a3ed6c3 (Commits, GitHub, GitLab) Commit: a3ed6c3132c760ca967a55f4f10a925f8be06838
Dependencies: Stopgaps:

Status badges

Description (last modified by gh-mjungmath)

At the current stage, the Pfaffian is computed using the definition by perfect matchings. This is tremendously demanding.

According to https://arxiv.org/abs/2008.04247, there is an algorithm similar to the Faddeev-LeVerrier? algorithm for the determinant running in at most O(pow(n,4)). Furthermore, it is easy to implement. This algorithm works for any torsion-free ring.

Change History (58)

comment:1 Changed 15 months ago by gh-mjungmath

  • Description modified (diff)

comment:2 Changed 15 months ago by gh-mjungmath

  • Branch set to u/gh-mjungmath/baer_faddeev_leverrier

comment:3 Changed 15 months ago by gh-mjungmath

  • Commit set to 468f23815ade42a2192b0a9cd378de8fdc594dcd
  • Status changed from new to needs_review

comment:4 Changed 15 months ago by git

  • Commit changed from 468f23815ade42a2192b0a9cd378de8fdc594dcd to 185e9158e3ff09c3944ed47e26a1ef05054c1599

Branch pushed to git repo; I updated commit sha1. New commits:

185e915Trac #30681: fast pfaffian

comment:5 Changed 15 months ago by gh-mjungmath

Last edited 15 months ago by gh-mjungmath (previous) (diff)

comment:6 Changed 15 months ago by gh-mjungmath

  • Status changed from needs_review to needs_info

comment:7 Changed 15 months ago by gh-mjungmath

  • Description modified (diff)

comment:8 follow-up: Changed 15 months ago by tscrim

What do you need info on? The determinant version? I think that would be better as a separate ticket.

I don't see the point in nailing IntegralDomains() in memory here. This function should not be called with the kind of frequency that would mandate that.

When the result is not a field, how does other such algorithms work in Sage? Do they convert the result back to the original ring? I feel like if I use a matrix over ZZ, I would expect a result in ZZ.

I would fully Cythonize _pf_bfl and write it as

    cdef _pf_bfl(self):
        cdef Py_ssize_t n = self._ncols
        cdef Py_ssize_t q = n // 2
        cdef Py_ssize_t i, k

        # apply J:
        cdef Matrix A = <Matrix> copy(self)
        for i in range(0, n, 2):
            A.swap_columns_c(i, i+1)  # Avoid some checks
            # A.set_col_to_multiple_of_col(i+1, i+1, -1)
            for k in range(n):
                A.set_unsafe(k, i+1, -A.get_unsage(k, i+1))

        cdef Matrix M = <Matrix> copy(A)

        # Baer-Faddeev-Leverrier algorithm:
        for k in range(1, q):
            c = -M.trace() / (2*k)
            # M = A * (M + c*1)
            # Add c along the diagonal
            for i in range(n):
                M.set_unsafe(i, i, M.get_unsafe(i, i) + c)
            M = A * M
        c = -M.trace() / (2*q)
        return (-1)**q * c

It would be really good if there was an in place version of M = A * M, but alas, I don't think there is. So we will have to have that temporary object. :/

comment:9 follow-up: Changed 14 months ago by gh-mjungmath

Thank you for the feedback. This Cython version is already much better.

One question that occurred to me is whether the division by integers always works in the (fraction) field. It should, but I am not sure whether Sage always converts the integers mathematically correct.

Secondly, the algorithm should work for QQ algebras as well. I would guess, you can even assume an algebra over an integral domain of characteristic zero. Can you confirm that? But even if that works, it also depends on the answer of my first question.

Obviously, each ring constitutes an algebra over itself. But apparently, Sage doesn't recognize that. So, if the answer to all my above questions is "yes", how can we establish a general check for the matrix's base ring being an algebra over an integral domain with characteristic zero? And how would a conversion look like?

Last edited 14 months ago by gh-mjungmath (previous) (diff)

comment:10 in reply to: ↑ 8 Changed 14 months ago by gh-mjungmath

Replying to tscrim:

When the result is not a field, how does other such algorithms work in Sage? Do they convert the result back to the original ring? I feel like if I use a matrix over ZZ, I would expect a result in ZZ.

Yes, it should. At least for ZZ, the current version does the job, in particular:

+                if R in _Fields:
+                    pf = self._pf_bfl()
+                else:
+                    F = R.fraction_field()
+                    pf = self._coerce_element(self.change_ring(F)._pf_bfl())

But I am not sure whether this setup always works out.

comment:11 follow-up: Changed 14 months ago by gh-mjungmath

The _pf_bfl method is invoked from a regular Python function. I am not a Cython expert, but shouldn't that be a cpdef function then?

Or do you suggest to rewrite pfaffian as well into Cython code?

comment:12 in reply to: ↑ 11 Changed 14 months ago by tscrim

Replying to gh-mjungmath:

The _pf_bfl method is invoked from a regular Python function. I am not a Cython expert, but shouldn't that be a cpdef function then?

No, because it is called within Cython code. You might need to change

-self.change_ring(F)._pf_bfl()
+(<Matrix> self.change_ring(F))._pf_bfl()

in the appropriate place so Cython knows the method exists.

Or do you suggest to rewrite pfaffian as well into Cython code?

No, that isn't worth it I think because the heavy lifting is done in the _pf_bfl method, which the end user will not see. (Although that is not to say that you shouldn't put some typing on things where appropriate.)

comment:13 in reply to: ↑ 9 Changed 14 months ago by tscrim

Replying to gh-mjungmath:

One question that occurred to me is whether the division by integers always works in the (fraction) field. It should, but I am not sure whether Sage always converts the integers mathematically correct.

When one of them is a Sage integer/rational/etc., then yes. If they are Python ints then they are handled like Python ints IIRC (although possibly like C ints). However, in response to your more general question, I believe the answer is yes.

Secondly, the algorithm should work for QQ algebras as well. I would guess, you can even assume an algebra over an integral domain of characteristic zero. Can you confirm that? But even if that works, it also depends on the answer of my first question.

I am not sure what you are asking me here. I think the answer is yes, probably, but I am not sure.

Obviously, each ring constitutes an algebra over itself. But apparently, Sage doesn't recognize that. So, if the answer to all my above questions is "yes", how can we establish a general check for the matrix's base ring being an algebra over an integral domain with characteristic zero? And how would a conversion look like?

Well, one option is to just try and and let it fail if it does something invalid. It just needs to be documented with what it is guaranteed to work for.

comment:14 Changed 14 months ago by git

  • Commit changed from 185e9158e3ff09c3944ed47e26a1ef05054c1599 to d41151f55294854879964cfb5b8f03eff4ed4266

Branch pushed to git repo; I updated commit sha1. New commits:

d41151fTrac #30681: first try BFL algorithm + Cythonized

comment:15 Changed 14 months ago by gh-mjungmath

  • Status changed from needs_info to needs_review

You mean like that?

comment:16 Changed 14 months ago by git

  • Commit changed from d41151f55294854879964cfb5b8f03eff4ed4266 to 7a6b9b1af3c55b2a43794415f4beb99e1197a324

Branch pushed to git repo; I updated commit sha1. New commits:

7a6b9b1Trac #30681: verbose if BFL fails + doc improved + import reverted

comment:17 Changed 14 months ago by git

  • Commit changed from 7a6b9b1af3c55b2a43794415f4beb99e1197a324 to dc44a5adf4a96c3fe53b1a1816b8668003a6bb3a

Branch pushed to git repo; I updated commit sha1. New commits:

dc44a5aTrac #30681: doc improved

comment:18 follow-up: Changed 14 months ago by tscrim

The / division will almost always convert to the fraction field, so there is no point in not first converting the initial matrix (in fact, this might even cause a crash). The other option would be doing // division, which then might return wrong results if the conditions set in the .. WARNING:: are not satisfied.

comment:19 in reply to: ↑ 18 ; follow-up: Changed 14 months ago by gh-mjungmath

Replying to tscrim:

The / division will almost always convert to the fraction field, so there is no point in not first converting the initial matrix (in fact, this might even cause a crash).

That was more or less my first thought, too. But think about examples like polynomial rings over ZZ where the element can be coerced into an element of QQ[x]. In that case, the fraction field setup is not necessary. In fact, if you already start with an QQ algebra that even has no fraction field implemented, the conversion to a fraction field will fail even though the algorithm would work. Thus, I would avoid going into the fraction field prior because it is not always the suitable choice.

Do you have an example where no conversion would cause a crash?

comment:20 in reply to: ↑ 19 Changed 14 months ago by tscrim

Replying to gh-mjungmath:

Replying to tscrim:

The / division will almost always convert to the fraction field, so there is no point in not first converting the initial matrix (in fact, this might even cause a crash).

That was more or less my first thought, too. But think about examples like polynomial rings over ZZ where the element can be coerced into an element of QQ[x]. In that case, the fraction field setup is not necessary. In fact, if you already start with an QQ algebra that even has no fraction field implemented, the conversion to a fraction field will fail even though the algorithm would work. Thus, I would avoid going into the fraction field prior because it is not always the suitable choice.

Division could be a subtle thing in Sage whether it computes the pushout first or the fraction field first. This might warrant some testing.

Do you have an example where no conversion would cause a crash?

Since it doesn't change the type of the matrix, the entries might be expected to be Integer but instead be set to Rational. The type mismatch would then lead to a crash.

comment:21 Changed 14 months ago by gh-mjungmath

Okay, what do you suggest? Something like that perhaps?

try:
    elt = R.zero() / 1 # try division by integers
    F = elt.parent() # get parent
    R(elt) # try conversion back
    temp = <Matrix> self.change_ring(F)
    ... 

If that succeeds, the matrix should be applicable for the BFL algorithm.

Last edited 14 months ago by gh-mjungmath (previous) (diff)

comment:22 Changed 14 months ago by gh-mjungmath

Different approach: the paper explains that you can simply go into the rationalization. Is Sage capable of that? I didn't find anything useful so far, but maybe I haven't looked closely enough.

comment:23 follow-up: Changed 14 months ago by tscrim

I am not quite sure what you mean by the rationalization. I am guessing you don't mean the fraction field.

Actually, one thing you could do is

        # Baer-Faddeev-Leverrier algorithm:
        R = self._parent._base
        for k in range(1, q):
            c = R(-M.trace() / (2*k))
            # add c along the diagonal
            for i in range(n):
                M.set_unsafe(i ,i, M.get_unsafe(i, i) + c)
            M = A * M
        c = R(-M.trace() / (2*q))
        return (-1)**q * c

It will be slightly slower for matrices over fields, but my guess is it wouldn't even be significant except in some possibly some really small examples. For non-fields, it would guarantee every element belongs to R after every division. However, this would then fail for Pfaffians over, e.g., ZZ because an intermediate computation needed to be in QQ.

I think the best approach is probably converting the matrix to the fraction field, then coercing the result back to R.

Last edited 14 months ago by tscrim (previous) (diff)

comment:24 in reply to: ↑ 23 Changed 14 months ago by tscrim

Replying to tscrim:

I am not quite sure what you mean by the rationalization. I am guessing you don't mean the fraction field.

Ah, perhaps you mean extension of scalars to be a QQ-algebra. If this is the case, then Sage doesn't have this feature yet. You can change rings to be over QQ, but that won't work for algebras over, e.g., ZZ['t']. Well, I guess perhaps you could keep doing this until the base ring is itself:

sage: ZZ.base_ring()                                                                                                  
Integer Ring

and then build new rings back. Although that feels a bit overkill right now as we can just go to the fraction field for this time with perhaps a TODO note.

comment:25 Changed 14 months ago by gh-mjungmath

Fraction fields and TODO, check. But what about QQ algebras in particular? This would be very beneficial for e.g. mixed forms.

Besides, what do you think about my proposal in comment 21?

Last edited 14 months ago by gh-mjungmath (previous) (diff)

comment:26 follow-up: Changed 14 months ago by tscrim

My first thought for QQ algebras would be to test the category, but then I am not so convinced this is a good approach. Your idea in comment:21 should work unless an implementation does something very special for zero. I think you are better off doing

try:
    F = parent(R.one() / 2)
    R(F.one())  # trivial check that we can go back
    temp = <Matrix> self.change_ring(F)
    ...

since that way it checks something non-trivial and that division is something that will actually be performed (so it will more quickly fail if you are working over GF(2)).

comment:27 in reply to: ↑ 26 ; follow-up: Changed 14 months ago by gh-mjungmath

Replying to tscrim:

My first thought for QQ algebras would be to test the category, but then I am not so convinced this is a good approach.

Why do you think this is not a good approach?

What about:

F = R.base_ring()
if QQ.is_subring(F):
    pf = self._pf_bfl()

This should work always.

Your idea in comment:21 should work unless an implementation does something very special for zero. I think you are better off doing

try:
    F = parent(R.one() / 2)
    R(F.one())  # trivial check that we can go back
    temp = <Matrix> self.change_ring(F)
    ...

since that way it checks something non-trivial and that division is something that will actually be performed (so it will more quickly fail if you are working over GF(2)).

I have a bad feeling about this approach. There is still plenty of room what can go wrong here, isn't it? I peeked a bit into the matrix code; take a closer look at with_added_multiple_of_row in matrix0.pyx, there you want to multiply a row by an integer. Perhaps this is something we can adapt:

F = Sequence([R.one()] + [1/(2*j) for j in range(1, k//2 + 1)]).universe()
if F in Rings():
    temp = <Matrix> self.change_ring(F)
    pf = R(temp._pf_bfl())

I am not completely sure how and whether this works.

Last edited 14 months ago by gh-mjungmath (previous) (diff)

comment:28 in reply to: ↑ 27 ; follow-up: Changed 14 months ago by tscrim

Replying to gh-mjungmath:

Replying to tscrim:

My first thought for QQ algebras would be to test the category, but then I am not so convinced this is a good approach.

Why do you think this is not a good approach?

Because you could have something over ZZ['q'].fraction_field(), or worse over QQ.category() (this was done to minimize the number of categories constructed when varying over prime fields).

What about:

F = R.base_ring()
if QQ.is_subring(F):
    pf = self._pf_bfl()

This should work always.

I am not entirely sure how much I trust that. It is more robust than I expected (I actually expected QQ.is_subring(ZZ['t'].fraction_field()) to fail). From this check, I am comfortable enough to use this.

Your idea in comment:21 should work unless an implementation does something very special for zero. I think you are better off doing

try:
    F = parent(R.one() / 2)
    R(F.one())  # trivial check that we can go back
    temp = <Matrix> self.change_ring(F)
    ...

since that way it checks something non-trivial and that division is something that will actually be performed (so it will more quickly fail if you are working over GF(2)).

I have a bad feeling about this approach. There is still plenty of room what can go wrong here, isn't it?

It would fail later on, but that is okay.

I peeked a bit into the matrix code; take a closer look at with_added_multiple_of_row in matrix0.pyx, there you want to multiply a row by an integer. Perhaps this is something we can adapt:

F = Sequence([R.one()] + [1/(2*j) for j in range(1, k//2 + 1)]).universe()
temp = <Matrix> self.change_ring(F)
pf = R(temp._pf_bfl())

There shouldn't be any difference between dividing by 2 or 4 or 6 (unless working over GF(3), but like I said above, this would fail later on, which is okay with me). Perhaps the QQ.is_subring() approach is best. At a certain point, we might be putting too much thought into this.

comment:29 in reply to: ↑ 28 Changed 14 months ago by gh-mjungmath

Replying to tscrim:

There shouldn't be any difference between dividing by 2 or 4 or 6 (unless working over GF(3), but like I said above, this would fail later on, which is okay with me). Perhaps the QQ.is_subring() approach is best. At a certain point, we might be putting too much thought into this.

Exactly. The algorithm should also work if the ring's characteristic is bigger than the matrix dimension. This is also guaranteed by this approach. But perhaps you are right, we are putting too much thought into this.

I think, this might be good:

if R in IntegralDomains():
    F = R.fraction_field()
    temp = <Matrix> self.change_ring(F)
else:
    F = R.base_ring()
    temp = <Matrix> copy(self)
if QQ.is_subring(F):
    pf = R(temp._pf_bfl())

This checks already make sure that the characteristic is zero even without using characteristic (which is good because the mixed form algebra for example does not have this method). Moreover, this guarantees the integer division to work.

Last edited 14 months ago by gh-mjungmath (previous) (diff)

comment:30 Changed 14 months ago by git

  • Commit changed from dc44a5adf4a96c3fe53b1a1816b8668003a6bb3a to 81a12346099c1a1d1bdf1970b0dd522a6e7cb05e

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

81a1234Trac #30681: fast pfaffian

comment:31 Changed 14 months ago by gh-mjungmath

Check this out. This should be fine. The QQ sub-ring check is good enough to guarantee a working algorithm. And the number of applicable rings is reasonable. Ready for (final) review.

Last edited 14 months ago by gh-mjungmath (previous) (diff)

comment:32 Changed 14 months ago by git

  • Commit changed from 81a12346099c1a1d1bdf1970b0dd522a6e7cb05e to e1cc07a62988aef49eaffa7d6dad261b263fdb10

Branch pushed to git repo; I updated commit sha1. New commits:

e1cc07aTrac #30681: improved doc

comment:33 Changed 14 months ago by git

  • Commit changed from e1cc07a62988aef49eaffa7d6dad261b263fdb10 to ced33c96adcd60d9bbbfe899368e8bb5da297db9

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

ced33c9Trac #30681: fast pfaffian

comment:34 Changed 14 months ago by gh-mjungmath

Squashed to one commit. The changes were negligible.

comment:35 Changed 14 months ago by gh-mjungmath

  • Status changed from needs_review to needs_work

comment:36 Changed 14 months ago by git

  • Commit changed from ced33c96adcd60d9bbbfe899368e8bb5da297db9 to a8cceac6028f2f71d136505ca0521bff59bf87e0

Branch pushed to git repo; I updated commit sha1. New commits:

a8cceacTrac #30681: base_ring check not enough for fraction fields, ZZ[x] as test added

comment:37 Changed 14 months ago by gh-mjungmath

  • Status changed from needs_work to needs_review

comment:38 Changed 14 months ago by git

  • Commit changed from a8cceac6028f2f71d136505ca0521bff59bf87e0 to 061df4147ff2504d797f6dd6bef90c4183078b86

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

061df41Trac #30681: base_ring check not robust

comment:39 Changed 14 months ago by gh-mjungmath

Turns out that is_subring is indeed very robust. Works for algebras such as mixed forms and scalar fields, too. Using base_ring was not a good idea, fails for e.g. fraction fields.

Sorry for the mess. I hope that should be it. If I should add any more tests/examples, let me know.

Last edited 14 months ago by gh-mjungmath (previous) (diff)

comment:40 Changed 14 months ago by tscrim

Minor details:

-        Computes the Pfaffian of ``self`` using the Baer-Faddeev-LeVerrier
-        algorithm. This method assumes that the base ring is an `\QQ`-algebra.
+        Computes the Pfaffian of ``self`` using the Baer-Faddeev-LeVerrier
+        algorithm.
+
+        .. WARNING::
+
+            This method assumes that the base ring is an `\QQ`-algebra.

Also Computes -> Compute for _pf_perfect_matchings and bad comma spacing:

-M.set_unsafe(i ,i, M.get_unsafe(i, i) + c)
+M.set_unsafe(i, i, M.get_unsafe(i, i) + c)

I would also put the bfl method as the first choice in the doc since it is the "default" choice.

comment:41 Changed 14 months ago by git

  • Commit changed from 061df4147ff2504d797f6dd6bef90c4183078b86 to 6247a7527d8428f6914c723ff9e57f6908818814

Branch pushed to git repo; I updated commit sha1. New commits:

6247a75Trac #30681: minor details

comment:42 Changed 14 months ago by gh-mjungmath

Done. Thank you. :)

Last edited 14 months ago by gh-mjungmath (previous) (diff)

comment:43 Changed 14 months ago by tscrim

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

Thank you.

comment:44 Changed 14 months ago by gh-mjungmath

Thank you for the review. :)

comment:45 Changed 14 months ago by git

  • Commit changed from 6247a7527d8428f6914c723ff9e57f6908818814 to d28f8f8bfa237b43c3b08695b498e4daa5e312c1
  • Status changed from positive_review to needs_review

Branch pushed to git repo; I updated commit sha1 and set ticket back to needs_review. New commits:

d28f8f8Trac #30681: authorship and date added

comment:46 Changed 14 months ago by git

  • Commit changed from d28f8f8bfa237b43c3b08695b498e4daa5e312c1 to 7211239193f0bc03a32c2d29d5d2878de380c7c1

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

3ffff71# Das ist eine Kombination aus 3 Commits.
7211239Trac #30681: fast pfaffian

comment:47 Changed 14 months ago by git

  • Commit changed from 7211239193f0bc03a32c2d29d5d2878de380c7c1 to 5bb279286d8ad62e1804a39bb17ad475d1b71d33

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

5bb2792Trac #30681: fast pfaffian

comment:48 Changed 14 months ago by gh-mjungmath

  • Status changed from needs_review to positive_review

Added authorship and date and screwed up the merge.

comment:49 Changed 14 months ago by git

  • Commit changed from 5bb279286d8ad62e1804a39bb17ad475d1b71d33 to b3b2a83fa2cbb217c3aa09585c66aad2a8809945
  • Status changed from positive_review to needs_review

Branch pushed to git repo; I updated commit sha1 and set ticket back to needs_review. New commits:

b3b2a83Trac #30681: cache functionality ensured

comment:50 Changed 14 months ago by gh-mjungmath

I just noticed that I forgot to fetch the cached result. Commit is pushed.

comment:51 Changed 14 months ago by tscrim

Yes, that is good to use, although I think the cache check should be first. Once that is moved, you can set a positive review on my behalf.

comment:52 Changed 14 months ago by git

  • Commit changed from b3b2a83fa2cbb217c3aa09585c66aad2a8809945 to 1fc9e3b543916d5a4360b3ed3b0d22e6fbd0b63b

Branch pushed to git repo; I updated commit sha1. New commits:

1fc9e3bTrac #30681: fetch cached result first

comment:53 Changed 14 months ago by gh-mjungmath

  • Status changed from needs_review to positive_review

Thanks.

comment:54 Changed 14 months ago by git

  • Commit changed from 1fc9e3b543916d5a4360b3ed3b0d22e6fbd0b63b to 35ca1bd77a9d74aa2bf5efc37efb5a5ac1996e68
  • Status changed from positive_review to needs_review

Branch pushed to git repo; I updated commit sha1 and set ticket back to needs_review. New commits:

35ca1bdTrac #30681: tautological test replaced due to caching

comment:55 Changed 14 months ago by git

  • Commit changed from 35ca1bd77a9d74aa2bf5efc37efb5a5ac1996e68 to a3ed6c3132c760ca967a55f4f10a925f8be06838

Branch pushed to git repo; I updated commit sha1. New commits:

a3ed6c3Trac 30681: comment removed

comment:56 Changed 14 months ago by gh-mjungmath

  • Status changed from needs_review to positive_review

comment:57 Changed 14 months ago by mkoeppe

  • Milestone changed from sage-9.2 to sage-9.3

comment:58 Changed 13 months ago by vbraun

  • Branch changed from u/gh-mjungmath/baer_faddeev_leverrier to a3ed6c3132c760ca967a55f4f10a925f8be06838
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.