Fast Pfaffian using Faddeev-LeVerrier

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

CC: @tscrim @mkoeppe @fchapoton @videlec

Component: linear algebra

Keywords: pfaffian

Author: Michael Jung

Branch/Commit: a3ed6c3

Reviewer: Travis Scrimshaw

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

[mjungmath]

Description changed:

--- 
+++ 
@@ -1,5 +1,5 @@
 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(n^4). Furthermore, it is easy to implement. This algorithm works for any torsion-free ring.
+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.

 I suppose there is no way to check the torsion of rings in Sage right now. Hence, I would suggest to restrict the algorithm to integral domains of characteristic zero by performing the computation in the fraction field if necessary and then converting the result back. But I am not certain whether division by integers always turns out working as intended in the fraction field.

[mjungmath]

Branch: u/gh-mjungmath/baer_faddeev_leverrier

[mjungmath]

Commit: 468f238

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 468f238 to 185e915

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

185e915

Trac #30681: fast pfaffian

[mjungmath]

Description changed:

--- 
+++ 
@@ -1,5 +1,3 @@
 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.
-
-I suppose there is no way to check the torsion of rings in Sage right now. Hence, I would suggest to restrict the algorithm to integral domains of characteristic zero by performing the computation in the fraction field if necessary and then converting the result back. But I am not certain whether division by integers always turns out working as intended in the fraction field.

[tscrim]

comment:8

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. :/

[mjungmath]

comment:9

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?

[mjungmath]

comment:10

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.

[mjungmath]

comment:11

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?

[tscrim]

comment:12

Replying to @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.)

[tscrim]

comment:13

Replying to @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.

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

d41151f

Trac #30681: first try BFL algorithm + Cythonized

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 185e915 to d41151f

[mjungmath]

comment:15

You mean like that?

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

7a6b9b1

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

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from d41151f to 7a6b9b1

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

dc44a5a

Trac #30681: doc improved

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 7a6b9b1 to dc44a5a

[tscrim]

comment:18

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.

[mjungmath]

comment:19

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?

[tscrim]

comment:20

Replying to @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.

[mjungmath]

comment:21

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.

[mjungmath]

comment:22

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.

[tscrim]

comment:23

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.

[tscrim]

comment:24

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.

[mjungmath]

comment:25

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?

[tscrim]

comment:26

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

[mjungmath]

comment:27

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.

[tscrim]

comment:28

Replying to @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.

[mjungmath]

comment:29

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.

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

81a1234

Trac #30681: fast pfaffian

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from dc44a5a to 81a1234

[mjungmath]

comment:31

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.

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 81a1234 to e1cc07a

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

e1cc07a

Trac #30681: improved doc

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

ced33c9

Trac #30681: fast pfaffian

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from e1cc07a to ced33c9

[mjungmath]

comment:34

Squashed to one commit. The changes were negligible.

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

a8cceac

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

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from ced33c9 to a8cceac

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

061df41

Trac #30681: base_ring check not robust

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from a8cceac to 061df41

[mjungmath]

comment:39

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.

[tscrim]

comment:40

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.

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 061df41 to 6247a75

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

6247a75

Trac #30681: minor details

[mjungmath]

comment:42

Done. Thank you. :)

[tscrim]

Reviewer: Travis Scrimshaw

[tscrim]

comment:43

Thank you.