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

[sagetrac-git]

Changed commit from 468f238 to 185e915

[sagetrac-git]

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?