8
$\begingroup$

The Pfaffian of an even-dimensional anti-symmetric matrix $A$ is defined as:

$$\mathrm{Pf}[A] = \frac{1}{2^{n}n!}\sum_{\pi\in S_{2n}}(-1)^{\pi} a_{i_{1}i_{2}}a_{i_{3}i_{4}}\cdots a_{i_{2n-1}i_{2n}}$$

where the sum is taken over all permutation $\pi$ is:

$$ \begin{equation} \pi = \begin{pmatrix} 1 & 2 & \cdots & 2n \\ i_{1} & i_{2} & \cdots & i_{2n} \end{pmatrix} \end{equation}$$

and $(-1)^{\pi}$ is the parity of the permutation $\pi \in S_{2n}$.

In general, summing over the permutations is a quite daunting task, as an example $S_{6}$ contains 720 terms.

I could not find any built-in function in Mathematica that can evaluate a Pfaffian, which amounts to summing over all permutations for a given $n$. Can anyone shed a light on how to code a Pfaffian in Mathematica? Thanks!

Improve this question
$\endgroup$
2

1 Answer 1

Reset to default
13
$\begingroup$

In the answer to: "https://mathematica.stackexchange.com/questions/125794/compute-numeric-pfaffians-of-matrices-efficiently" there is a pointer to the code below. I hope this will be helpful.

Note also that there is a resource function for Pfaffians at: https://resources.wolframcloud.com/FunctionRepository/resources/Pfaffian/.

The following contains functions to calculate the Pfaffian for real and complex skew-symmetric matrices. Note that the functions don't check (yet) if the input matrix is really skew-symmetric

Functions for skew-symmetric Parlett-Reid tridiagonalization and the corresponding Pfaffian routine

Compute the LTL decomposition of a skew-symmetric matrix using the Parlett-Reid algorithm. The function return T, L and P, where T is a tridiagonal matrix, L a unit lower triangulat matrix and P a permutation matrix, such that P A P^T = L T L^T . Note that this function is not needed for the Pfaffian computation, but is only provided fro demontration purposes.

SkewLTL[Mat_] := Module[{A, L, Pv, N, ip},
  A = Mat;
  N = Dimensions[A][[1]];
  L = IdentityMatrix[N];
  Pv = Range[N];
  For[i = 1, i < N - 1, i++,
   (*find out the maximum entry in the column i, starting from row i+1*)
   ip = i + 
     Position[Abs[A[[i + 1 ;;, i]]], Max[Abs[A[[i + 1 ;;, i]]]]][[1, 1]];
   (*if the maximum entry is not at i+1, permute the matrix so that it is*)
   If[i + 1 != ip, 
    (*Interchange rows and columns in A*)
    A[[{i + 1, ip}, ;;]] = A[[{ip, i + 1}, ;;]]; 
    A[[;; , {i + 1, ip}]] = A[[;; , {ip, i + 1}]];
    (*interchange rows in L; this amounts to accumulating the product of 
    individual Gauss eliminations and permutations*)
    L[[{i + 1, ip}, 1 ;; i]] = L[[{ip, i + 1}, 1 ;; i]];
    (*Accumulate the total permutation matrix*)
    Pv[[{i + 1, ip}]] = Pv[[{ip, i + 1}]];
    ];
   (*Build the Gauss vector*)
   L[[i + 2 ;;, i + 1]] = A[[i + 2 ;;, i]]/A[[i + 1, i]];
   (*Row and column i are eliminated*)
   A[[i + 2 ;;, i]] = 0; A[[i, i + 2 ;;]] = 0;
   (*Update the remainder of the matrix using an outer product skew-
   symmetric update. Note that
   column and row i+1 are not affected by the update*)
   A[[i + 2 ;; , i + 2 ;; ]] += 
    Transpose[{L[[i + 2 ;;, i + 1]]}] . {A[[i + 2 ;;, i + 1]]} - 
     Transpose[{A[[i + 2 ;;, i + 1]]}] . { L[[i + 2 ;;, i + 1]]};
   ];
  Return[{A, L, SparseArray[{i_, i_} -> 1, {N, N}][[Pv]]}]
  ]

Compute the Pfaffian of a skew-symmetric matrix using the LTL decomposition

PfaffianLTL[Mat_] := Module[{A, N, ip, pfaff},
  A = Mat;
  N = Dimensions[A][[1]];
  If[OddQ[N], Return[0]];
  pfaff = 1;
  For[i = 1, i < N - 1, i += 2,
   (*find out the maximum entry in the column i, starting from row i+1*)
   ip = i + 
     Position[Abs[A[[i + 1 ;;, i]]], Max[Abs[A[[i + 1 ;;, i]]]]][[1, 1]];
   (*if the maximum entry is not at i+1, permute the matrix so that it is*)
   If[i + 1 != ip, 
    (*Interchange rows and columns in A*)
    A[[{i + 1, ip}, ;;]] = A[[{ip, i + 1}, ;;]]; 
    A[[;; , {i + 1, ip}]] = A[[;; , {ip, i + 1}]];
    (*interchange contributes det(P)=-1*)
    pfaff = -pfaff;
    ];
   (*Multiply with every other entry on the diagonal*)
   pfaff = pfaff*A[[i, i + 1]];
   (*Build the Gauss vector*)
   A[[i + 2 ;;, i]] = A[[i + 2 ;;, i]]/A[[i + 1, i]];
   (*Update the remainder of the matrix using an outer product skew-
   symmetric update. Note that
   column and row i+1 are not affected by the update*)
   A[[i + 2 ;; , i + 2 ;; ]] += (# - Transpose[#]) &@
    Outer[Times, A[[i + 2 ;;, i]], A[[i + 2 ;;, i + 1]]]
   (* The above is much faster than this construct for me: 
   Transpose[{A[[i+2;;,i]]}].{A[[i+2;;,i+1]]}-Transpose[{A[[i+2;;,i+1]]}].{ 
   A[[i+2;;,i]]};*)
   ];
  Return[pfaff*A[[N - 1, N]]]
  ]

Functions for the Householder tridiagonalization and the corresponding Pfaffian routine

The following contains functions to calculate the Pfaffian for real and complex skew-symmetric matrices. Note that the functions don't check (yet) if the inpiut matrix is really skew-symmetric

Function to compute a so-called Householder vector v for x, i.e. a vector such that (1-2/(v^H v) v v^H) * x is a multiple of the unit vector e_1.

HouseholderVectorReal[x_] := Module[{temp, tempfac, normx},
  normx = Norm[x];
  If[normx == 0, 
   Return[{UnitVector[Dimensions[x][[1]], 1], 0, 0}],
   If[x[[1]] > 0,
    tempfac = normx,
    tempfac = -normx];
   temp = x;
   temp[[1]] += tempfac;
   Return[{Normalize[temp], 2, -tempfac} ]]]

HouseholderVectorComplex[x_] := Module[{temp, tempfac, normx},
  normx = Norm[x];
  If[normx == 0, 
   Return[{UnitVector[Dimensions[x][[1]], 1], 0, 0}],
   tempfac = E^(I Arg[x[[1]]]) Norm[x]; 
   temp = x;
   temp[[1]] += tempfac;
   Return[{Normalize[temp], 2, -tempfac} ]]]

Now the functions that do the full tridagonalization that is at the heart of the Pfaffian calculation. For an input matrix A, they return {T,Q} such that A=Q T Q^T . In the real case, this should be the same as what is returned from the HessenbergDecomposition, for the complex case there is no Mathematica equivalent. Note that these functions are not needed for the Pfaffian calculation, they are here merely for testing.

SkewTridiagonalize[Mat_] := 
 If[MatrixQ[Mat, NumberQ[#] && ! MatchQ[#, _Complex] &], 
  SkewTridiagonalizeReal[Mat], SkewTridiagonalizeComplex[Mat]]

SkewTridiagonalizeReal[Mat_] := Module[{A, Q, v, beta, alpha},
  A = Mat;
  Q = IdentityMatrix[Dimensions[A][[1]]];
  For[i = 1, i < Dimensions[A][[1]] - 1, i++, 
   (*Compute the Householder vector*)
   {v, beta, alpha} = HouseholderVectorReal[A[[i + 1 ;; , i]]];
   (*eliminate the entries in row and column i*) 
   A[[i + 1, i]] = alpha; A[[i, i + 1]] = -alpha; A[[i + 2 ;;, i]] = 0; 
   A[[i, i + 2 ;;]] = 0;
   (*update the matrix*)
   w = beta* A[[i + 1 ;; , i + 1 ;;]] . v; 
   A[[i + 1 ;; , i + 1 ;; ]] += Transpose[{v}] . {w} - Transpose[{w}] . { v};
   (*accumulate the Householder reflections into the full transformation*)
   y = Q[[1 ;;, i + 1 ;;]] . v; 
   Q[[1 ;; , i + 1 ;;]] -= beta*Transpose[{y}] . {v}]; 
  Return[{A, Q}]]

SkewTridiagonalizeComplex[Mat_] := Module[{A, Q, v, beta, alpha},
  A = Mat;
  Q = IdentityMatrix[Dimensions[A][[1]]];
  For[i = 1, i < Dimensions[A][[1]] - 1, i++, 
   (*Compute the Householder vector*)
   {v, beta, alpha} = HouseholderVectorComplex[A[[i + 1 ;; , i]]]; 
   (*eliminate the entries in row and column i*) 
   A[[i + 1, i]] = alpha; A[[i, i + 1]] = -alpha; A[[i + 2 ;;, i]] = 0; 
   A[[i, i + 2 ;;]] = 0;
   (*update the matrix*)
   w = beta* A[[i + 1 ;; , i + 1 ;;]] . Conjugate[v]; 
   A[[i + 1 ;; , i + 1 ;; ]] += Transpose[{v}] . {w} - Transpose[{w}] . { v};
   (*accumulate the Householder reflections into the full transformation*)
   y = Q[[1 ;;, i + 1 ;;]] . v; 
   Q[[1 ;; , i + 1 ;;]] -= beta*Transpose[{y}] . Conjugate[{v}]]; 
  Return[{A, Q}]]

Functions to compute the Pfaffian of a real or complex skew-symmetric matrix.

PfaffianHReal[Mat_] := Module[{N, A, v, beta, alpha, pfaff},
  A = Mat;
  pfaff = 1;
  N = Dimensions[A][[1]];
  If[OddQ[N], Return[0]];
  For[i = 1, i < N - 1, i += 2,
   (*Compute the Householder vector*) 
   {v, beta, alpha} = HouseholderVectorReal[A[[i + 1 ;; , i]]];
   (*multiply with off-
   diagonal entry and determinant of Householder reflection*) 
   pfaff *= If[beta == 0, 1, -1]*(-alpha);
   (*update the matrix*)
   w = beta* A[[i + 1 ;; , i + 1 ;;]] . v; 
   A[[i + 1 ;; , i + 1 ;; ]] += Transpose[{v}] . {w} - Transpose[{w}] . { v}]; 
  Return[pfaff*A[[N - 1, N]]]]

PfaffianHComplex[Mat_] := Module[{N, A, v, beta, alpha, pfaff},
  A = Mat;
  pfaff = 1;
  N = Dimensions[A][[1]];
  If[OddQ[N], Return[0]];
  For[i = 1, i < N - 1, i += 2, 
   (*Compute the Householder vector*) 
   {v, beta, alpha} = HouseholderVectorComplex[A[[i + 1 ;; , i]]]; 
   (*multiply with off-
   diagonal entry and determinant of Householder reflection*) 
   pfaff *= If[beta == 0, 1, -1]*(-alpha);
   (*update the matrix*)
   w = beta* A[[i + 1 ;; , i + 1 ;;]] . Conjugate[v]; 
   A[[i + 1 ;; , i + 1 ;; ]] += Transpose[{v}] . {w} - Transpose[{w}] . { v}]; 
  Return[pfaff*A[[N - 1, N]]]]

PfaffianH[Mat_] := 
 If[MatrixQ[Mat, NumberQ[#] && ! MatchQ[#, _Complex] &], PfaffianHReal[Mat], 
  PfaffianHComplex[Mat]]

Finally another function for testing: Compute the Pfaffian (in the real case) from the HessenbergDecomposition

PfaffianHessenberg::Mat = 
  "Pfaffian computation with Hessenberg decomposition only works for real \
matrices";

PfaffianHessenberg[Mat_] := 
 If[MatrixQ[Mat, NumberQ[#] && ! MatchQ[#, _Complex] &], 
  Module[{H, Q}, {Q, H} = HessenbergDecomposition[Mat]; 
   Return[Det[Q]*Product[H[[i, i + 1]], {i, 1, Dimensions[H][[1]], 2}]]], 
  Message[PfaffianHessenberg::Mat]]

Pfaffian routine combining all of the approaches above

Options[Pfaffian] = {Method -> "ParlettReid"}

{Method -> "ParlettReid"}

Pfaffian::Method = 
  "Unrecognized option `1` (must be either \"ParlettReid\", \"Householder\" \
or \"Hessenberg\")";

Pfaffian[Mat_, OptionsPattern[]] := Switch[OptionValue[Method],
  "ParlettReid", PfaffianLTL[Mat], "Householder", PfaffianH[Mat], 
  "Hessenberg", PfaffianHessenberg[Mat], _, 
  Message[Pfaffian::Method, OptionValue@Method]]

Tests

Now some tests

First for real matrices

A = RandomReal[{0, 1}, {8, 8}]; A = A - A\[Transpose];

Pfaffian[A]

-0.311919

Chop[Det[A] - Pfaffian[A]^2]

0

Pfaffian[A, Method -> "Householder"]

-0.311919

Pfaffian[A, Method -> "Hessenberg"]

-0.311919

{T, L, P} = SkewLTL[A];

P . A . P\[Transpose] - L . T . L\[Transpose] // MatrixForm // Chop



{T, Q} = SkewTridiagonalize[A];

A - Q . T . Q\[Transpose] // MatrixForm // Chop

Then for complex matrices

B = RandomComplex[{0, 1 + I}, {8, 8}]; B = B - B\[Transpose];

Pfaffian[B]

-0.247022 + 0.895042 I

Chop[Det[B] - Pfaffian[B]^2]

0

Pfaffian[B, Method -> "Householder"]

-0.247022 + 0.895042 I

PfaffianHessenberg[B]

PfaffianHessenberg::Mat: Pfaffian computation with Hessenberg decomposition only works for real matrices

{T, L, P} = SkewLTL[B];

P . B . P\[Transpose] - L . T . L\[Transpose] // MatrixForm // Chop


{T, Q} = SkewTridiagonalize[B];

B - Q . T . Q\[Transpose] // MatrixForm // Chop
Improve this answer
$\endgroup$
2
  • $\begingroup$ I was going to suggest something along the lines of Total[Signature[#]*Apply[Times, A[[Sequence@@#]] & /@ Partition[#,2]] & /@ Permutations[Range[2n]]]/(2^n*n!) but it seems to be far inferior to what's shown here. $\endgroup$
    – Roman
    Commented Sep 2, 2022 at 10:00
  • 1
    $\begingroup$ Damn, this answer is very impressive, my +1 sir! $\endgroup$
    – Hans Olo
    Commented Sep 3, 2022 at 9:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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