A simple division-free algorithm for computing Pfaffians

TL;DR

This paper introduces a straightforward, division-free algorithm for computing Pfaffians that leverages matrix multiplication and truncation, offering efficiency especially for sparse matrices and enabling extraction of related polynomials.

Contribution

The paper presents a novel, simple algorithm for Pfaffian computation that avoids division and adapts the Bird determinant algorithm for this purpose.

Findings

Algorithm has complexity O(nM(n)) for 2n x 2n matrices.

Efficient extraction of characteristic and Pfaffian characteristic polynomials.

Suitable for sparse matrices with optimized dense-sparse multiplication.

Abstract

We present a very simple algorithm for computing Pfaffians which uses no division operations. Essentially, it amounts to iterating matrix multiplication and truncation. Its complexity, for a $2n\times 2n$ matrix, is $O(nM(n))$, where $M(n)$ is the cost of matrix multiplication. In case of a sparse matrix, $M(n)$ is the cost of the dense-sparse matrix multiplication. The algorithm is an adaptation of the Bird algorithm for determinants. We show how to extract, with practically no additional work, the characteristic polynomial and the Pfaffian characteristic polynomial from these algorithms.

Equations

$$\mu(X)=\left[\begin{array}[]{cccc}-\sum_{i=2}^{n}x_{ii}&x_{12}&\ldots&x_{1n}\\ 0&-\sum_{i=3}^{n}x_{ii}&\ldots&x_{2n}\\ \ldots&\ldots&\ldots&\ldots\\ 0&0&\ldots&-\sum_{i=n+1}^{n}x_{ii}\end{array}\right].$$

$$\mathop{\rm pf}\nolimits A[\alpha]=\sum\big{\{}(-1)^{j}a_{1j}\mathop{\rm pf}\nolimits A[\alpha\setminus 1j]\mid j\in[2..2n]\big{\}}$$

$$x^{p,p-1}_{ij}=\sum\big{\{}\mathop{\rm pf}\nolimits A[i\alpha j]\cdot\mathop{\rm pf}\nolimits B[\alpha]\mid\alpha\in S_{2p-2}(\beta_{i})\big{\}},$$

$$x^{p,p}_{ij}=\sum\big{\{}\mathop{\rm pf}\nolimits A[i\alpha]\cdot\mathop{\rm pf}\nolimits B[\alpha j]\mid\alpha\in S_{2p-1}(\beta_{i})\big{\}},$$

$$x_{11}^{n,n}=\mathrm{pf}A[1{\beta }_1]\cdot \mathrm{pf}B[{\beta }_{1}1]=-\mathrm{pf}A[[\mathrm{1..2}n]]\cdot \mathrm{pf}B[[\mathrm{1..2}n]]=-\mathrm{pf}A\cdot \mathrm{pf}B.$$

$$x_{ij}^{1,0}=\mathrm{pf}A[ij]\cdot \mathrm{pf}B[\epsilon ]=a_{ij}$$

$$x^{1,1}_{ij}=\sum\big{\{}\mathop{\rm pf}\nolimits A[ik]\cdot\mathop{\rm pf}\nolimits B[kj]\mid k\in S_{1}(\beta_{i})\big{\}}=\sum_{k=i+1}^{2n}a_{ik}b_{kj}$$

$$x_{ij}^{p,p}$$

$$\displaystyle=\sum\big{\{}\mathop{\rm pf}\nolimits A[i\alpha k]\cdot\mathop{\rm pf}\nolimits B[\alpha]b_{kj}\mid\alpha\in S_{2p-2}(\beta_{i}),k\in\beta_{i}\big{\}}.$$

$$\mathop{\rm pf}\nolimits B[\alpha j]=-\mathop{\rm pf}\nolimits B[j\alpha]=-\sum\big{\{}b_{jk}\mathop{\rm pf}\nolimits B[\alpha\setminus k]\cdot\mathop{\rm sgn}\nolimits\left({\alpha\atop k(\alpha\setminus k)}\right)\biggm{|}k\in\alpha\big{\}}$$

$$-\sum\big{\{}\mathop{\rm pf}\nolimits A[ik(\alpha\setminus k)]\cdot\mathop{\rm pf}\nolimits B[\alpha\setminus k]b_{jk}\mid\alpha\in S_{2p-1}(\beta_{i}),k\in\alpha\big{\}}.$$

$$\sum$$

$$=$$

$$\{$$

$$=\{(ik({\alpha }^{\prime }\setminus k),{\alpha }^{\prime }\setminus k)\mid {\alpha }^{\prime }\in S_{2p-1}(\beta ),k\in {\alpha }^{\prime }\}.$$

$$d_{ii}^{p,p}=-\sum \{x_{kk}^{p,p}\mid k\in {\beta }_i\}.$$

$$d^{p,p}_{ii}=-\sum\big{\{}\mathop{\rm pf}\nolimits A[k\alpha]\cdot\mathop{\rm pf}\nolimits B[\alpha k]\mid k\in\beta_{i},\alpha\in S_{2p-1}(\beta_{k})\big{\}}.$$

$$d^{p,p}_{ii}=\sum\big{\{}\mathop{\rm pf}\nolimits A[\alpha]\cdot\mathop{\rm pf}\nolimits B[\alpha]\mid\alpha\in S_{2p}(\beta_{i})\big{\}}.$$

$$x_{ij}^{p+1,p}=d_{ii}^{p,p}{a}_{ij}+\sum \{x_{ik}^{p,p}{a}_{kj}\mid k\in {\beta }_i\}.$$

$$x_{ij}^{p+1,p}=$$

$$\displaystyle+\sum\big{\{}\mathop{\rm pf}\nolimits A[i\alpha]\cdot\mathop{\rm pf}\nolimits B[\alpha k]a_{kj}\mid\alpha\in S_{2p-1}(\beta_{i}),k\in\beta_{i}\big{\}}.$$

$$\sum\big{\{}\mathop{\rm pf}\nolimits A[i\alpha j]\cdot\mathop{\rm pf}\nolimits B[\alpha]\mid\alpha\in S_{2p}(\beta_{i})\big{\}}.$$

$$\mathrm{pf}A[i\alpha j]$$

$$\displaystyle=a_{ij}\mathop{\rm pf}\nolimits A[\alpha]-\sum\big{\{}a_{jk}\mathop{\rm pf}\nolimits A[i(\alpha\setminus k)]\cdot\mathop{\rm sgn}\nolimits\left({i\alpha\atop ki(\alpha\setminus k)}\right)\biggm{|}k\in\alpha\big{\}}.$$

$$\sum$$

$$\displaystyle+\sum\big{\{}a_{jk}\mathop{\rm pf}\nolimits A[i(\alpha\setminus k)]\cdot\mathop{\rm pf}\nolimits B[k(\alpha\setminus k)]\mid\alpha\in S_{2p}(\beta_{i}),k\in\alpha\big{\}}.$$

$$\sum$$

$$\displaystyle=\sum\big{\{}\mathop{\rm pf}\nolimits A[i(\alpha\setminus k)]\cdot\mathop{\rm pf}\nolimits B[k(\alpha\setminus k)]a_{kj}\mid\alpha\in S_{2p}(\beta),k\in\alpha\big{\}}$$

$$\{$$

$$=\{(i({\alpha }^{\prime }\setminus k),k({\alpha }^{\prime }\setminus k))\mid {\alpha }^{\prime }\in S_{2p}(\beta ),k\in {\alpha }^{\prime }\}.$$

$$\det (\lambda I-A)={\lambda }^n+c_1{\lambda }^{n-1}+c_2{\lambda }^{n-2}+\cdots +c_{n-1}\lambda +c_n$$

$$x_{ii}^{p}=-(-1)^{p}\sum\big{\{}\det A[i\alpha]\mid\alpha\in S_{p-1}([i+1..n])\big{\}};$$

$$c_{p}=\sum\big{\{}\det(-A[\alpha])\mid\alpha\in S_{p}([1..n])\big{\}}=-\sum\big{\{}x_{ii}^{p}\mid i\in[1..n]\big{\}}=-\mathop{\rm tr}F^{p}_{A}.$$

$$\mathrm{pf}(\lambda B-A)$$

Full text

A simple division-free algorithm for computing Pfaffians

Adam J. Przeździecki

Adam Przeździecki, Warsaw University of Life Sciences—SGGW, Warsaw, Poland

adam_przezdziecki@sggw.edu.pl

Abstract. We present a very simple algorithm for computing Pfaffians which uses no division operations. Essentially, it amounts to iterating matrix multiplication and truncation. Its complexity, for a $2n\times 2n$ matrix, is $O(nM(n))$, where $M(n)$ is the cost of matrix multiplication. In case of a sparse matrix, $M(n)$ is the cost of the dense-sparse matrix multiplication.

The algorithm is an adaptation of the Bird algorithm for determinants. We show how to extract, with practically no additional work, the characteristic polynomial and the Pfaffian characteristic polynomial from these algorithms.

Mathematics Subject Classification. 15A15.

Keywords. Combinatorial problems, Pfaffian, Division-free algorithms, Characteristic polynomial.

1. Introduction

Let $X=[x_{ij}]$ be an $n\times n$ matrix. Bird [2] defined

[TABLE]

Given $A=[a_{ij}]$, define $F^{1}_{A}=A$, and inductively, $F^{p+1}_{A}=\mu(F^{p}_{A})A$. Compared to Bird’s original notation, we increase the upper index $p$ by $1$ so that $p$ becomes equal to the degree of the entries of $F^{p}_{A}$ viewed as homogeneous polynomials in the $a_{ij}$. He proves

Theorem 1.1** (Bird).**

The matrix $F^{n}_{A}$ is everywhere zero except for its leading entry, which equals $(-1)^{n-1}\det A$.

Let $A$, $B$ be two $2n\times 2n$ skew-symmetric matrices. Define $F^{1,0}_{AB}=A$, and inductively, $F^{p,p}_{AB}=\mu(F^{p,p-1}_{AB})B$ and $F^{p+1,p}_{AB}=\mu(F^{p,p}_{AB})A$. We prove

Theorem 1.2**.**

The matrix $F^{n,n}_{AB}$ is everywhere zero except for its leading entry, which equals $-\mathop{\rm pf}\nolimits A\cdot\mathop{\rm pf}\nolimits B$.

Thus if we are interested in the Pfaffian of a single matrix $A$, it is enough to set $B$ to any matrix whose Pfaffian is $1$ and change the sign. In Section 3 we discuss three examples of such matrices and describe how to extract the characteristic polynomial from the Bird algorithm and the Pfaffian characteristic polynomial from our modification.

Clearly, the complexity of computing the determinant or the Pfaffian by means of Theorems $1$ or $2$ is $O(nM(n))$, where $M(n)$ is the complexity of matrix multiplication. We agree with Bird’s claim that “it is difficult to imagine a simpler procedure for computing the determinant”. Recently, Bär [1] published a similarly simple algorithm for the Pfaffian and the Pfaffian characteristic polynomial, but his algorithm is not entirely division-free as it requires division by integers up to half of the dimension of the matrix. Bär remarks that his algorithm can be accelerated to $O(\sqrt{n}M(n))$ by using ideas from Preparata and Sarwate [7], but this is done at a significant cost to the simplicity of the algorithm.

2. Proof of Theorem 1.2

The result can be proved by an adaptation of Rote’s [8] or Bird’s [2] argument. We choose the latter. For the most part we follow Bird’s notation. Given a word $\alpha$ in the alphabet $[1..2n]$ and a matrix $\displaystyle A=[a_{ij}]_{i,j=1}^{2n}$, let $A[\alpha]$ denote the matrix $\displaystyle[a_{\alpha_{i}\alpha_{j}}]_{i,j=1}^{k}$ where $k$ is the length of $\alpha$. Note that if $A$ is skew-symmetric then so is $A[\alpha]$ and for any $i$, $j$ we have $\mathop{\rm pf}\nolimits A[ij]=a_{ij}$. An odd permutation of letters in $\alpha$ changes the sign of $\mathop{\rm pf}\nolimits A[\alpha]$. The proof presented below depends only on the Laplace-type expansion for Pfaffians, (see Fulton and Pragacz [3, Equation D.1, p. 116]):

[TABLE]

where $\alpha\setminus\beta$ is the word $\alpha$ with the symbols in $\beta$ removed.

We derive explicit formulas for the matrices $F^{p,p-1}_{AB}$ and $F^{p,p}_{AB}$ in terms of the Pfaffians of submatrices of $A$ and $B$. Let $S_{p}(\alpha)$ denote the set of subsequences of $\alpha$ of length $p$. We claim that, for skew-symmetric $2n\times 2n$ matrices $A$ and $B$, the $ij$-th entry $x^{p,p-1}_{ij}$ of $F^{p,p-1}_{AB}$ and $x^{p,p}_{ij}$ of $F^{p,p}_{AB}$ are

[TABLE]

[TABLE]

where $\beta_{i}=[i+1..2n]$ and $i\alpha j$, $i\alpha$, $\alpha j$ denote the appropriate concatenations of words. Assuming (2.3), we see that $S_{2n-1}(\beta_{i})$ is nonempty only for $i=1$, in which case it has only one element, $\beta_{1}$, and therefore

[TABLE]

The change of sign is caused by the odd permutation of letters in $\beta_{1}1\mapsto[1..2n]$. As in [2] we see that $x^{n,n}_{1j}=0$ for $j>1$ since $\beta_{1}j$ contains $j$ twice. This shows that (2.3) implies Theorem 1.2.

We prove (2.2) and (2.3) by induction. For $p=1$ we agree that the Pfaffian of a $0\times 0$ matrix is $1$ and we verify (2.2):

[TABLE]

where $\varepsilon$ is the empty word. Verifying (2.3), we see that $S_{1}(\beta_{i})=\{i+1,...,2n\}$ and

[TABLE]

is the $ij$-th entry of the matrix product of the above diagonal part of $A$ by $B$, as required.

We use (2.2) and (2.3) as the induction hypothesis and prove their analogues for $p$ increased by $1$.

Proof that (2.2) implies (2.3). The diagonal of $\mu(F^{p,p-1}_{AB})$ is null since $\mathop{\rm pf}\nolimits A[i\alpha j]$ in (2.2) is null for $i=j$.

For $i<j$ we use (2.2) to expand the $ij$-th entry of $F^{p,p}_{AB}=\mu(F^{p,p-1}_{AB})B$:

[TABLE]

In order to show that (2.4) equals the right hand side of (2.3) we use (2.1) to expand the second Pfaffian in (2.3) along its first row/column and obtain

[TABLE]

where $\mathop{\rm sgn}\nolimits\left({\alpha\atop k(\alpha\setminus k)}\right)$ denotes the sign of the permutation which brings $k$ to the front. Thus (2.3) equals

[TABLE]

To prove that (2.4) equals (2.5) it is sufficient to show

[TABLE]

for every word $\beta$ without repeated letters. The change of sign is due to $b_{kj}=-b_{jk}$.

For $\alpha$ of even length, we have $\mathop{\rm pf}\nolimits A[i\alpha k]=\mathop{\rm pf}\nolimits A[ik\alpha]$ and if $k\in\alpha$ then $\mathop{\rm pf}\nolimits A[i\alpha k]=0$; therefore, (2.6) reduces to

[TABLE]

Take $\alpha^{\prime}\in S_{2p-1}(\beta)$ and $k\in\alpha^{\prime}$. Then $\alpha^{\prime}\setminus k\in S_{2p-2}(\beta)$ and $k\in\beta\setminus\alpha$. Conversely, take $\alpha\in S_{2p-2}(\beta)$ and $k\in\beta\setminus\alpha$. Then $k$ can be inserted in a unique position in $\alpha$ to give a word $\alpha^{\prime}\in S_{2p-1}(\beta)$. Furthermore, $\alpha=\alpha^{\prime}\setminus k$.

Proof that (2.3) for $p$ implies (2.2) for $p+1$. The $i$-th diagonal element of $\mu(F^{p,p}_{AB})$ is

[TABLE]

Applying (2.3) as induction hypothesis, we obtain

[TABLE]

Since $S_{2p}(\beta_{i})=\{k\alpha\mid k\in\beta_{i},\alpha\in S_{2p-1}(\beta_{k})\}$ and $\mathop{\rm pf}\nolimits B[\alpha k]=-\mathop{\rm pf}\nolimits B[k\alpha]$ for $\alpha$ of odd length, we have

[TABLE]

Using (2.8) we compute the $ij$-th entry of $F^{p+1,p}_{AB}=\mu(F^{p,p}_{AB})A$:

[TABLE]

We use (2.8) and the induction hypothesis (2.3) to expand the above as

[TABLE]

We intend to prove that (2.9) equals the right hand side of (2.2) for $p$ increased by $1$, which is

[TABLE]

For $\alpha$ of even length we have $\mathop{\rm pf}\nolimits A[i\alpha j]=-\mathop{\rm pf}\nolimits A[ji\alpha]$ and we use (2.1) to expand the first Pfaffian in (2.2) along its first row/column to obtain

[TABLE]

Therefore (2.10) is equal to

[TABLE]

To prove that (2.9) equals (2.11) it is sufficient to show

[TABLE]

for every word $\beta$ without repeated letters. Note the changes of signs $\mathop{\rm pf}\nolimits B[\alpha k]=-\mathop{\rm pf}\nolimits B[k\alpha]$ on the left, and $a_{kj}=-a_{jk}$ on the right.

If $k\in\alpha$ then $\mathop{\rm pf}\nolimits B[k\alpha]=0$ and therefore (2.12) reduces to

[TABLE]

This holds by the same argument as in the proof of (2.7).

3. The characteristic and Pfaffian-characteristic polynomials

Both the Bird algorithm and our adaptation to Pfaffians allow extraction of the characteristic polynomials at negligible additional cost.

Proposition 3.1**.**

The coefficients of the characteristic polynomial

[TABLE]

are given by $c_{p}=-\mathop{\rm tr}F^{p}_{A}$.

Proof.

By [2, Equation (1)] the $i$-th diagonal entry of $F^{p}_{A}$ is

[TABLE]

note that the index $p$ is increased by $1$ relative to Bird’s original notation. We have

[TABLE]

∎

In the case of the Pfaffian polynomial

[TABLE]

the choice of the skew-symmetric matrix $B=[b_{ij}]$ which should play the role of the identity matrix is not obvious. Below we cite three examples of such matrices $B$ which are most often met in the literature or for other reasons interesting.

[TABLE]

The Pfaffian of each of these matrices is $1$. The matrix $B_{0}$ is sparse and it is the easiest one to guess; it is considered for example by Bär [1] and Rote [8]. From the combinatorial point of view, which is closest to the author, the matrix $B_{1}$ is more interesting since for every $\alpha\in S_{2p}([1..2n])$ we have $\mathop{\rm pf}\nolimits B_{1}[\alpha]=1$. For this reason $B_{1}$ is considered by Rote [8] and by Iwata [4] who also uses $B_{2}=-B_{1}^{-1}$. Knus, Merkurjev, Rost and Tignol [5] and Krivoruchenko [6] work with $\mathop{\rm pf}\nolimits B\cdot\mathop{\rm pf}\nolimits(\lambda B^{-1}-A)$ for general, invertible $B$.

We have the following lemma for the Pfaffian of the sum of $B_{2}$ and another skew-symmetric matrix.

Lemma 3.2** (Stembridge [9, Lemma 4.2]).**

If $B_{2}$ is as above, then

[TABLE]

The closest analogy for Proposition 3.1 is

Proposition 3.3**.**

The coefficients of the Pfaffian characteristic polynomial

[TABLE]

are given by $c_{p}=(-1)^{p+1}\mathop{\rm tr}F^{p,p}_{AB_{1}}$.

Note that the algorithm uses $B_{1}$ while in the definition of the polynomial we have $B_{2}=-B_{1}^{-1}$.

Proof.

Lemma 3.2 implies

[TABLE]

Equation (2.3) for the pair of matrices $A$, $B_{1}$ and $\mathop{\rm pf}\nolimits B_{1}[\alpha]=1$ yields

[TABLE]

hence

[TABLE]

∎