Division-Free Algorithms
for the Determinant and the Pfaﬃan:
Algebraic and Combinatorial ApproachesGünter Rote
Institut für Informatik, Freie Universität Berlin, Takustraße 9, D-14195 Berlin
rote@inf.fu-berlin.de1

Introduction

The most common algorithm for computing the determinant of an n×n matrix A
is Gaussian elimination and needs O(n3 ) additions, subtractions, multiplications,and divisions. On the other hand, the explicit deﬁnition of the determinant as
the sum of n! products,
sign π · a1π(1) a2π(2) . . . anπ(n)(1)
det A =πshows that the determinant can be computed without divisions. The summation
is taken over the set of all permutations π of n elements. Avoiding divisions seems
attractive when working over a commutative ring which is not a ﬁeld, for example
when the entries are integers, polynomials, or rational or even more complicatedexpressions. Such determinants arise in combinatorial problems, see [11].
We will describe an O(n4 ) algorithm that works without divisions, and wewill look at this algorithm from a combinatorial and an algebraic viewpoint.
We will also consider the Pfaﬃan of a skew-symmetric matrix, a quantity
closely related to the determinant. The results are in many ways analogous tothose for the determinant, but the algebraic aspect of the algorithms is not
explored yet.
This survey is yet another article which highlights the close connection between linear algebra and cycles, paths, and matchings in graphs [1,17,25]. Muchof this material is based on the papers of Mahajan and Vinay [14,15] and of
Mahajan, Subramanya, and Vinay [13].

2

Algorithms for the Determinant

2.1

Alternate Expressions for the Determinant

A direct evaluation of the expression (1) requires n!(n − 1) multiplications (plusa smaller number of additions and subtractions, which will be ignored in thesequel). One can try to rewrite this expression and factorize it in a clever way
such that common subexpressions need to be evaluated only once. Expanding
H. Alt (Ed.): Computational Discrete Mathematics, LNCS 2122, pp. 119–135, 2001.
c Springer-Verlag Berlin Heidelberg 2001


120

Günter Rote

one row or columns leads to a recursive procedure which takes T (n) = n · (1 +T (n − 1)) ≈ n!(e − 1) multiplications. Splitting the matrix
into two equal partsn(1 + B(n/2) + B(n/2))and applying Laplace expansion needs B(n) = n/2
multiplications. This is between 2n (for n > 3) and 22n , which is still exponential.
It will turn out that one has to ﬁrst extend the expression by additional
redundant terms. Only this will allow to reduce the time to evaluate the expression.
2.2

Cycle Decompositions

We ﬁrst rewrite the product a1π(1) a2π(2) . . . anπ(n) in terms of the cycle structureof the permutation π. Every permutation corresponds to a partition of the set{1, . . . , n} into disjoint cycles. For example, the cycle decomposition
C = (8 7 4)(6)(5 3 1 2) = (1 2 5 3)(4 8 7)(6)(2)

corresponds to the following permutation.
i:12345678π=
π(i) : 2 5 1 8 3 6 4 7
We can get a canonical form of a cycle decomposition C by letting each cycle startat its smallest element, which is called the head of the cycle and is underlined in(2), and by rearranging the cycles in increasing order of their heads. Formally,a cycle (c1 , c2 , . . . , ci ) of length i can be deﬁned as a sequence of i ≥ 1 distinct
elements from the ground set {1, . . . , n}, where c1 is the smallest element andis called the head of C. The weight of C is weight(C) = ac1 c2 ac2 c3 . . . aci c1 . Thisis the product of the weight of its arcs if C is regarded as a cycle in a graph.
A cycle decomposition C is a sequence C1 , . . . , Ck of disjoint cycles withincreasing heads, covering every element 1, . . . , n exactly once. The weight of C
is the product of the weights of its cycles, and the sign of a cycle decomposition
C with k cycles is deﬁned assign C = sign π := (−1)n+k .
One can check that this coincides with the usual deﬁnition of the sign and parityof a permutation.
So we can rewrite (1) as follows:
det A =
sign C · weight(C)(3)
cycle decompositions C

2.3

Clow Sequences

If we relax the requirement that all elements of a cycle are distinct, we arriveat the concept of a closed ordered walk or clow for short, in the terminology of
Mahajan and Vinay [14]. A clow (c1 , c2 , . . . , ci ) of length i is a sequence of i ≥ 1

Division-Free Algorithms for the Determinant and the Pfaﬃan121

numbers where the head c1 of C is the unique smallest element: c1 < c2 , . . . , ci .
The weight of C is again weight(C) = ac1 c2 ac2 c3 . . . aci c1 . (Valiant [24], who
introduced this concept, called this a c1 -loop.)
A clow sequence C is a sequence C1 , . . . , Ck of clows with a strictly increasing
sequence of heads:
C = (c1 , . . . ), (c2 , . . . ), . . . , (ck , . . . ),(4)

with c1 < c2 < · · · < ck . The weight of C is the product of the weights of its
clows, its length is the sum of the lengths of the clows, and the sign of C isjust (−1)n+k . It is clear that a cycle decomposition is a special case of a clow
sequence for which the notions of weight and sign coincide with the one given
above.
We can replace the summation over the set of cycle decompositions in formula (3) by the summation over the larger class of clow sequences.
Theorem 1.
det A =

sign C · weight(C)(5)

clow sequences C of length n

Proof. It must be shown that all “bad” clow sequences which are not cycle
decompositions cancel because the corresponding weights occur equally oftenwith a positive and a negative sign. This is done by pairing up the bad clow
sequences such that clow sequences which are matched have opposite signs. This
proof technique is typical of many combinatorial proof for theorems in linear
algebra [20,18,23,25]. The following lemma is stated in Valiant [24]; its proof isgiven in full detail in [14,15].
Lemma 1. Let C be a clow sequence that contains repeated elements. Then we
can ﬁnd a clow sequence C¯ with the same weight and opposite sign. Moreover, thismapping between C and C¯ is an involution, i. e., if we construct the corresponding¯ we obtain C again.
other clow sequence for C,
Proof. Let C = C1 , . . . , Ck . We successively add the clows Ck , Ck−1 , . . . until
a repetition occurs. Assume that Cm+1 , . . . , Ck is a set of disjoint cycles but
Cm , Cm+1 , . . . , Ck contains a repeated element, where 1 ≤ m ≤ k. Let Cm =(c1 , . . . , ci ), and let cj be the ﬁrst element in this clow which is either (A) equal
to a previous element c in the clow (1 <  < j) or (B) equal to an element in
one of the cycles Cm+1 , . . . , Ck . Precisely one of these two cases will occur.
In case (A), we remove the cycle (c =cj , c +1 , . . . , cj−1 ) from the clow Cmand turn it into a separate cycle, which we insert it into C in the proper position
after cyclically shifting its head to the front. In case (B), we insert the cycle into
the clow Cm after the element cj . One easily checks that the operations in cases(a) and (b) are inverses of each other; they change the number of clows by one,and hence they invert the sign. This concludes the proof of the lemma and of
Theorem 1.
✷

122

2.4

Günter Rote

Dynamic Programming According to Length

In some sense, (1) and (3) are the most economical formulas for the determinantbecause all superﬂuous terms have been canceled and no further reduction is possible. However, we shall now see that the redundant form (5) is more amenable toan eﬃcient algorithm. The idea is to construct clow sequences incrementally, arc
by arc: We use dynamic programming to systematically compute the following
quantities:[l, c, c0 , s] is the sum of all partial clow sequences with length l, ending
in the current vertex c, where the head of the current clow is c0 , and thesign s = ±1 gives the parity of the number of ﬁnished clows so far.
A partial clow sequence is an initial piece of the sequence given in (4): It consistsof a number of completed clows, and one incomplete clow. For completing theclow sequence, we only have to know the head and the current last vertex of theincomplete clow. In the end, we need to know the parity of the total number ofclows; and therefore we have to store (−1)k , where k is the number of completed
clows so far.
A dynamic programming algorithm can best be described by a dynamic programming graph. Our graph has O(n3 ) nodes corresponding to the quantities[l, c, c0 , s], with 1 ≤ l ≤ n, 1 ≤ c0 ≤ c ≤ n, and s = ±1.
A partial clow sequence can be grown either by extending the current clow orby closing the current clow and starting a new one. Correspondingly, the dynamic
programming graph has two kinds of outgoing arcs from the node [l, c, c0 , s]: arcs
to nodes [l + 1, c , c0 , s] of cost acc , for all c > c0 , and arcs of cost acc0 , to allnodes [l + 1, c0 , c0 , −s], for all c0 > c0 . The latter type of arcs will start a newclow with head c0 and no edges yet. We consider all nodes [0, c0 , c0 , 1] as startnodes and add two artiﬁcial end nodes [n, n + 1, n + 1, s] with s = ±1. The sumof all path lengths from start nodes to end nodes, taken with appropriate sign
according to s, is the desired quantity in (5), i. e., the determinant.
The sum of all path weights in this acyclic graph can be computed in time
proportional to the number of arcs of this graph, which is O(n4 ): there are O(n)arcs leaving each of O(n3 ) nodes. For example, one can calculate node values in
the order or increasing l values. The storage requirement in this case is O(n2 )because only the nodes with two adjacent l values need to be stored at any time.
Mahajan and Vinay [14] used this algorithm to show that the computationof the determinant belongs to the complexity class GapL, a complexity classin which summation over path weights in an acyclic graphs can be done, but,very loosely speaking, subtraction is allowed only once at the very end. One
can simplify the algorithm by collapsing corresponding nodes [l, c, c0 , +1] and[l, c, c0 , −1] with opposite sign ﬁelds into one node [l, c, c0 ] and introducing “negative” arc weights: the second type of arcs have then cost −acc0 instead of acc0 .
2.5

Adding a Vertex at a Time: Combinatorial Approach

The previous algorithm considers all partial clow sequences in order of increasing
length l. A diﬀerent factorization is possible: we can ﬁrst consider the clow with

Division-Free Algorithms for the Determinant and the Pfaﬃan123

head c1 = 1, if such a clow is present, and the remaining clow sequence, which isa sequence on the ground set {2, . . . , n}. For expressing the clows with head 1 wepartition the matrix A into its ﬁrst row r, its ﬁrst column s, and the remaining
matrix M .


a11 r
A=
s M
The sum of all clows of length l can be written as follows:
length l = 0: 1
length l = 1: a11nlength l = 2:
a1i ai1 = rs
i=2

length l = 3:..
.
length l:n na1i aij aj1 = rM s

i=2 j=2

rM l−2 s

The term 1 for length l = 0 corresponds to the “empty” clow consisting of noedges. It accounts for the case when the head of the ﬁrst clow is larger than 1.
We denote by qi the sum of signed weights of all clow-sequences of length n−1−iwhich don’t contain the element 1. Here we deﬁne the sign of a clow sequencewith k clows as (−1)n−1+k . Then we can write:
det A = −1 · q−1+ a11 · q0 + rs · q1 + · · · + rM i−1 s · qi + · · ·(6)The minus sign accounts for the change of the sign convention when going from
q−1
over a ground set of n−1 elements to the original ground set with n elements.
In the other terms this minus sign is canceled because we have one additional
clow.


By Lemma 1, the quantity q−1
must be zero: q−1
is the signed sum of allclow sequences of length n on the ground set {2, . . . , n}. Such sequences must
contain repeated elements and hence all terms cancel. So we get
det A = a11 · q0 + rs · q1 + · · · + rM i−1 s · qi + · · · + rM n−2 s · qn−1.

We intend to compute the values qi recursively, so we need to consider clow
sequences of all possible lengths. Let us denote by qi the sum of signed weights
of all clow-sequences of length n − i (with the original deﬁnition of the sign).
So q0 = det A, and q0 = det M . By the above argument, we get the followingformula for computing q0 , q1 , q2 , . . . , qn from the sequence q0 , q1 , . . . , qn−1.

124Günter Rote
0   qn−1
qn−1  
0 


 
 qn−2qn−2  

0 

 
  ..  ..  .
.
..  
 . =   q  q2  rM n−4 s rM n−5 s · · · a11 −1 0   2    q 1
 q1  rM n−3 s rM n−4 s · · · rs a11 −1 
q0
n−2
n−3
s · · · rM s rs a11
rM
q0
s rM
qn





−1
a11
rs
..
.

0
−1
a11
..
.

···
···
···
..
.

0
0
0
..
.

0
0
0
..
.

(7)

This matrix multiplication is nothing but a convolution between the sequenceq0 , q1 , . . . and the sequence −1, a11 , rs, rM s, . . . , which can most conveniently
be expressed as a polynomial multiplication:
qn λn + qn−1 λn−1 + · · · + q1 λ + q0
= (−λ + a11 + rsλ−1 + rM sλ−2 + · · · + rM i sλ−i−1 + · · · )× (qn−1
λn−1 + qn−2
λn−2 + · · · + q1 λ + q0 )(8)

Actually, these are not polynomials, but something like Laurent series in λ−1 .
In explicit form, this means

 

−1
0
··· 0
0 0 0 ··· qnqn−1   a11−1··· 0
0 0 0 · · · 
 qn−1  qn−2   rs
a11
··· 0
0 0 0 · · ·  qn−2  ..  

. 
..
..
.
.
.
.
.
.
..
..
..
.. .. . .   .. 
 .  
.
.


 


 q1  q2  rM n−4 s rM n−5 s · · · −1
0 0 0 · · ·


=


 q1  rM n−3 s rM n−4 s · · · a11 −1 0 0 · · ·  q0 


 


 q0  rM n−2 s rM n−3 s · · · rs a11 −1 0 · · ·  0 


 


 0  rM n−1 s rM n−2 s · · · rM s rs a11 −1 · · ·  0 

 



 0   rM n s rM n−1 s · · · rM 2 s rM s rs a11 · · ·
..

 

.
..
..
..
..
..
.. .. . .
..
.
.
.
.
.
.
.
. .
This is an extended version of (7) which reveals the “hidden parts” of (8). For
example, the line corresponding to q−1 = 0 yields the relation

0 = rs · q0 + rM s · q1 + · · · + rM i s · qi + · · · + rM n−1 s · qn−1,

which may be an interesting identity but is of little use for computing anything.
Mahajan and Vinay [14] have analyzed the recursive algorithm based on (7)more carefully, and they have shown that the quantities qi , after omitting irrel
evant terms like q−1
from (6), actually represent a subset of the clow-sequencesthat are characterized by the so-called preﬁx property.
2.6

The Characteristic Polynomial

We have associated the polynomial
PA (λ) = qn λn + qn−1 λn−1 + · · · + q1 λ + q0

Division-Free Algorithms for the Determinant and the Pfaﬃan

125

to the matrix A, and we have shown how PA (λ) can be computed from thecorresponding polynomial PM (λ) of the (n − 1) × (n − 1) submatrix M :PA (λ) = (−λ + a11 + rsλ−1 + rM sλ−2 + · · · + rM i sλ−i−1 + · · · ) · PM (λ) (9)This allows us to inductively compute polynomials for larger and larger submatrices, starting from P(ann ) (λ) = ann − λ, until we ﬁnally obtain PA (λ) and thedeterminant det A = q0 .
We have seen that the inductive approach of adding a vertex at a time naturally leads us to consider not only clow sequences of length n as in (5), butalso clow sequences of all other lengths. We have also found it convenient to
use the corresponding sums of weights as coeﬃcients of a polynomial PA (λ). It
turns out that this polynomial is nothing but the characteristic polynomial ofthe matrix A:
PA (λ) := det(A − λI) = qn λn + qn−1 λn−1 + · · · + q1 λ + q0

(10)

The following well-known formula generalizes formula (3) for the determinant
to all coeﬃcients of the characteristic polynomial.
qi =sign C · weight(C)(11)
families C of disjoint cycles with length n − iHere we consider not partitions into cycles, but families C of disjoint cycles witha total of n − i elements. The sign of such a family with k cycles is deﬁned as(−1)n−i+k . Formula (11) comes from the fact that a cycle family which covers
n − i vertices leaves i vertices uncovered, at which one can place the diagonal
entries −λ when the determinant of A − λI is expanded.
The following theorem states that the quantities qi deﬁned in the previous
section are indeed the coeﬃcients of the characteristic polynomial (10) as weknow it. It is proved in the same way as Theorem 1, see [14].
Theorem 2. The coeﬃcients qi (i = 0, 1, . . . , n) of the characteristic polynomial of A can be expressed assign C · weight(C),
(12)
qi =
clow sequences C of length n − iwhere the sign of a clow sequence of length n − i with k cycles is deﬁned as✷(−1)n−i+k .
The dynamic programming algorithm of Section 2.4 can easily be extended
to compute all coeﬃcients qi : Simply adding artiﬁcial terminal nodes [n − i, n +1, n + 1], with arcs of the correct sign, will do the job.
2.7

Adding a Vertex at a Time: Algebraic Approach

We will now give an algebraic proof of (7). The determinant of a bordered matrix


a11 r
(13)
A=
s M

126

Günter Rotecan be expressed asdet A = det M · (a11 − rM −1 s)(14)

if the matrix M is invertible. This follows from the fact that det M is the cofactorof a11 , and hence (A−1 )11 = det M/ det A. The element (1, 1) of the inverse A−1
can be determined by block Gauss-Jordan elimination of (13), starting with thelower right block M . This gives (A−1 )11 = (a11 − rM −1 s)−1 which leads to (14).
Now we wish to express PA (λ) = det(A − λI) in terms of PM (λ) = det(M −λI). Applying (14) to the matrix A − λI gives a11 − λrPA (λ) = |A − λI| = sM − λI = |M − λI| · (a11 − λ − r(M − λI)−1 s)= PM (λ) · (ann − λ + r(λI − M )−1 s) =: PM (λ) · F (λ)
For simplifying the expression F (λ), which is the multiplier which turns PM into
PA , we rewrite (λI − M )−1 as follows.
(λI − M )−1 = λ−1 (I − λ−1 M )−1= λ−1 (I + λ−1 M + λ−2 M 2 + · · · )= λ−1 I + λ−2 M + λ−3 M 2 + · · ·This gives
F (λ) = a11 − λ + r(λI − M )−1 s= −λ + a11 + rsλ−1 + rM sλ−2 + rM 2 sλ−3 + rM 3 sλ−4 + · · · ,(15)which is the same factor as in (8) and (9).
Another proof of (7), based on the Cayley-Hamilton Theorem, is due to
P. Beame, see [3].
2.8

Who Invented the Incremental Algorithm?

Berkowitz [3] used the algorithm based on the recursion (7) to derive a parallel
procedure for evaluating the determinant. He called this algorithm Samuelson’s
method, giving a reference to Samuelson’s original paper [19] and to the textbook
of Faddeev and Faddeeva [8]. Valiant [24], who gave the ﬁrst combinatorial interpretation of (7), referred to the algorithm as the Samuelson-Berkowitz algorithm.
However, Samuelson’s algorithm for computing the characteristic polynomial isan O(n3 ) algorithm that uses divisions. This is also the way in which it is described in [8, §45]. So the above algorithm can deﬁnitely not be attributed toSamuelson. However, a passage towards the end of Section §45 in [8] contains thesystem (7). Faddeyev and Faddeyeva state this as a possible stricter justiﬁcation of Samuelson’s method. The mention brieﬂy that this formula follows from

Division-Free Algorithms for the Determinant and the Pfaﬃan127

the consideration of bordered determinants. This calculation has been sketchedabove in Section 2.7.
The fact that (7) appears in the section about Samuelson’s method in thebook [8] may be the reason for the erroneous attribution. It seems that themethod ought to be credited to Faddeyev and Faddeyeva. However, being numerical analysts, they would probably not have considered the implied O(n4 )procedure as a good “method”. Berkowitz [3] was apparently the ﬁrst one to use
the recursive formula (7) algorithmically.

3

The Pfaﬃan

The determinant of a skew-symmetric matrix (aij = −aji ) is the square ofanother expression, which is called the Pfaﬃan [12,16]. For example,
 0 a12 a13 a14  −a12 0 a23 a24 2 −a13 −a23 0 a34  = (a12 a34 − a13 a24 + a14 a23 ) −a14 −a24 −a34 0 
The determinant of a skew-symmetric matrix of odd order is always 0. Formally,
the Pfaﬃan of a skew-symmetric matrix A with an even number of rows andcolumns can be deﬁned as follows:
pf A =
sign M · weight(M)(16)
perfect matchings M

Here, a perfect matching M with k = n/2 edges is written as
M = (i1 , j1 ), (i2 , j2 ), . . . , (ik , jk ),(17)

where, by convention, il < jl for each l. The sign of the matching M is deﬁned
as the sign of
1 2 3 4 ... n − 1 n
i1 j1 i2 j2 . . . ik jkwhen this is regarded as a permutation of {1, . . . , n}. The weight of M isai1 j1 ai2 j2 . . . aik jk . One can check that the sign of M does not depend on theorder in which the edges are given. Moreover, exchanging a pair il , jl changes
the sign but not the signed weight since aji = −aij . This means that we are infact free to choose any convenient permutation π as a representation of a given
matching M.
Because of their close connection with matchings, Pfaﬃans are of interestin combinatorics, see [12,21]. For example, for certain graphs, including planargraphs, it is possible to count the number of perfect matchings in polynomial
time by using Pfaﬃans. Knuth [10] gives a short history of Pfaﬃans, and argues
that Pfaﬃans are in some way more fundamental than determinants, to which
they are closely related (see Section 4).

128Günter Rote

The Pfaﬃan can be computed in O(n3 ) steps by an elimination procedurethat is similar to Gaussian elimination for symmetric matrices. Alternatively, thesquare root of the determinant gives the Pfaﬃan up to the correct sign. Bothapproaches may be undesirable if divisions or square roots should be avoided.
3.1

Alternating Cycle Decompositions

The combinatorial approaches to the determinant of Section 2 can be applied to
Pfaﬃans as well. We consider the product of the Pfaﬃans of two skew-symmetricmatrices A and B. Expanding the deﬁnitions (16) leads to an overlay of two
matchings MA and MB . The union of MA with MB decomposes into a disjoint
union of alternating cycles with a total of n edges:
MA  MB = C = C1 ∪ C2 ∪ · · ·An alternating cycle is a cycle of even length whose edges alternate between
MA and MB . Edges which are contained both in MA and MB must be kept
separate the union MA  MB ; for each “common” edge, MA  MB contains
two (distinguishable) parallel edges, forming an alternating cycle of length 2.
An alternating cycle can be traversed in two orientations. We specify theorientation by taking the vertex with the smallest number as the “ﬁrst vertex”
or head of the cycle and insisting that the ﬁrst edge should belong to MA .
Taking the union of two matchings brings us back into the realm of cycle
decompositions which we studied in Section 2. Now, the weight of an alternating
cycle C = (c1 , c2 , . . . , ci ) with c1 < c2 , . . . , ci isweight(C) = ac1 c2 bc2 c3 ac3 c4 bc4 c5 . . . aci−1 ci bci c1 .

(18)

As in Section 2.2, an alternating cycle family C is a set of disjoint alternating
cycles. Its length is the total number of elements in its cycles. If the length is n,
we speak of an alternating cycle decomposition. The weight of C is the product of
the weights of its cycles, and the the sign of an alternating cycle decomposition
C with k cycles is deﬁned as sign C = (−1)k .
Theorem 3.
pf A · pf B =sign C · weight(C)alternating cycle decompositions C

Proof. The Pfaﬃan has been deﬁned as a weighted sum of matchings; so theleft side is a sum over all pairs of matchings (MA , MB ), where MA is regarded as a matching for A and MB is regarded as a matching for B. Now,
there is a unique correspondence between pairs of matchings and alternating
cycle decompositions. The only thing that has to be checked is that the signs
are correct. Suppose that MA  MB decomposes into k disjoint cycles C =(c1 , c2 , . . . , ci ), (d1 , d2 , . . . , dj ), . . . . We may conveniently represent MA and MB
by the permutations
1 2 3 4 ··· i − 1 i i + 1 i + 2 ··· i + j i + j + 1 ··· nπA =c1 c2 c3 c4 · · · ci−1 ci
d1 d2 · · · dj···Division-Free Algorithms for the Determinant and the Pfaﬃan

and
πB =



1 2 3 4 ··· i − 1 i i + 1 i + 2 ··· i + j
c2 c3 c4 c5 · · · ci c1
d2 d3 · · · d1

i + j + 1 ··· n
···

129


,

respectively. This gives sign MA · sign MB = sign πA · sign πB = sign((πA )−1 ◦πB ) = (−1)k , since the permutation (πA )−1 ◦ πB decomposes into k even-length
cycles.
✷As an easy consequence, we obtain the relation between the determinant and
the Pfaﬃan that was mentioned at the beginning of this section.
Corollary 1.

det A = (pf A)2

Proof. When we set A = B in Theorem 3, we almost get the determinant formula (3): the diﬀerence between alternating cycle decompositions and cycle decompositions disappears except for the requirement that an alternating cycle
must alway have even length.
Changing the orientation of an odd cycle in the expansion of the determinantof a skew-symmetric matrix reverses the weight to its negative. Thus it is easy
to establish a sign-reversing bijection among all cycle decompositions containing
odd cycles, for example by reversing the odd cycle with the largest head. So thecycle decompositions containing odd cycles cancel.
✷
If we want to use Theorem 3 to get an expression for the Pfaﬃan of a single
matrix A, we must set B to a matrix B0 whose Pfaﬃan is 1. So we select the
skew-symmetric matrix


0 1
−1 0





0
1




−1
0
B0 = 
(19)



.
..





0 1
−1 0
that represents the matching M0 = (1, 2), (3, 4), . . . , (n − 1, n), which we use asa reference matching. So we immediately get the following alternate expression
for the Pfaﬃan.
Corollary 2.
pf A =

sign C · weight(C)✷M0 -alternating cycle decompositions C

Here, the general concepts of alternating cycles, cycles families, etc., specialize
as follows: An M0 -alternating cycle (c1 , c2 , . . . , ci ) is a cycle of even length i. In
addition, we must have c3 = c2 − 1 if c2 is even and c3 = c2 + 1 if c2 is odd, andthe same relation holds between c4 and c5 , . . . , ci and c1 . Its weight is given by(18) with respect to B0 , i. e., b2j−1,2j = 1 and b2j,2j−1 = −1. All other notions
are unchanged.

130

3.2

Günter Rote

Alternating Clow Sequences

Generalizing Theorem 1 in Section 2.3, Mahajan, Subramanya, and Vinay [13]have extended the notion of an alternating cycle cover to an alternating clow
sequence. This gives the following theorem.
Theorem 4.
pf A · pf B =sign C · weight(C)(20)

sign C · weight(C)(21)

alternating clow sequences C of length n

pf A =

M0 -alternating clow sequences C of length nHere, an alternating clow (c1 , c2 , . . . , ci ) of length i is deﬁned like a clow, withthe additional restriction that i must be even. The condition on the clow headis unchanged: c1 < c2 , . . . , ci . The weight of alternating clows is deﬁned likein (18) for alternating cycles. An M0 -alternating clow (c1 , c2 , . . . , ci ) must haveevery other edge (c2 , c3 ), (c4 , c5 ), . . . , (ci , c1 ) belonging to M0 . Alternatingclow sequences, their length, weight, and sign are deﬁned just as for clows andalternating cycle families. (Our notation and terminology diﬀers from [13].)
The proof of (20), and of (21) which is a special case of it, follows the linesof the proof of Lemma 1: a sign-reversing bijection among “bad” alternatingclow sequences which contain repeated elements is established. The only change
occurs in case (A), when a clow Ck forms a sub-cycle (c =cj , c +1 , . . . , cj−1 ). If
that sub-cycle has odd length, we cannot split it from the current clow because it
does not form an alternating cycle. In this case, we simply reverse the orientation
and replace this part of Ck by (c =cj , cj−1 , . . . , c +1 ). The weight of Ck changes
sign, and it can also be checked easily that the resulting mapping between badalternating clow sequences remains an involution.
✷
Mahajan, Subramanya, and Vinay [13] considered only the case (21) of asingle matrix, and they gave a slightly diﬀerent proof.
3.3

Incremental Algorithms for the Pfaﬃan

A dynamic programming algorithm for evaluating (21) by summing over allpartial alternating clow sequences in order of increasing length can easily beformulated in analogy to Section 2.4, see [13]. However, we shall instead formulate another algorithm along the lines of the recursive algorithm based on (7) in
Section 2.5.
Since Pfaﬃans only make sense for matrices of even order, we partition thegiven skew-symmetric matrix by splitting two rows and columns from it.

0 a12 rT
A =  −a12 0 −s 
−rT s M

Division-Free Algorithms for the Determinant and the Pfaﬃan

131

Enumeration of alternating clow sequences with head c1 = 1 according to length
yields the power series
G(λ) = −λ2 + a12 + rB0 sλ−2 + rB0 M B0 sλ−4 + rB0 M B0 M B0 sλ−6 + · · ·(22)
B0 is the “skew-symmetric unit matrix” (19) of appropriate dimension, here(n − 2) × (n − 2). The coeﬃcient of λ−i+2 , for i ≥ 2, represents the negative
sum of all weights of all alternating clow sequences with head 1 and length i.
Note that a clow starting with vertex 1 must terminate with the edge (2, 1), andhence the factor b21 = −1 is alway present in the above terms. We deﬁne thePfaﬃan-characteristic polynomial of a skew-symmetric n × n matrix A asP̃A (λ) := q̃n λn + q̃n−2 λn−2 + · · · + q̃2 λ2 + q̃0 ,where the quantities q̃i denote the sum of signed weights of all M0 -alternatingclow-sequences of length n − i. This deﬁnition can also be used for odd n. Similarly as in Theorem 2, q̃i is equal to the sum of signed weights of all M0 alternating cycle decompositions of length n − i. It can also be interpreted thesum of signed weights of certain matchings with (n − i)/2 edges as in (16), with
an appropriate deﬁnition of signs.
Note that an alternating clow cannot have head 2. So we split alternating clow
sequences into the (possibly empty) clow with head 1 and the remaining clow
sequence, which consists of elements {3, . . . , n} only, in analogy to Section 2.5.
We obtain the recursive equationP̃A (λ) = G(λ) · P̃M (λ),(23)

which immediately leads to a recursive division-free algorithm for computingP̃A (λ) and the Pfaﬃan pf A = q̃0 .

4Open Questions

The relation


det A = pf

0 A
−AT 0



shows that determinants are just special cases of Pfaﬃans that correspond to
bipartite graphs. In this sense, the Pfaﬃan is a more basic notion, which would
deserve a more thorough understanding, despite the traditional prevalence of thedeterminant in the curriculum and in applications. We mention a few questions
that arise from the results which we have surveyed.
4.1

The Pfaﬃan-Characteristic Polynomial:
Algebraic Interpretation

As in the case of the determinant, our algorithm has naturally lead to a polynomial P̃A (λ) of degree n that is associated to a skew-symmetric matrix A of132

Günter Roteorder n. In contrast to the case of the determinant, we have no idea what analgebraic interpretation of this polynomial might be, or how one might prove
the relation (23) by algebraic means. It isnot generally the case that the characteristic polynomial is the square of P̃A ( (λ)), as the relation q0 = q̃02 between
their constant terms would suggest.
In Theorems 3 and 4, it has turned out that it appears more natural todeal with the product of Pfaﬃans of two diﬀerent matrices A and B: after all,Corollary 2 is only a special case of Theorem 3 which is by no means easier toformulate, and the same relation holds between (20) and (21) in Theorem 4.
Many identities about Pfaﬃans involve products of two Pfaﬃans that are takenfrom submatrices of a single given matrix. One might speculate that taking two
diﬀerent matrices and their submatrices could lead to more general theorems
and at the same time more transparent proofs is some other cases as well.
Instead of substituting the matrix B0 from (19) into (20) of Theorem 4, one
could also use the matrix
0 1 1 1 ··· 1
−1 0 1 1 · · · 1
−1 −1 0 1 · · · 1
B1 = −1 −1 −1 0 · · · 1 , .. .. .. .. . . .. 
 . . . .
. .−1 −1 −1 −1 · · · 0

which might be called the skew-symmetric analog of the all-ones matrix andhas pf B1 = 1. This results in a diﬀerent expression for pf A than (21) in Theorem 4 and to another algorithm. That algorithm is in a way simpler than theincremental algorithm given in Section 3.3, since it increments the matrix by
one row and column at a time instead of two. However, it uses more computations. The matrix B1 also leads to a diﬀerent concept of a Pfaﬃan-characteristic
polynomial.
There is an analogous algorithm for directly computing the product of two
Pfaﬃans according to equation (20) of Theorem 4. It seems that the Pfaﬃancharacteristic polynomial ought to be associated to two matrices instead of one.
The incremental algorithm corresponding to equation (20) leads to a polynomial
P̃A,B (λ) that is related to the characteristic polynomial by the equation
P̃A,A (λ) = PA (λ),
for skew-symmetric matrices A.
Is there a way to introduce a second matrix B also into the computation ofdeterminants in the same natural way? This might be related to the generalized
eigenvalue problem that requires the solution of the characteristic equation
det(A − λB) = 0,
for an arbitrary matrix B instead of the usual identity matrix. The above characteristic equation is equal to the characteristic equation of AB −1 if B is invertible.
This problem is important for some applications. Can any of the algorithms of
Section 2 be adapted to compute the coeﬃcients of this polynomial?

Division-Free Algorithms for the Determinant and the Pfaﬃan

4.2

133

Complexity

All division-free methods for computing the determinant or the Pfaﬃan that
we have described take O(n4 ) steps. This has already been mentioned for the
dynamic programming algorithm of Section 2.4. In the recursive algorithms (7)and (23) of Sections 2.5 and 3.3, the step of the recursion from the (n−1)×(n−1)
matrix M to the n × n matrix A takes O(n3 ) time for computing the coeﬃcientsof the factor F (λ) or G(λ) respectively, and O(n2 ) steps for multiplying the twopolynomials, if this is done in the straightforward way. This amounts to a totalof O(n4 ) steps.
The bottleneck is thus the computation of the iterated matrix-vector products in (15) and (22): for F (λ) in (15), we successively compute the row vectors
r, rM , rM 2 , . . . , and multiply these by the column vector s. When the matrix
A is sparse, this computation will be faster: If A has m nonzero elements, eachmatrix-vector product will take only O(m) steps instead of O(n2 ) (assumingm ≥ n). This leads to an algorithm with a running time of O(mn2 ). Thus, whenm is small, this algorithm may be preferable to ordinary Gaussian eliminationeven when divisions are not an issue. Due to its simplicity, the iterative algorithm gains an additional advantage over Gaussian elimination or other direct
methods, which must be careful about the choice of pivots in order to avoid
numerical problems or to exploit sparsity.
Sparsity will also improve the dynamic programming algorithm of Section 2.4
from O(n4 ) to O(mn2 ) steps. It is a challenging problem to devise a division-freeO(n3 ) method for computing the determinant.
The algorithms discussed in this paper can also be used to design good parallel algorithms for determinants and Pfaﬃans [3,14,13].
Despite their advantages, the methods described in this article are not wellknown. A straightforward implementation of the incremental algorithm of Section 2.5 in Maple started to beat the built-in Maple procedure for the characteristic polynomial for random integer matrices of order 20, and was clearlyfaster already for matrices of order 40. For a simple 5 × 5 matrix with rational
1
entries like (aij ) = ( xi +x), our algorithm for the characteristic polynomial wasj
two times faster than Maple’s algorithm for computing the determinant only.
4.3

Alternative Algorithms

In the introduction it was mentioned that division-free algorithms might be useful for evaluating determinants of integers. The most straightforward algorithm
would be to carry out Gaussian elimination with rational arithmetic. This approach may result in an explosive growth in the lengths of the numbers unlessfractions are reduced after every step. However, to avoid the excessive growthof numbers, it is not necessary to reduce fractions by their greatest common
divisor: one can perform pivots involving divisions for which the result is known
to be integral. This form of integral pivoting is attributed to Camille Jordan [5,p. 69], see also [6,2,9,7]. Edmonds [6] has shown that the numbers do not explode134

Günter Roteand hence, this algorithm runs in polynomial time. The same approach can beextended to Pfaﬃans [9].
Thus, for evaluating determinants or Pfaﬃans of integers, this approach issuperior to our division-free algorithms because it takes only O(n3 ) operations.
However, when dealing with multivariate polynomials, even “integer division”
can be awkward, and division-free algorithms may prove to be worthwhile.
Strassen [22] has given a general recipe for converting an algorithm that uses
divisions to obtain an integral ﬁnal result into an algorithm without divisions. Divisions are simulated by expanding expressions into truncated power series. Thisis somehow reminiscent of the recursion (9) which also involves power series thatcan be truncated because the result is a polynomial. According to [4], Strassen’s
approach leads to an algorithm with O(n5 ) steps in the case of Gaussian elimination. It would be interesting to see if there is any relation of this algorithm tothe algorithms presented here.

References
1. Martin Aigner, Lattice paths and determinants. In: Computational Discrete Mathematics, ed. Helmut Alt, (this volume), Lecture Notes Comput. Sci., Vol. 2122 2001,
pp. 1–12.
2. Erwin H. Bareiss, Sylvester’s identity and multistep integer-preserving Gaussianelimination. Math. Comput. 22 (1968), 565–578.
3. Stuart J. Berkowitz, On computing the determinant in small parallel time using asmall number of processors. Inf. Process. Lett. 18 (1984), 147–150.
4. Allan Borodin, Joachim von zur Gathen, John Hopcroft, Fast parallel matrix and
GCD computations. Inf. Control 52 (1982), 241–256
5. E. Durant, Solution numérique des équations algébriques, tome II: systèmes desplusieurs équations. Masson & Cie., Paris 1961.
6. Jack Edmonds, Systems of distinct representatives and linear algebra. J. Res. Nat.
Bur. Standards, Sect. B 71 (1967), 241–245.
7. Jack Edmonds, J.-F. Maurras, Note sur les Q-matrices d’Edmonds. RAIRO, Rech.
opér. 31 (1997), 203–209.
8. D. K. Faddeyev, V. N. Faddeyeva, Vyčislitel nye metody lineı̆noı̆ algebry (in Russian), Moscow, 1960. English translation: D. K. Faddeev, V. N. Faddeeva, Numerical Methods of Linear Algebra. Freeman, San Francisco 1963. German translation:
D. K. Faddejew, W. N. Faddejewa, Numerische Methoden der linearen Algebra,
several editions since 1964.
9. G. Galbiati and Franceso Maﬃoli, On the computation of pfaﬃans. Discr. Appl.
Math. 51 (1994), 269–275.
10. Donald E. Knuth, Overlapping Pfaﬃans. Electron. J. Comb. 3 (1996), No. 2, article R5, 13 pp. Printed version: J. Comb. 3 (1996), No. 2, 147–159.
11. Christian Krattenthaler, Advanced determinant calculus. Séminaire Lotharingien
de Combinatoire B42q (1999), 67 pp.
12. László Lovász, M. D. Plummer, Matching Theory. Ann. Discr. Math., Vol. 29.
North-Holland Mathematics Studies, Vol. 121. Amsterdam 1986.
13. Meena Bhaskar Mahajan, P R Subramanya, V Vinay, A combinatorial algorithm
for Pfaﬃans. In: Computing and combinatorics. Proc. 5th annual international
conference. (COCOON ’99), Tokyo, July 1999, ed. Takao Asano et al., LectureDivision-Free Algorithms for the Determinant and the Pfaﬃan

135

Notes Comput. Sci. 1627, Springer-Verlag, pp. 134–143 (1999). Extended version:
DIMACS Technical Report 99-39, Rutgers University, July 1999.
14. Meena Bhaskar Mahajan, V Vinay, Determinant: Combinatorics, algorithms, andcomplexity. Chicago J. Theor. Comput. Sci., Vol. 1997, Article no. 1997-5, 26 pp.
15. Meena Mahajan, V. Vinay, Determinant: Old algorithms, new insights. SIAM J.
Discrete Math. 12 (1999), 474–490.
16. Thomas Muir, A Treatise on the Theory of Determinants. MacMillan and Co.,
London 1882; repr. Dover, New York 1960.
17. Günter Rote, Path problems in graphs. In: Computational graph theory, ed. Gottfried Tinhofer et al., Computing Suppl. 7, 155–189, Springer-Verlag, Wien 1990.
18. D. E. Rutherford, The Cayley-Hamilton theorem for semi-rings. Proc. Roy. Soc.
Edinburgh, Sect. A 66 (1961–64), 211–215 (1964).
19. Paul A. Samuelson, A method of determining explicitly the coeﬃcients of thecharacteristic equation. Ann. Math. Statist. 13 (1942), 424–429.
20. Dennis Stanton, Dennis White, Constructive combinatorics. Springer-Verlag, NewYork 1986.
21. John R. Stembridge, Nonintersecting paths, pfaﬃans, and plane partitions. Adv.
Math. 83 (1990), 96–113.
22. Volker Strassen, Vermeidung von Divisionen. J. reine angew. Math. 264 (1973),
184–202.
23. Howard Straubing, A combinatorial proof of the Cayley-Hamilton theorem. Discrete Math. 43 (1983), 273–279.
24. Leslie G. Valiant, Why is Boolean complexity theory diﬃcult? In: Boolean Function
Complexity, ed. M. S. Paterson, LMS Lecture Notes Series, Vol. 169, Cambridge
Univ. Press, 1992, pp. 84–94.
25. Doron Zeilberger, A combinatorial approach to matrix algebra. Discrete Math. 56(1985), 61–72.