5
$\begingroup$

I have the following code, that computes the Pfaffian of an even dimensional anti-symmetric matrix via direct row expansion:

drop[m_, parts__List] /; Length[{parts}] <= ArrayDepth[m] := (m[[##1]] &) @@  MapThread[ Complement, {Range[Dimensions[m, Length[{parts}]]], {parts}}];
pf[matrix_] := If[Dimensions[matrix][[1]] === 2, matrix[[1, 2]], Sum[(-1)^i matrix[[1, i]] pf[drop[matrix, {1, i}, {1, i}]], {i, 2, Dimensions[matrix][[1]]}]];

We can generate example anti-symmetric matrices like:

size=10;
mat = Table[ If[i > j, i (i - j), j (i - j)], {i, 1, size}, {j, 1, size}];

For size=10;, the function returns the Pfaffian quickly

pf[mat]

-10

However, if I set size=20, the calculation takes forever and never finishes. Considering that the square of the Pfaffian is the determinant Det[mat]==pf[mat]^2, and determinants of larger matrices are still computed quickly by Mathematica:

size=20;
mat = Table[ If[i > j, i (i - j), j (i - j)], {i, 1, size}, {j, 1, size}];
Det[mat]

400

I wonder if it is possible to write an algorithm that computes Pfaffians of larger matrices quickly as well? (I really need the overall sign information for the actual Pfaffians I want to compute, so I cannot fall back to computing determinants and square roots.) Thanks for any suggestion!

$\endgroup$
3
$\begingroup$

Turns out, very efficient implementations have been written by M. Wimmer in the work

https://arxiv.org/abs/1102.3440

In the section "Ancillary Files" to the right, click on "(120 additional files not shown)" and scroll down until you find the line "pfapack/mathematica/pfaffian.nb". This file contains a routine called PfaffianLTL[] which computes the Pfaffian incredibly fast (even for size=200 it takes just a few seconds).

$\endgroup$
  • 1
    $\begingroup$ I was going to suggest something along the lines of using the built-in LUDecomposition, but hadn't worked out the math, yet. :) $\endgroup$ – rcollyer Sep 7 '16 at 19:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.