Opened 2 years ago
Closed 2 years ago
#30681 closed enhancement (fixed)
Fast Pfaffian using FaddeevLeVerrier
Reported by:  Michael Jung  Owned by:  

Priority:  major  Milestone:  sage9.3 
Component:  linear algebra  Keywords:  pfaffian 
Cc:  Travis Scrimshaw, Matthias Köppe, Frédéric Chapoton, Vincent Delecroix  Merged in:  
Authors:  Michael Jung  Reviewers:  Travis Scrimshaw 
Report Upstream:  N/A  Work issues:  
Branch:  a3ed6c3 (Commits, GitHub, GitLab)  Commit:  a3ed6c3132c760ca967a55f4f10a925f8be06838 
Dependencies:  Stopgaps: 
Description (last modified by )
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 FaddeevLeVerrier? algorithm for the determinant running in at most O(pow(n,4)). Furthermore, it is easy to implement. This algorithm works for any torsionfree ring.
Change History (58)
comment:1 Changed 2 years ago by
Description:  modified (diff) 

comment:2 Changed 2 years ago by
Branch:  → u/ghmjungmath/baer_faddeev_leverrier 

comment:3 Changed 2 years ago by
Commit:  → 468f23815ade42a2192b0a9cd378de8fdc594dcd 

Status:  new → needs_review 
comment:4 Changed 2 years ago by
Commit:  468f23815ade42a2192b0a9cd378de8fdc594dcd → 185e9158e3ff09c3944ed47e26a1ef05054c1599 

comment:6 Changed 2 years ago by
Status:  needs_review → needs_info 

comment:7 Changed 2 years ago by
Description:  modified (diff) 

comment:8 followup: 10 Changed 2 years ago by
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) # BaerFaddeevLeverrier 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 followup: 13 Changed 2 years ago by
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?
comment:10 Changed 2 years ago by
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 inZZ
.
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 followup: 12 Changed 2 years ago by
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 Changed 2 years ago by
Replying to ghmjungmath:
The
_pf_bfl
method is invoked from a regular Python function. I am not a Cython expert, but shouldn't that be acpdef
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 Changed 2 years ago by
Replying to ghmjungmath:
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 int
s then they are handled like Python int
s IIRC (although possibly like C int
s). However, in response to your more general question, I believe the answer is yes.
Secondly, the algorithm should work for
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 2 years ago by
Commit:  185e9158e3ff09c3944ed47e26a1ef05054c1599 → d41151f55294854879964cfb5b8f03eff4ed4266 

Branch pushed to git repo; I updated commit sha1. New commits:
d41151f  Trac #30681: first try BFL algorithm + Cythonized

comment:16 Changed 2 years ago by
Commit:  d41151f55294854879964cfb5b8f03eff4ed4266 → 7a6b9b1af3c55b2a43794415f4beb99e1197a324 

Branch pushed to git repo; I updated commit sha1. New commits:
7a6b9b1  Trac #30681: verbose if BFL fails + doc improved + import reverted

comment:17 Changed 2 years ago by
Commit:  7a6b9b1af3c55b2a43794415f4beb99e1197a324 → dc44a5adf4a96c3fe53b1a1816b8668003a6bb3a 

Branch pushed to git repo; I updated commit sha1. New commits:
dc44a5a  Trac #30681: doc improved

comment:18 followup: 19 Changed 2 years ago by
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 followup: 20 Changed 2 years ago by
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 Changed 2 years ago by
Replying to ghmjungmath:
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 ofQQ[x]
. In that case, the fraction field setup is not necessary. In fact, if you already start with an
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 2 years ago by
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.
comment:22 Changed 2 years ago by
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 followup: 24 Changed 2 years ago by
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
# BaerFaddeevLeverrier 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 nonfields, 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
.
comment:24 Changed 2 years ago by
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 2 years ago by
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?
comment:26 followup: 27 Changed 2 years ago by
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 nontrivial 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 followup: 28 Changed 2 years ago by
Replying to tscrim:
My first thought for
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 nontrivial 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.
comment:28 followup: 29 Changed 2 years ago by
Replying to ghmjungmath:
Replying to tscrim:
My first thought for
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 nontrivial 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
inmatrix0.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 Changed 2 years ago by
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 theQQ.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.
comment:30 Changed 2 years ago by
Commit:  dc44a5adf4a96c3fe53b1a1816b8668003a6bb3a → 81a12346099c1a1d1bdf1970b0dd522a6e7cb05e 

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:
81a1234  Trac #30681: fast pfaffian

comment:31 Changed 2 years ago by
Check this out. This should be fine. The QQ
subring check is good enough to guarantee a working algorithm. And the number of applicable rings is reasonable. Ready for (final) review.
comment:32 Changed 2 years ago by
Commit:  81a12346099c1a1d1bdf1970b0dd522a6e7cb05e → e1cc07a62988aef49eaffa7d6dad261b263fdb10 

Branch pushed to git repo; I updated commit sha1. New commits:
e1cc07a  Trac #30681: improved doc

comment:33 Changed 2 years ago by
Commit:  e1cc07a62988aef49eaffa7d6dad261b263fdb10 → ced33c96adcd60d9bbbfe899368e8bb5da297db9 

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:
ced33c9  Trac #30681: fast pfaffian

comment:35 Changed 2 years ago by
Status:  needs_review → needs_work 

comment:36 Changed 2 years ago by
Commit:  ced33c96adcd60d9bbbfe899368e8bb5da297db9 → a8cceac6028f2f71d136505ca0521bff59bf87e0 

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

comment:37 Changed 2 years ago by
Status:  needs_work → needs_review 

comment:38 Changed 2 years ago by
Commit:  a8cceac6028f2f71d136505ca0521bff59bf87e0 → 061df4147ff2504d797f6dd6bef90c4183078b86 

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:
061df41  Trac #30681: base_ring check not robust

comment:39 Changed 2 years ago by
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.
comment:40 Changed 2 years ago by
Minor details:
 Computes the Pfaffian of ``self`` using the BaerFaddeevLeVerrier  algorithm. This method assumes that the base ring is an `\QQ`algebra. + Computes the Pfaffian of ``self`` using the BaerFaddeevLeVerrier + 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 2 years ago by
Commit:  061df4147ff2504d797f6dd6bef90c4183078b86 → 6247a7527d8428f6914c723ff9e57f6908818814 

Branch pushed to git repo; I updated commit sha1. New commits:
6247a75  Trac #30681: minor details

comment:43 Changed 2 years ago by
Reviewers:  → Travis Scrimshaw 

Status:  needs_review → positive_review 
Thank you.
comment:45 Changed 2 years ago by
Commit:  6247a7527d8428f6914c723ff9e57f6908818814 → d28f8f8bfa237b43c3b08695b498e4daa5e312c1 

Status:  positive_review → needs_review 
Branch pushed to git repo; I updated commit sha1 and set ticket back to needs_review. New commits:
d28f8f8  Trac #30681: authorship and date added

comment:46 Changed 2 years ago by
Commit:  d28f8f8bfa237b43c3b08695b498e4daa5e312c1 → 7211239193f0bc03a32c2d29d5d2878de380c7c1 

comment:47 Changed 2 years ago by
Commit:  7211239193f0bc03a32c2d29d5d2878de380c7c1 → 5bb279286d8ad62e1804a39bb17ad475d1b71d33 

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:
5bb2792  Trac #30681: fast pfaffian

comment:48 Changed 2 years ago by
Status:  needs_review → positive_review 

Added authorship and date and screwed up the merge.
comment:49 Changed 2 years ago by
Commit:  5bb279286d8ad62e1804a39bb17ad475d1b71d33 → b3b2a83fa2cbb217c3aa09585c66aad2a8809945 

Status:  positive_review → needs_review 
Branch pushed to git repo; I updated commit sha1 and set ticket back to needs_review. New commits:
b3b2a83  Trac #30681: cache functionality ensured

comment:50 Changed 2 years ago by
I just noticed that I forgot to fetch the cached result. Commit is pushed.
comment:51 Changed 2 years ago by
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 2 years ago by
Commit:  b3b2a83fa2cbb217c3aa09585c66aad2a8809945 → 1fc9e3b543916d5a4360b3ed3b0d22e6fbd0b63b 

Branch pushed to git repo; I updated commit sha1. New commits:
1fc9e3b  Trac #30681: fetch cached result first

comment:54 Changed 2 years ago by
Commit:  1fc9e3b543916d5a4360b3ed3b0d22e6fbd0b63b → 35ca1bd77a9d74aa2bf5efc37efb5a5ac1996e68 

Status:  positive_review → needs_review 
Branch pushed to git repo; I updated commit sha1 and set ticket back to needs_review. New commits:
35ca1bd  Trac #30681: tautological test replaced due to caching

comment:55 Changed 2 years ago by
Commit:  35ca1bd77a9d74aa2bf5efc37efb5a5ac1996e68 → a3ed6c3132c760ca967a55f4f10a925f8be06838 

Branch pushed to git repo; I updated commit sha1. New commits:
a3ed6c3  Trac 30681: comment removed

comment:56 Changed 2 years ago by
Status:  needs_review → positive_review 

comment:57 Changed 2 years ago by
Milestone:  sage9.2 → sage9.3 

comment:58 Changed 2 years ago by
Branch:  u/ghmjungmath/baer_faddeev_leverrier → a3ed6c3132c760ca967a55f4f10a925f8be06838 

Resolution:  → fixed 
Status:  positive_review → closed 
Branch pushed to git repo; I updated commit sha1. New commits:
Trac #30681: fast pfaffian