The Faddeev-LeVerrier algorithm and the Pfaffian

Abstract

We adapt the Faddeev-LeVerrier algorithm for the computation of characteristic polynomials to the computation of the Pfaffian of a skew-symmetric matrix. This yields a very simple, easy to implement and parallelize algorithm of computational cost $\mathrm{O}(n^{\beta +1})$ where n is the size of the matrix and $\mathrm{O}(n^{\beta })$ is the cost of multiplying $n\times n$-matrices, $\beta \in [2,2.37286)$. We compare its performance to that of other algorithms and show how it can be used to compute the Euler form of a Riemannian manifold using computer algebra.

Introduction

The computation of the determinant of a large $n\times n$-matrix A is a challenging task because naively using the formula $\det (A)={\sum }_{\sigma }\mathrm{sign}(\sigma )a_{1\sigma (i)}\cdots a_{n\sigma (n)}$ would require $\mathrm{O}((n-1)n!)$ many multiplications since the are n! permutations σ of $\{1,\dots ,n\}$. This is practical only for very small matrices.

The Gauss algorithm allows to transform a matrix in triangular form with $\mathrm{O}(n^3)$ many operations. The determinant can then be read off. However, the Gauss algorithm includes divisions by matrix entries, which requires a pivoting strategy to avoid numerical instabilities. More seriously, it cannot be carried out if the matrix takes entries in a commutative ring in which we cannot divide.

It is therefore desirable to find a fast division-free algorithm. The classical Faddeev-LeVerrier algorithm provides precisely this; its computational complexity is of the order $\mathrm{O}(n^{\beta +1})$ and it avoids any division by a matrix entry. Here $\mathrm{O}(n^{\beta })$ is the computational cost of the multiplication of two $n\times n$-matrices. The precise value of β is unknown but it does lie in the interval $[2,2.37286)$, compare the discussion at the beginning of Section 5.

Moreover, the Faddeev-Leverrier algorithm is very simple and easy to implement. One essentially has to carry out n matrix multiplications; hence using any software which parallelizes matrix multiplication will automatically parallelize the whole algorithm. In Section 2 we will recall the algorithm and present the probably shortest possible derivation. It even simplifies the approach presented in [12]. However, no claim of originality is made for this part. See [3], [9], [10], [11], [23] for various extensions of the algorithm.

The Faddeev-LeVerrier algorithm is division free in the sense that no divisions by matrix entries are required. But we do have to divide by integers. For this reason we allow matrices with entries in general commutative $\mathbb{Q}$-algebras R. We could as well allow any torsion-free commutative ring R because such an R embeds into its rationalization $\mathbb{Q}{\otimes }_{\mathbb{Z}}R$ via $r\mapsto 1\otimes r$. So the algorithm can be used if $R=\mathbb{Z}$, $R=\mathbb{Q}$, $R=\mathbb{R}$, $R=\mathbb{C}$, or if R is any field of characteristic 0. If R has torsion, e.g. if R is a finite field, then different algorithms are required, see e.g. [4].

Our interest lies in the fast computation of the Pfaffian. The Pfaffian is defined for skew-symmetric matrices and is a polynomial in its entries which squares to the determinant. The Pfaffian may be less prominent than the determinant but it has important applications in physics, in combinatorics, and in geometry. As to physics, see the introduction of [8], [24] and the references therein and the application to the topological charge of a disordered nanowire in [24]. In combinatorics, the Pfaffian of the directed adjacency matrix of a suitable graph yields the number of perfect matchings. In Section 6 we will see an application to differential geometry.

As for determinants, the computation of a Pfaffian using the definition to be given in (7) would be way too slow. There are established algorithms of order $\mathrm{O}(n^3)$ which transform the matrix in a normal form from which the Pfaffian can be read off. This is similar to computing the determinant using Gauss elimination. These algorithms are well suited for the numerical treatment of real or complex matrices and, again, cannot be applied to matrices with entries in commutative rings.

We adapt the Faddeev-LeVerrier algorithm to compute the Pfaffian in Section 4 and obtain an algorithm of order $\mathrm{O}(n^{\beta +1})$. Again, it is very simple, easy to implement and can easily be parallelized. The derivation of this algorithm requires some theoretical material about Pfaffians which we provide in Section 3 since it seems hard to find a good presentation in the literature. The algorithm is again division free in the sense that no divisions by matrix entries are required.

In Section 5 we compare the performance of our algorithm to that of other established algorithms. Based on an earlier version of this paper, our algorithm has been implemented in SageMath 9.3 which allows for an easy comparison with the previous implementation using the definition of the Pfaffian.

The main application we have in mind is to compute the Pfaffian of matrices taking entries in algebras of symbolic expressions as typically used in computer algebra. In Section 6 we show how the algorithm can be used to compute the Euler form of a Riemannian manifold using SageMath. In this case, the entries of the matrix are formal differential forms of mixed even degree. Thus we are dealing with a commutative $\mathbb{Q}$-algebra in which we cannot divide. In particular, algorithms based on transforming the matrix in a canonical form do not apply. Moreover, multiplication in this algebra is very costly as it involves expensive simplification routines. Our algorithm appears to provide the best known approach to deal with this kind of scenario.

Acknowledgment. The author wants to thank Zachary Hamaker for pointing out the complexity results for matrix multiplication and Darij Grinberg and an unknown referee for many comments and suggestions on how to improve the presentation. Moreover, he is grateful for financial support by SPP 2026 funded by Deutsche Forschungsgemeinschaft.