Hafnian of a matrix

Question

Is there a Mathematica built-in function for the calculation of the Hafnian of a matrix? I am interested in the function as it appears naturally in the calculation of moments of a multivariate Gaussian distribution.

Answer 1

Mathematica

Use the recursive formula.

Try it online!

haf[matrix_?MatrixQ] := 
  haf[matrix] = Module[{n, h, j, copy}, n = Length[matrix];
    If[n == 2, Return[matrix[[1, 2]]]];
    h = 0;
    For[j = 2, j <= n, j++, If[matrix[[1, j]] == 0, Continue[]];
     copy = Delete[matrix, {{1}, {j}}];
     copy = Transpose[copy];
     copy = Delete[copy, {{1}, {j}}];
     h += matrix[[1, j]]*haf[copy];];
    h];

Answer 2

Using the match function by Mr.Wizard.

match[x : {_, _}] := {{x}}
match[{a_, b__}] := 
 Join @@ Table[{{a, x}, ##} & @@@ match[{b}~DeleteCases~x], {x, {b}}]
match[n_?EvenQ] := match@Range[n]
match[n_?OddQ] := 
 DeleteCases[match@Range[n + 1], {___, n + 1, ___}, {2}]

Code for Hafnian:

haf[m_?MatrixQ] := 
 With[{d = Dimensions[m]}, 
  If[And[Divide @@ d == 1, And @@ EvenQ[d]], 
   Module[{mtc = match[Length[m]]},
    Total[Times @@@ (Extract[m, #] & /@ mtc)]
    ], {}]]

Testing:

haf[Array[a, {6, 6}]]

a[1, 6] a[2, 5] a[3, 4] + a[1, 5] a[2, 6] a[3, 4] + a[1, 6] a[2, 4] a[3, 5] + a[1, 4] a[2, 6] a[3, 5] + a[1, 5] a[2, 4] a[3, 6] + a[1, 4] a[2, 5] a[3, 6] + a[1, 6] a[2, 3] a[4, 5] + a[1, 3] a[2, 6] a[4, 5] + a[1, 2] a[3, 6] a[4, 5] + a[1, 5] a[2, 3] a[4, 6] + a[1, 3] a[2, 5] a[4, 6] + a[1, 2] a[3, 5] a[4, 6] + a[1, 4] a[2, 3] a[5, 6] + a[1, 3] a[2, 4] a[5, 6] + a[1, 2] a[3, 4] a[5, 6]

Using test cases from @138 Aspen

mat1 = {{-1, 1, 1, -1, 0, 0, 1, -1}, {1, 0, 1, 0, -1, 0, -1, -1}, {1, 
    1, -1, 1, -1, -1, 0, -1}, {-1, 0, 1, -1, -1, 1, -1, 
    0}, {0, -1, -1, -1, -1, 0, 0, -1}, {0, 0, -1, 1, 0, 0, 1, 
    1}, {1, -1, 0, -1, 0, 1, 1, 0}, {-1, -1, -1, 0, -1, 1, 0, 1}};
mat2={{-1, 1, 0, 1, 0, -1, 0, 0, -1, 1, -1, 1, 0, -1}, {1, -1, 1, -1, 1, 
   1, -1, 0, -1, 1, 1, 0, 0, -1}, {0, 1, 1, 1, -1, 1, -1, -1, 0, 
   0, -1, 0, -1, -1}, {1, -1, 1, -1, 1, 0, 1, 1, -1, -1, 0, 0, 1, 
   1}, {0, 1, -1, 1, 0, 1, 0, 1, -1, -1, 1, 1, 0, -1}, {-1, 1, 1, 0, 
   1, 1, -1, 0, 1, -1, -1, -1, 1, -1}, {0, -1, -1, 1, 0, -1, -1, -1, 
   0, 1, -1, 0, 1, -1}, {0, 0, -1, 1, 1, 0, -1, 0, 0, -1, 0, 0, 0, 
   1}, {-1, -1, 0, -1, -1, 1, 0, 0, 1, 1, 0, 1, -1, 0}, {1, 1, 
   0, -1, -1, -1, 1, -1, 1, 1, 1, 0, 1, 0}, {-1, 1, -1, 0, 1, -1, -1, 
   0, 0, 1, -1, 0, -1, 0}, {1, 0, 0, 0, 1, -1, 0, 0, 1, 0, 0, 1, 1, 
   1}, {0, 0, -1, 1, 0, 1, 1, 0, -1, 1, -1, 1, 1, -1}, {-1, -1, -1, 
   1, -1, -1, -1, 1, 0, 0, 0, 1, -1, -1}};
haf[mat1]
haf[mat2]

yields 4 and 83 respectively.

I look forward to learning from other approaches.

As @JM alludes it soons becomes intractable (number perfect matches $\frac{(2n)!}{n! 2^n}$ , as can be be seen by partitioning 2n into pairs and accounting for n pairs and 2 ways to order each pair).