Pfaffian skew-symmetric using armadillo

Question

I have been trying to compute the pfaffian of skew-symmetric matrices. Currently, I am using the armadillo library. But the time complexity is exponential. Can you please suggest some efficient ways?

Here is the function using the armadillo library;

double pfaffian(arma::mat Test, int N){

  if(N == 0){
    return 1.0f;
  }
  else  if (N == 2){
    return Test(0, N2-1);
  } 
  else{
    double pfa = 0.0f;
    for(int j = 2; j <= N; j++){
        arma::mat temp = Test;

        cout << "j: " << j << endl;
        temp.shed_row(0);
        temp.shed_col(0);

        temp.shed_row(j-2);
        temp.shed_col(j-2);

        pfa += pow((-1),j) * Test(0,j-1)
                * pfaffian(temp, N-2);
    }
    return pfa;
   }
}

Answer 1

The Pfaffian has the property $\newcommand{\pf}{\operatorname{pf}}$ $$ \pf(OAO^T)=\det(O)\pf(A) $$ If $O$ is a Householder reflector, then $\det(O)=-1$. The reduction to Hessenberg form requires reflectors to reduce the contents of column $i$ from row $i+1$ down to have zeros from row $i+2$ downwards, and this from column $1$ to $2N-2$. Thus in total the signs of the reflectors cancel.

As these orthogonal similarity transformation preserve the skew-symmetry, the resulting Hessenberg matrix is tri-diagonal. The Pfaffian of that is just the product over the odd-index entries of the super-diagonal. $$ \pf\pmatrix{ 0&a_1&0&\ldots\\ -a_1&0&a_2&0&\ldots\\ 0&-a_2&0&a_3&0&\ldots\\ &&\ddots&&\ddots\\ &&0&-a_{2N-2}&0&a_{2N-1}\\ &&&0&-a_{2N-1}&0 }=a_1a_3\cdots a_{2N-3}a_{2N-1} $$ The even-indexed entries can be removed using row and their transposed column operations where the transformation matrix has determinant $1$, without changing the odd-index entries.

(I'm not sure why the wikipedia page does not mention this, as the properties leading to this are clearly exhibited.)

The Hessenberg transform has complexity $O(N^3)$, clearly less than the $O((2N)!)$ of the recursive formula.