# Compute numeric Pfaffians of matrices efficiently?

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!

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).
• I was going to suggest something along the lines of using the built-in LUDecomposition, but hadn't worked out the math, yet. :) – rcollyer Sep 7 '16 at 19:23