The Faddeev-LeVerrier algorithm and the Pfaffian
Introduction
The computation of the determinant of a large -matrix A is a challenging task because naively using the formula would require many multiplications since the are n! permutations σ of . This is practical only for very small matrices.
The Gauss algorithm allows to transform a matrix in triangular form with 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 and it avoids any division by a matrix entry. Here is the computational cost of the multiplication of two -matrices. The precise value of β is unknown but it does lie in the interval , 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 -algebras R. We could as well allow any torsion-free commutative ring R because such an R embeds into its rationalization via . So the algorithm can be used if , , , , 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 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 . 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 -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.
Section snippets
Faddeev-LeVerrier algorithm for the characteristic polynomial
Let R be a commutative -algebra. In particular, this covers the case , , and . In Section 6 we will consider a more sophisticated algebra. Let where denotes the set of all -matrices with entries in R. We want to compute the coefficients of its characteristic polynomial , in particular its determinant.
Definition of the Pfaffian
Denote the space of skew-symmetric -matrices by . Let be even and let . The Pfaffian of is defined as where the sum is taken over all perfect matchings P of . Here a perfect matching is a partition of the set into a disjoint union of m sets with two elements each, . If we use the convention that each pair is ordered by then the sign of P is defined as the sign of the
The Faddeev-LeVerrier algorithm for the Pfaffian
Now we are ready to adapt the Faddeev-LeVerrier algorithm to compute the Pfaffian.
Performance
There are many perfect matchings of . Thus a direct implementation of the Pfaffian based on (7) would require many multiplications which is unpractical unless the matrix is very small.
Obviously, the multiplication of two -matrices can be performed at computational cost . Surprisingly, Strassen in [22] found a way to perform matrix multiplication at cost with . The optimal exponent β is unknown to date. It clearly satisfies and after
An application in differential geometry
In this last section we want to illustrate the usefulness of the algorithm for symbolic computation. Our version of the Faddeev-LeVerrier algorithm is division free in the sense that we never have to divide by an entry of the skew-symmetric matrix A whose Pfaffian we are computing. This not only avoids numerical instabilities but also allows us to consider matrices with entries in general commutative -algebras in which division may not be possible. We will make use of this now.
Let M be an
Declaration of Competing Interest
None declared.
References (24)
On computing the determinant in small parallel time using a small number of processors
Inf. Process. Lett.
(1984)- et al.
Numeric and symbolic evaluation of the Pfaffian of general skew-symmetric matrices
Comput. Phys. Commun.
(2011) On the modified Leverrier-Faddeev algorithm
Linear Algebra Appl.
(1982)A modified Leverrier-Faddeev algorithm for matrices with multiple eigenvalues
Linear Algebra Appl.
(1980)- et al.
On Faddeev-Leverrier's methods for the computation of the characteristic polynomial of a matrix and of eigenvectors
Linear Algebra Appl.
(1993) - et al.
The combinatorial approach yields an NC algorithm for computing Pfaffians
Discrete Appl. Math.
(2004) - et al.
An improved parallel processor bound in fast matrix inversion
Inf. Process. Lett.
(1978) - et al.
A new extension of Leverrier's algorithm
Linear Algebra Appl.
(1993) - et al.
A refined laser method and faster matrix multiplication
Geometric Algebra
(1957)