A59:01

A59:02

A59:03

A59:04 A59:05

A59:06

A59:07 A59:08 A59:09 A59:10

A59:11

A59:12 A59:13 A59:14 A59:15 A59:16 A59:17 A59:18 A59:19 A59:20 A59:21 A59:22 A59:23 A59:24 A59:25 A59:26

A59:27

A59:28 A59:29 A59:30 A59:31 A59:32 A59:33

Determinants and Pfaffians: Algebraic and Combinatorial Approaches

Günter Rote

Institut für Informatik, Freie Universität Berlin, Takustraße 9, D-14195 Berlin,

                   Germany, email: rote@inf.fu-berlin.de

1 Introduction

The most common algorithm for computing the determinant of an n×n matrix A is Gaussian elimination and needs O(n^{3}) additions, subtractions, multiplications, and divisions. On the other hand, the explicit definition of the determinant as the sum of n! products,

det A = ∑_{π} sign π · a_{1π(1)}a_{2π(2)} … a_{nπ(n)}  (1)

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 field, for example when the entries are integers, polynomials, or rational or even more complicated expressions. Such determinants arise in combinatorial problems, see [11].

We will describe an O(n^{4}) algorithm that works without divisions, and we will look at this algorithm from a combinatorial and an algebraic viewpoint.

We will also consider the Pfaffian of skew-symmetric matrix, a quantity closely related to the determinant. The results are in many ways analogous to those 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]. Much of 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 (plus a smaller number of additions and subtractions, which will be ignored in the sequel). 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 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

A60:01 A60:02 A60:03 A60:04 A60:05

A60:06

A60:07 A60:08 A60:09

A60:10

A60:11

A60:12

A60:13 A60:14 A60:15 A60:16 A60:17 A60:18 A60:19 A60:20 A60:21 A60:22 A60:23

A60:24

A60:25 A60:26 A60:27

A60:28

A60:29

A60:30 A60:31 A60:32 A60:33

and applying Laplace expansion needs B(n) = \binom{n}{⌊n/2⌋}(1+B(⌊n/2⌋)+B(⌈n/2⌉)) multiplications. This is between 2^{n} (for n > 3) and 2^{2n}, which is still exponential.

It will turn out that one has to first extend the expression by additional redundant terms. Only this will allow to reduce the time to evaluate the expression.

2.2 Cycle Decompositions

We first rewrite the product a_{1π(1)}a_{2π(2)} … a_{nπ(n)} in terms of the cycle structure of 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) = (\underline{1} 2 5 3)(\underline{4} 8 7)(\underline{6})  (2)

corresponds to the following permutation.

π = (\begin{matrix}i : & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 \\ π(i) : & 2 & 5 & 1 & 8 & 3 & 6 & 4 & 7 \end{matrix})

We can get a canonical form of a cycle decomposition C by letting each cycle start at 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 (c_{1}, c_{2}, …, c_{i}) of length i can be defined as a sequence of i ≥ 1 distinct elements from the ground set {1, …, n}, where c_{1} is the smallest element and is called the head of C. The weight of C is weight(C) = a_{c_{1}c_{2}}a_{c_{2}c_{3}} … a_{c_{i}c_{1}}. This is 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 C_{1}, …, C_{k} of disjoint cycles with increasing 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 defined as

sign C = sign π := (−1)^{n+k}.

One can check that this coincides with the usual definition of the sign and parity of a permutation.

So we can rewrite (1) as follows:

det A = ∑_{cycle decompositions C} sign C · weight(C)  (3)

2.3 Clow Sequences

If we relax the requirement that all elements of a cycle are distinct, we arrive at the concept of a closed ordered walk or clow for short, in the terminology of Mahajan and Vinay [14]. A clow (c_{1}, c_{2}, …, c_{i}) of length i is a sequence of i ≥ 1 numbers where the head c_{1} of C is the unique smallest element: c_{1} < c_{2}, …, c_{i}.

A61:01 A61:02 A61:03 A61:04

A61:05

A61:06 A61:07 A61:08 A61:09 A61:10 A61:11 A61:12

A61:13

A61:14

A61:15 A61:16 A61:17 A61:18 A61:19 A61:20 A61:21

A61:22 A61:23 A61:24 A61:25

A61:26 A61:27 A61:28 A61:29 A61:30 A61:31 A61:32 A61:33 A61:34 A61:35 A61:36 A61:37 A61:38

The weight of C is again weight(C) = a_{c_{1}c_{2}}a_{c_{2}c_{3}} … a_{c_{i}c_{1}}. (Valiant [24], who introduced this concept, called this a c_{1}-loop.)

A clow sequence C is a sequence C_{1}, …, C_{k} of clows with a strictly increasing sequence of heads:

C = (\underline{c_{1}}, …), (\underline{c_{2}}, …), …, (\underline{c_{k}}, …),  (4)

with c_{1} < c_{2} < ⋯ < c_{k}. 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 is just (−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 = ∑_{clow sequences C of length n} sign C · weight(C)  (5)

Proof. It must be shown that all “bad” clow sequences which are not cycle decompositions cancel because the corresponding weights occur equally often with 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 is given in full detail in [14, 15].

Lemma 1. Let C be a clow sequence that contains repeated elements. Then we can find a clow sequence C̄ with the same weight and opposite sign. Moreover, this mapping between C and C̄ is an involution, i. e., if we construct the corresponding other clow sequence for C̄, we obtain C again.

Proof. Let C = C_{1}, …, C_{k}. We successively add the clows C_{k}, C_{k−1}, … until a repetition occurs. Assume that C_{m+1}, …, C_{k} is a set of disjoint cycles but C_{m}, C_{m+1}, …, C_{k} contains a repeated element, where 1 ≤ m ≤ k. Let C_{m} = (c_{1}, …, c_{i}), and let c_{j} be the first 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 C_{m+1}, …, C_{k}. Precisely one of these two cases will occur.

In case (A), we remove the cycle (c_{ℓ}=c_{j}, c_{ℓ+1}, …, c_{j−1}) from the clow C_{m} and 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 C_{m} after the element c_{j}. 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. □

A62:01

A62:02 A62:03 A62:04 A62:05 A62:06 A62:07

A62:08 A62:09 A62:10

A62:11 A62:12 A62:13 A62:14 A62:15 A62:16 A62:17 A62:18 A62:19 A62:20 A62:21 A62:22 A62:23 A62:24 A62:25 A62:26 A62:27 A62:28 A62:29 A62:30 A62:31 A62:32 A62:33 A62:34 A62:35 A62:36 A62:37 A62:38 A62:39 A62:40

A62:41

A62:42 A62:43

2.4 Dynamic Programming According to Length

In some sense, (1) and (3) are the most economical formulas for the determinant because all superfluous terms have been canceled and no further reduction is possible. However, we shall now see that the redundant form (5) is more amenable to an efficient algorithm. The idea is to construct clow sequences incrementally, arc by arc: We use dynamic programming to systematically compute the following quantities:

[l, c, c_{0}, 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 c_{0}, and the sign s = ±1 gives the parity of the number of finished clows so far.

A partial clow sequence is an initial piece of the sequence given in (4): It consists of a number of completed clows, and one incomplete clow. For completing the clow sequence, we only have to know the head and the current last vertex of the incomplete clow. In the end, we need to know the parity of the total number of clows; 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(n^{3}) nodes corresponding to the quantities [l, c, c_{0}, s], with 1 ≤ l ≤ n, 1 ≤ c_{0} ≤ c ≤ n, and s = ±1.

A partial clow sequence can be grown either by extending the current clow or by 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, c_{0}, s]: arcs to nodes [l + 1, c′, c_{0}, s] of cost a_{cc′}, for all c′ > c_{0}, and arcs of cost a_{cc_{0}}, to all nodes [l + 1, c′_{0}, c′_{0}, −s], for all c′_{0} > c_{0}. The latter type of arcs will start a new clow with head c′_{0} and no edges yet. We consider all nodes [0, c_{0}, c_{0}, 1] as start nodes and add two artificial end nodes [n, n + 1, n + 1, s] with s = ±1. The sum of 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(n^{4}): there are O(n) arcs leaving each of O(n^{3}) nodes. For example, one can calculate node values in the order or increasing l values. The storage requirement in this case is O(n^{2}) 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 computation of the determinant belongs to the complexity class GapL, a complexity class in 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, c_{0}, +1] and [l, c, c_{0}, −1] with opposite sign fields into one node [l, c, c_{0}] and introducing “negative” arc weights: the second type of arcs have then cost −a_{cc_{0}} instead of a_{cc_{0}}.

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 different factorization is possible: we can first consider the clow with

A63:01 A63:02 A63:03 A63:04

A63:05

A63:06

A63:07 A63:08 A63:09 A63:10 A63:11

A63:12

A63:13 A63:14 A63:15 A63:16 A63:17 A63:18 A63:19

A63:20

A63:21 A63:22 A63:23 A63:24

head c_{1} = 1, if such a clow is present, and the remaining clow sequence, which is a sequence on the ground set {2, …, n}. For expressing the clows with head 1 we partition the matrix A into its first row r, its first column s, and the remaining matrix M.

A = ( \begin{matrix}a_{11} & r \\ s & M \end{matrix} )

The sum of all clows of length l can be written as follows:

length l = 0: 1
length l = 1: a_{11}
length l = 2: ∑_{i=2}^{n} a_{1i}a_{i1} = rs
length l = 3: ∑_{i=2}^{n} ∑_{j=2}^{n} a_{1i}a_{ij}a_{j1} = rMs
⋮
length l:  rM^{l−2}s

The term 1 for length l = 0 corresponds to the “empty” clow consisting of no edges. It accounts for the case when the head of the first clow is larger than 1. We denote by q′_{i} the sum of signed weights of all clow-sequences of length n−1−i which don’t contain the element 1. Here we define the sign of a clow sequence with k clows as (−1)^{n−1+k}. Then we can write:

det A = −1 · q′_{−1} + a_{11} · q′_{0} + rs · q′_{1} + ⋯ + rM^{i−1}s · q′_{i} + ⋯  (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 all clow 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 = a_{11} · q′_{0} + rs · q′_{1} + ⋯ + rM^{i−1}s · q′_{i} + ⋯ + rM^{n−2}s · q′_{n−1}.

We intend to compute the values q′_{i} recursively, so we need to consider clow sequences of all possible lengths. Let us denote by q_{i} the sum of signed weights of all clow-sequences of length n − i (with the original definition of the sign). So q_{0} = det A, and q′_{0} = det M. By the above argument, we get the following

A64:01

A64:02 A64:03 A64:04

A64:05

A64:06

A64:07

A64:08 A64:09

A64:10 A64:11

A64:12

A64:13 A64:14 A64:15 A64:16 A64:17

A64:18

A64:19

A64:20

formula for computing q_{0}, q_{1}, q_{2}, …, q_{n} from the sequence q′_{0}, q′_{1}, …, q′_{n−1}.

( q_{n} ; q_{n−1} ; q_{n−2} ; ⋮ ; q_{2} ; q_{1} ; q_{0}) = (\begin{matrix}−1 & 0 & ⋯ & 0 & 0 & 0 \\ a_{11} & −1 & ⋯ & 0 & 0 & 0 \\ rs & a_{11} & ⋯ & 0 & 0 & 0 \\ ⋮ & ⋮ & ... & ⋮ & ⋮ & ⋮ \\ rM^{n−4}s & rM^{n−5}s & ⋯ & a_{11} & −1 & 0 \\ rM^{n−3}s & rM^{n−4}s & ⋯ & rs & a_{11} & −1 \\ rM^{n−2}s & rM^{n−3}s & ⋯ & rMs & rs & a_{11} \end{matrix}) (q′_{n−1} ; q′_{n−2} ; ⋮ ; q′_{2} ; q′_{1} ; q′_{0})  (7)

This matrix multiplication is nothing but a convolution between the sequence q′_{0}, q′_{1}, … and the sequence −1, a_{11}, rs, rMs, …, which can most conveniently be expressed as a polynomial multiplication:

q_{n}λ^{n} + q_{n−1}λ^{n−1} + ⋯ + q_{1}λ + q_{0}

= (−λ + a_{11} + rsλ^{−1} + rMsλ^{−2} + ⋯ + rM^{i}sλ^{−i−1} + ⋯)

× (q′_{n−1}λ^{n−1} + q′_{n−2}λ^{n−2} + ⋯ + q′_{1}λ + q′_{0}) (8)

Actually, these are not polynomials, but something like Laurent series in λ^{−1}. In explicit form, this means

( q_{n} ; q_{n−1} ; q_{n−2} ; ⋮ ; q_{2} ; q_{1} ; q_{0} ; 0 ; 0 ; ⋮) = (\begin{matrix}−1 & 0 & ⋯ & 0 & 0 & 0 & 0 & ⋯ \\ a_{11} & −1 & ⋯ & 0 & 0 & 0 & 0 & ⋯ \\ rs & a_{11} & ⋯ & 0 & 0 & 0 & 0 & ⋯ \\ ⋮ & ⋮ & ... & ⋮ & ⋮ & ⋮ & ⋮ & ... \\ rM^{n−4}s & rM^{n−5}s & ⋯ & −1 & 0 & 0 & 0 & ⋯ \\ rM^{n−3}s & rM^{n−4}s & ⋯ & a_{11} & −1 & 0 & 0 & ⋯ \\ rM^{n−2}s & rM^{n−3}s & ⋯ & rs & a_{11} & −1 & 0 & ⋯ \\ rM^{n−1}s & rM^{n−2}s & ⋯ & rMs & rs & a_{11} & −1 & ⋯ \\ rM^{n}s & rM^{n−1}s & ⋯ & rM^{2}s & rMs & rs & a_{11} & ⋯ \\ ⋮ & ⋮ & ... & ⋮ & ⋮ & ⋮ & ⋮ & ... \end{matrix}) (q′_{n−1} ; q′_{n−2} ; ⋮ ; q′_{1} ; q′_{0} ; 0 ; 0 ; ⋮)

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 · q′_{0} + rMs · q′_{1} + ⋯ + rM^{i}s · q′_{i} + ⋯ + rM^{n−1}s · q′_{n−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 q_{i}, after omitting irrelevant terms like q′_{−1} from (6), actually represent a subset of the clow-sequences that are characterized by the so-called prefix property.

2.6 The Characteristic Polynomial

We have associated the polynomial

P_{A}(λ) = q_{n}λ^{n} + q_{n−1}λ^{n−1} + ⋯ + q_{1}λ + q_{0}

A65:01 A65:02

A65:03

A65:04 A65:05 A65:06 A65:07 A65:08 A65:09 A65:10 A65:11 A65:12

A65:13

A65:14 A65:15

A65:16

A65:17 A65:18 A65:19 A65:20 A65:21 A65:22 A65:23 A65:24

A65:25 A65:26

A65:27

A65:28 A65:29

A65:30 A65:31 A65:32

A65:33

A65:34

A65:35

to the matrix A, and we have shown how P_{A}(λ) can be computed from the corresponding polynomial P_{M}(λ) of the (n − 1) × (n − 1) submatrix M:

P_{A}(λ) = (−λ + a_{11} + rsλ^{−1} + rMsλ^{−2} + ⋯ + rM^{i}sλ^{−i−1} + ⋯) · P_{M}(λ) (9)

This allows us to inductively compute polynomials for larger and larger submatrices, starting from P_{(a_{nn})}(λ) = a_{nn} − λ, until we finally obtain P_{A}(λ) and the determinant det A = q_{0}.

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), but also clow sequences of all other lengths. We have also found it convenient to use the corresponding sums of weights as coefficients of a polynomial P_{A}(λ). It turns out that this polynomial is nothing but the characteristic polynomial of the matrix A:

P_{A}(λ) := det(A − λI) = q_{n}λ^{n} + q_{n−1}λ^{n−1} + ⋯ + q_{1}λ + q_{0}  (10)

The following well-known formula generalizes formula (3) for the determinant to all coefficients of the characteristic polynomial.

q_{i} = ∑_{families C of disjoint cycles with length n − i} sign C · weight(C)  (11)

Here we consider not partitions into cycles, but families C of disjoint cycles with a total of n − i elements. The sign of such a family with k cycles is defined 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 q_{i} defined in the previous section are indeed the coefficients of the characteristic polynomial (10) as we know it. It is proved in the same way as Theorem 1, see [14].

Theorem 2. The coefficients q_{i} (i = 0, 1, …, n) of the characteristic polynomial of A can be expressed as

q_{i} = ∑_{clow sequences C of length n − i} sign C · weight(C),  (12)

where the sign of a clow sequence of length n − i with k cycles is defined as (−1)^{n−i+k}. □

The dynamic programming algorithm of Section 2.4 can easily be extended to compute all coefficients q_{i}: Simply adding artificial 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

A = ( \begin{matrix}a_{11} & r \\ s & M \end{matrix} )  (13)

A66:01

A66:02

A66:03 A66:04 A66:05 A66:06 A66:07 A66:08

A66:09

A66:10 A66:11

A66:12

A66:13

A66:14

A66:15 A66:16 A66:17

A66:18

A66:19 A66:20 A66:21 A66:22 A66:23 A66:24 A66:25 A66:26 A66:27 A66:28 A66:29

can be expressed as

det A = det M · (a_{11} − rM^{−1}s)  (14)

if the matrix M is invertible. This follows from the fact that det M is the cofactor of a_{11}, and hence (A^{−1})_{11} = det M/ det A. The element (1, 1) of the inverse A can be determined by block Gauss-Jordan elimination of (13), starting with the lower right block M. This gives (A^{−1})_{11} = (a_{11} −rM^{−1}s)^{−1} which leads to (14).

Now we wish to express P_{A}(λ) = det(A − λI) in terms of P_{M}(λ) = det(M − λI). Applying (14) to the matrix A − λI gives

P_{A}(λ) = |A − λI| = | \begin{matrix}a_{11} − λ & r \\ s & M − λI \end{matrix} |
= |M − λI| · (a_{11} − λ − r(M − λI)^{−1}s)
= P_{M}(λ) · (a_{nn} − λ + r(λI − M)^{−1}s) =: P_{M}(λ) · F (λ)

For simplifying the expression F (λ), which is the multiplier which turns P_{M} into P_{A}, 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

\begin{aligned}F (λ) & = a_{11} − λ + r(λI − M)^{−1}s \\ & = −λ + a_{11} + rsλ^{−1} + rMsλ^{−2} + rM^{2}sλ^{−3} + rM^{3}sλ^{−4} + ⋯, \end{aligned} (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 first combinatorial interpretation of (7), referred to the algorithm as the Samuelson-Berkowitz algorithm. However, Samuelson’s algorithm for computing the characteristic polynomial is an O(n^{3}) algorithm that uses divisions. This is also the way in which it is described in [8, §45]. So the above algorithm can definitely not be attributed to Samuelson. However, a passage towards the end of Section §45 in [8] contains the system (7). Faddeyev and Faddeyeva state this as a possible stricter justification of Samuelson’s method. The mention briefly that this formula follows from

A67:01 A67:02 A67:03 A67:04 A67:05 A67:06 A67:07 A67:08

A67:09

A67:10 A67:11

A67:12

A67:13 A67:14 A67:15

A67:16

A67:17

A67:18

A67:19 A67:20

A67:21

A67:22 A67:23 A67:24 A67:25 A67:26 A67:27 A67:28 A67:29 A67:30 A67:31 A67:32 A67:33

the consideration of bordered determinants. This calculation has been sketched above in Section 2.7.

The fact that (7) appears in the section about Samuelson’s method in the book [8] may be the reason for the erroneous attribution. It seems that the method ought to be credited to Faddeyev and Faddeyeva. However, being numerical analysts, they would probably not have considered the implied O(n^{4}) procedure as a good “method”. The first algorithmic use of (7) as a recursive formula is apparently due to Berkowitz [3].

3 The Pfaffian

The determinant of a skew-symmetric matrix (a_{ij} = −a_{ji}) is the square of another expression, which is called the Pfaffian [12, 16]. For example,

| \begin{matrix}0 & a_{12} & a_{13} & a_{14} \\ −a_{12} & 0 & a_{23} & a_{24} \\ −a_{13} & −a_{23} & 0 & a_{34} \\ −a_{14} & −a_{24} & −a_{34} & 0 \end{matrix} | = (a_{12}a_{34} − a_{13}a_{24} + a_{14}a_{23})^{2}

The determinant of a skew-symmetric matrix of odd order is always 0. Formally, the Pfaffian of a skew-symmetric matrix A with an even number of rows and columns can be defined as follows:

pf A = ∑_{perfect matchings M} sign M · weight(M)  (16)

Here, a perfect matching M with k = n/2 edges is written as

M = (\underline{i_{1}}, j_{1}), (\underline{i_{2}}, j_{2}), …, (\underline{i_{k}}, j_{k}),  (17)

where, by convention, i_{l} < j_{l} for each l. The sign of the matching M is defined as the sign of

(\begin{matrix}1 & 2 & 3 & 4 & … & n − 1 & n \\ i_{1} & j_{1} & i_{2} & j_{2} & … & i_{k} & j_{k} \end{matrix})

when this is regarded as a permutation of {1, …, n}. The weight of M is a_{i_{1}j_{1}}a_{i_{2}j_{2}} … a_{i_{k}j_{k}}. One can check that the sign of M does not depend on the order in which the edges are given. Moreover, exchanging a pair i_{l}, j_{l} changes the sign but not the signed weight since a_{ji} = −a_{ij}. This means that we are in fact free to choose any convenient permutation π as a representation of a given matching M.

Because of their close connection with matchings, Pfaffians are of interest in combinatorics, see [12, 21]. For example, for certain graphs, including planar graphs, it is possible to count the number of perfect matchings in polynomial time by using Pfaffians. Knuth [10] gives a short history of Pfaffians, and argues that Pfaffians are in some way more fundamental than determinants, to which they are closely related (see Section 4).

A68:01 A68:02 A68:03 A68:04

A68:05

A68:06 A68:07 A68:08 A68:09 A68:10

A68:11

A68:12 A68:13 A68:14 A68:15 A68:16 A68:17 A68:18 A68:19 A68:20 A68:21

A68:22

A68:23 A68:24 A68:25 A68:26 A68:27

A68:28

A68:29

A68:30 A68:31 A68:32 A68:33 A68:34 A68:35 A68:36 A68:37

A68:38

The Pfaffian can be computed in O(n^{3}) steps by an elimination procedure that is similar to Gaussian elimination for symmetric matrices. Alternatively, the square root of the determinant gives the Pfaffian up to the correct sign. Both approaches 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 Pfaffians as well. We consider the product of the Pfaffians of two skew-symmetric matrices A and B. Expanding the definitions (16) leads to an overlay of two matchings M_{A} and M_{B}. The union of M_{A} with M_{B} decomposes into a disjoint union of alternating cycles with a total of n edges:

M_{A} ⊎ M_{B} = C = C_{1} ∪ C_{2} ∪ ⋯

An alternating cycle is a cycle of even length whose edges alternate between M_{A} and M_{B}. Edges which are contained both in M_{A} and M_{B} must be kept separate the union M_{A} ⊎ M_{B}; for each “common” edge, M_{A} ⊎ M_{B} contains two (distinguishable) parallel edges, forming an alternating cycle of length 2.

An alternating cycle can be traversed in two orientations. We specify the orientation by taking the vertex with the smallest number as the “first vertex” or head of the cycle and insisting that the first edge should belong to M_{A}.

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 = (c_{1}, c_{2}, …, c_{i}) with c_{1} < c_{2}, …, c_{i} is

weight(C) = a_{c_{1}c_{2}}b_{c_{2}c_{3}}a_{c_{3}c_{4}}b_{c_{4}c_{5}} … a_{c_{i−1}c_{i}}b_{c_{i}c_{1}}.  (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 defined as sign C = (−1)^{k}.

Theorem 3.

pf A · pf B = ∑_{alternating cycle decompositions C} sign C · weight(C)

Proof. The Pfaffian has been defined as a weighted sum of matchings; so the left side is a sum over all pairs of matchings (M_{A}, M_{B}), where M_{A} is regarded as a matching for A and M_{B} 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 M_{A} ⊎ M_{B} decomposes into k disjoint cycles C = (c_{1}, c_{2}, …, c_{i}), (d_{1}, d_{2}, …, d_{j}), . . .. We may conveniently represent M_{A} and M_{B} by the permutations

π_{A} = \binom{1 2 3 4 ⋯ i − 1 i i + 1 i + 2 ⋯ i + j i + j + 1 ⋯ n}{c_{1} c_{2} c_{3} c_{4} ⋯ c_{i−1} c_{i} d_{1} d_{2} ⋯ d_{j} ⋯}

A69:01

A69:02

A69:03 A69:04 A69:05

A69:06 A69:07 A69:08 A69:09

A69:10 A69:11 A69:12 A69:13 A69:14 A69:15 A69:16 A69:17 A69:18

A69:19 A69:20 A69:21

A69:22 A69:23 A69:24 A69:25

A69:26

A69:27 A69:28 A69:29 A69:30 A69:31 A69:32

and

π_{B} = \binom{1 2 3 4 ⋯ i − 1 i i + 1 i + 2 ⋯ i + j i + j + 1 ⋯ n}{c_{2} c_{3} c_{4} c_{5} ⋯ c_{i} c_{1} d_{2} d_{3} ⋯ d_{1} ⋯} ,

respectively. This gives sign M_{A} · sign M_{B} = 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 Pfaffian 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 difference 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 determinant of 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 the cycle decompositions containing odd cycles cancel. □

If we want to use Theorem 3 to get an expression for the Pfaffian of a single matrix A, we must set B to a matrix B_{0} whose Pfaffian is 1. So we select the skew-symmetric matrix

B_{0} = (\begin{matrix}0 & 1 \\ −1 & 0 \\ & & 0 & 1 \\ & & −1 & 0 \\ & & & & . & . & . \\ & & & & & & & 0 & 1 \\ & & & & & & & −1 & 0 \end{matrix})  (19)

that represents the matching M_{0} = (1, 2), (3, 4), …, (n − 1, n), which we use as a reference matching. So we immediately get the following alternate expression for the Pfaffian. Corollary 2.

pf A = ∑_{M_{0}-alternating cycle decompositions C} sign C · weight(C)  □

Here, the general concepts of alternating cycles, cycles families, etc., specialize as follows: An M_{0}-alternating cycle (c_{1}, c_{2}, …, c_{i}) is a cycle of even length i. In addition, we must have c_{3} = c_{2} − 1 if c_{2} is even and c_{3} = c_{2} + 1 if c_{2} is odd, and the same relation holds between c_{4} and c_{5}, …, c_{i} and c_{1}. Its weight is given by (18) with respect to B_{0}, i. e., b_{2j−1,2j} = 1 and b_{2j,2j−1} = −1. All other notions are unchanged.

A70:01

A70:02 A70:03 A70:04

A70:05

A70:06

A70:07

A70:08 A70:09 A70:10 A70:11 A70:12 A70:13 A70:14 A70:15 A70:16 A70:17 A70:18 A70:19 A70:20 A70:21 A70:22 A70:23

A70:24 A70:25

A70:26

A70:27 A70:28 A70:29 A70:30 A70:31 A70:32 A70:33

A70:34

3.2 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 = ∑_{alternating clow sequences C of length n} sign C · weight(C) (20)

pf A = ∑_{M_{0}-alternating clow sequences C of length n} sign C · weight(C)  (21)

Here, an alternating clow (c_{1}, c_{2}, …, c_{i}) of length i is defined like a clow, with the additional restriction that i must be even. The condition on the clow head is unchanged: c_{1} < c_{2}, …, c_{i}. The weight of alternating clows is defined like in (18) for alternating cycles. An M_{0}-alternating clow (c_{1}, c_{2}, …, c_{i}) must have every other edge (c_{2}, c_{3}), (c_{4}, c_{5}), …, (c_{i}, c_{1}) belonging to M_{0}. Alternating clow sequences, their length, weight, and sign are defined just as for clows and alternating cycle families. (Our notation and terminology differs from [13].)

The proof of (20), and of (21) which is a special case of it, follows the lines of the proof of Lemma 1: a sign-reversing bijection among “bad” alternating clow sequences which contain repeated elements is established. The only change occurs in case (A), when a clow C_{k} forms a sub-cycle (c_{ℓ}=c_{j}, c_{ℓ+1}, …, c_{j−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 C_{k} by (c_{ℓ}=c_{j}, c_{j−1}, …, c_{ℓ+1}). The weight of C_{k} changes sign, and it can also be checked easily that the resulting mapping between bad alternating clow sequences remains an involution. □

Mahajan, Subramanya, and Vinay [13] considered only the case (21) of a single matrix, and they gave a slightly different proof.

3.3 Incremental Algorithms for the Pfaffian

A dynamic programming algorithm for evaluating (21) by summing over all partial alternating clow sequences in order of increasing length can easily be formulated 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 Pfaffians only make sense for matrices of even order, we partition the given skew-symmetric matrix by splitting two rows and columns from it.

A = ( \begin{matrix}0 & a_{12} & r \\ −a_{12} & 0 & −s^{T} \\ −r^{T} & s & M \end{matrix} )

A71:01 A71:02

A71:03

A71:04 A71:05 A71:06 A71:07 A71:08 A71:09

A71:10

A71:11 A71:12 A71:13 A71:14 A71:15 A71:16 A71:17 A71:18 A71:19 A71:20

A71:21

A71:22 A71:23

A71:24

A71:25 A71:26

A71:27 A71:28 A71:29 A71:30 A71:31

A71:32

A71:33 A71:34

Enumeration of alternating clow sequences with head c_{1} = 1 according to length yields the power series

G(λ) = −λ^{2} + a_{12} + rB_{0}sλ^{−2} + rB_{0}MB_{0}sλ^{−4} + rB_{0}MB_{0}MB_{0}sλ^{−6} + ⋯(22)

B_{0} is the “skew-symmetric unit matrix” (19) of appropriate dimension, here (n − 2) × (n − 2). The coefficient 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), and hence the factor b_{21} = −1 is alway present in the above terms. We define the Pfaffian-characteristic polynomial of a skew-symmetric n × n matrix A as

P̃_{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 M_{0}-alternating clow-sequences of length n − i. This definition can also be used for odd n. Similarly as in Theorem 2, q̃_{i} is equal to the sum of signed weights of all M_{0}- alternating cycle decompositions of length n − i. It can also be interpreted the sum of signed weights of certain matchings with (n − i)/2 edges as in (16), with an appropriate definition 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 equation

P̃_{A}(λ) = G(λ) · P̃_{M}(λ),  (23)

which immediately leads to a recursive division-free algorithm for computing P̃_{A}(λ) and the Pfaffian pf A = q̃_{0}.

4 Open Questions

The relation

det A = pf (\begin{matrix}0 & A \\ −A^{T} & 0 \end{matrix})

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

4.1 The Pfaffian-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

A72:01 A72:02 A72:03 A72:04 A72:05 A72:06 A72:07 A72:08 A72:09 A72:10 A72:11 A72:12 A72:13 A72:14 A72:15

A72:16 A72:17 A72:18 A72:19 A72:20 A72:21 A72:22 A72:23 A72:24 A72:25 A72:26 A72:27

A72:28

A72:29 A72:30 A72:31 A72:32

A72:33

A72:34 A72:35 A72:36 A72:37

order n. In contrast to the case of the determinant, we have no idea what an algebraic 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}(\sqrt{(}λ)), as the relation q_{0} = q̃_{0}^{2} between their constant terms would suggest.

In Theorems 3 and 4, it has turned out that it appears more natural to deal with the product of Pfaffians of two different matrices A and B: after all, Corollary 2 is only a special case of Theorem 3 which is by no means easier to formulate, and the same relation holds between (20) and (21) in Theorem 4. Many identities about Pfaffians involve products of two Pfaffians that are taken from submatrices of a single given matrix. One might speculate that taking two different 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 B_{0} from (19) into (20) of Theorem 4, one could also use the matrix

B_{1} = (\begin{matrix}0 & 1 & 1 & 1 & ⋯ & 1 \\ −1 & 0 & 1 & 1 & ⋯ & 1 \\ −1 & −1 & 0 & 1 & ⋯ & 1 \\ −1 & −1 & −1 & 0 & ⋯ & 1 \\ ⋮ & ⋮ & ⋮ & ⋮ & ... & ⋮ \\ −1 & −1 & −1 & −1 & ⋯ & 0 \end{matrix}) ,

which might be called the skew-symmetric analog of the all-ones matrix and has pf B_{1} = 1. This results in a different expression for pf A than (21) in Theorem 4 and to another algorithm. That algorithm is in a way simpler than the incremental 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 B_{1} also leads to a different concept of a Pfaffian-characteristic polynomial.

There is an analogous algorithm for directly computing the product of two Pfaffians according to equation (20) of Theorem 4. It seems that the Pfaffiancharacteristic 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}(λ) = P_{A}(λ),

for skew-symmetric matrices A.

Is there a way to introduce a second matrix B also into the computation of determinants 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 coefficients of this polynomial?

A73:01

A73:02 A73:03 A73:04 A73:05 A73:06 A73:07 A73:08 A73:09 A73:10 A73:11 A73:12 A73:13 A73:14 A73:15 A73:16 A73:17 A73:18 A73:19 A73:20 A73:21 A73:22 A73:23 A73:24 A73:25 A73:26 A73:27 A73:28 A73:29 A73:30 A73:31 A73:32

A73:33

A73:34 A73:35 A73:36 A73:37 A73:38 A73:39 A73:40 A73:41 A73:42

4.2 Complexity

All division-free methods for computing the determinant or the Pfaffian that we have described take O(n^{4}) 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(n^{3}) time for computing the coefficients of the factor F (λ) or G(λ) respectively, and O(n^{2}) steps for multiplying the two polynomials, if this is done in the straightforward way. This amounts to a total of O(n^{4}) 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, each matrix-vector product will take only O(m) steps instead of O(n^{2}) (assuming m ≥ n). This leads to an algorithm with a running time of O(mn^{2}). Thus, when m is small, this algorithm may be preferable to ordinary Gaussian elimination even 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(n^{4}) to O(mn^{2}) steps. It is a challenging problem to devise a division-free O(n^{3}) method for computing the determinant.

The algorithms discussed in this paper can also be used to design good parallel algorithms for determinants and Pfaffians [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 clearly faster already for matrices of order 40. For a simple 5 × 5 matrix with rational entries like (a_{ij}) = (\frac{1}{x_{i}+x_{j}}), our algorithm for the characteristic polynomial was 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 unless fractions are reduced after every step. However, to avoid the excessive growth of 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 ex-

A74:01 A74:02 A74:03 A74:04 A74:05 A74:06 A74:07 A74:08 A74:09 A74:10 A74:11 A74:12 A74:13 A74:14

A74:15

A74:16 A74:17 A74:18 A74:19 A74:20 A74:21 A74:22 A74:23 A74:24 A74:25 A74:26 A74:27 A74:28 A74:29 A74:30 A74:31 A74:32 A74:33 A74:34 A74:35 A74:36 A74:37 A74:38 A74:39 A74:40 A74:41 A74:42 A74:43 A74:44 A74:45 A74:46 A74:47

plode and hence, this algorithm runs in polynomial time. The same approach can be extended to Pfaffians [9].

Thus, for evaluating determinants or Pfaffians of integers, this approach is superior to our division-free algorithms because it takes only O(n^{3}) 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 final result into an algorithm without divisions. Divisions are simulated by expanding expressions into truncated power series. This is somehow reminiscent of the recursion (9) which also involves power series that can be truncated because the result is a polynomial. According to [4], Strassen’s approach leads to an algorithm with O(n^{5}) steps in the case of Gaussian elimination. It would be interesting to see if there is any relation of this algorithm to the algorithms presented here.

References

1. Martin Aigner, Lattice paths and determinants. In: Computational Discrete Mathematics, ed. Helmut Alt, (this volume), June 2000, pp. 1–12.

2. Erwin H. Bareiss, Sylvester’s identity and multistep integer-preserving Gaussian elimination. Math. Comput. 22 (1968), 565–578.

3. Stuart J. Berkowitz, On computing the determinant in small parallel time using a small 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 numérique des équations algébriques, tome II: systèmes des plusieurs é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. opér. 31 (1997), 203–209.

8. D. K. Faddeyev, V. N. Faddeyeva, Vyč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 Maffioli, On the computation of pfaffians. Discr. Appl. Math. 51 (1994), 269–275.

10. Donald E. Knuth, Overlapping Pfaffians. 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. Séminaire Lotharingien de Combinatoire B42q (1999), 67 pp.

12. László Lová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 Pfaffians. In: Computing and combinatorics. Proc. 5th annual international conference. (COCOON ’99), Tokyo, July 1999, ed. Takao Asano et al., Lecture Notes Comput. Sci. 1627, Springer-Verlag, pp. 134–143 (1999). Extended version: DIMACS Technical Report 99-39, Rutgers University, July 1999.

A75:01 A75:02 A75:03 A75:04 A75:05 A75:06 A75:07 A75:08 A75:09 A75:10 A75:11 A75:12 A75:13 A75:14 A75:15 A75:16 A75:17 A75:18 A75:19 A75:20 A75:21 A75:22 A75:23 A75:24 A75:25

14. Meena Bhaskar Mahajan, V Vinay, Determinant: Combinatorics, algorithms, and complexity. 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. Gü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 coefficients of the characteristic equation. Ann. Math. Statist. 13 (1942), 424–429.

20. Dennis Stanton, Dennis White, Constructive combinatorics. Springer-Verlag, New York 1986.

21. John R. Stembridge, Nonintersecting paths, pfaffians, 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 difficult? 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.