The Butera-Pernici algorithm (1)

Introduction

We know that there is no general sub-exponential algorithm for computing the permanent of a square matrix. But we may very reasonably ask -- might there be a faster, possibly even polynomial-time algorithm, for some specific classes of matrices? For example, a sparse matrix will have most terms of the permanent zero -- can this be somehow leveraged for a better algorithm?

The answer seems to be a qualified "yes". In particular, if a matrix is banded, so that most diagonals are zero, then a very fast algorithm can be applied. This algorithm is described in an online article by Paolo Butera and Mario Pernici called Sums of permanental minors using Grassmann algebra. Accompanying software (Python programs) is available at github. This software has been rewritten for the SageMath system, and you can read about it in the documentation. The algorithm as described by Butera and Pernici, and as implemented in Sage, actually produces a generating function.

Our intention here is to investigate a simpler version, which computes the permanent only.

Basic outline

Let \(M\) be an \(n\times n\) square matrix, and consider the polynomial ring on \(n\) variables \(x_1,x_2,\ldots,x_n\). Each row of the matrix will correspond to an element of this ring; in particular row \(i\) will correspond to

\[ \sum_{j=1}^nm_{ij}a_j=m_{i1}a_1+m_{i2}a_2+\cdots+m_{in}a_n. \]

Suppose further that all the generating elements \(x_i\) are nilpotent of order two, so that \(x_i^2=0\).

Now if we take all the row polynomials and multiply them, each term of the product will have order \(n\). But by nilpotency, all terms which contain a repeated element will vanish. The result will be only those terms which contain each generator exactly once, of which there will be \(n!\). To obtain the permanent all that is required is to set \(x_i=1\) for each generator.

Here's an example in Sage.

1sage: R.<a,b,c,d,e,f,g,h,i,x1,x2,x3> = PolynomialRing(QQbar)
2sage: M = matrix([[a,b,c],[d,e,f],[g,h,i]])
3sage: X = matrix([[x1],[x2],[x3]])
4sage: MX = M*X
5
6[a*x1 + b*x2 + c*x3]
7[d*x1 + e*x2 + f*x3]
8[g*x1 + h*x2 + i*x3]

To implement nilpotency, it's easiest to reduce modulo the ideal defined by \(x_i^2=0\) for all \(i\). So we take the product of those row elements, and reduce:

1sage: I = R.ideal([x1^2, x2^2, x3^2])
2sage: pr = MX[0,0]*MX[1,0]*MX[2,0]
3sage: pr.reduce(I)
4
5c*e*g*x1*x2*x3 + b*f*g*x1*x2*x3 + c*d*h*x1*x2*x3 + a*f*h*x1*x2*x3 + b*d*i*x1*x2*x3 + a*e*i*x1*x2*x3

Finally, set each generator equal to 1:

1sage: pr.reduce(I).subs({x1:1, x2:1, x3:1})
2
3c*e*g + b*f*g + c*d*h + a*f*h + b*d*i + a*e*i

and this is indeed the permanent for a general \(3\times 3\) matrix.

Some experiments

Let's experiment now with the matrices we've seen in a previous post, which contain three consecutive super-diagonals of ones, and the rest zero.

Such a matrix is easy to set up in Sage:

 1sage: n = 10
 2sage: v = n*[0]; v[1:4] = [1,1,1]
 3sage: M = matrix.circulant(v)
 4sage: M
 5
 6[0 1 1 1 0 0 0 0 0 0]
 7[0 0 1 1 1 0 0 0 0 0]
 8[0 0 0 1 1 1 0 0 0 0]
 9[0 0 0 0 1 1 1 0 0 0]
10[0 0 0 0 0 1 1 1 0 0]
11[0 0 0 0 0 0 1 1 1 0]
12[0 0 0 0 0 0 0 1 1 1]
13[1 0 0 0 0 0 0 0 1 1]
14[1 1 0 0 0 0 0 0 0 1]
15[1 1 1 0 0 0 0 0 0 0]

Similarly we can define the polynomial ring:

1sage: R = PolynomialRing(QQbar,x,n)
2sage: R.inject_variables()
3
4Defining x0, x1, x2, x3, x4, x5, x6, x7, x8, x9

And now the polynomials corresponding to the rows:

 1sage: MX = M*matrix(R.gens()).transpose()
 2sage: MX
 3
 4[x1 + x2 + x3]
 5[x2 + x3 + x4]
 6[x3 + x4 + x5]
 7[x4 + x5 + x6]
 8[x5 + x6 + x7]
 9[x6 + x7 + x8]
10[x7 + x8 + x9]
11[x0 + x8 + x9]
12[x0 + x1 + x9]
13[x0 + x1 + x2]

If we multiply them, we will end up with a huge expression, far too long to display:

1sage: pr = prod(MX[i,0] for i in range(n))
2sage: len(pr.monomials)
3
414103

We could reduce this by the ideal, but that would be slow. Far better to reduce after each separate multiplication:

1sage: I = R.ideal([v^2 for v in R.gens()])
2sage: p = R.one()
3sage: for i in range(n):
4              p = p*MX[i,0]
5              p = p.reduce(I)
6
7sage: p.subs({v:1 for v in R.gens())
8125

The answer is almost instantaneous. We can repeat the above list of commands starting with different values of n; for example with n=20 the result is 15129, as we expect.

This is not yet optimal; for n=20 on my machine the final loop takes about 7.8 seconds. Butera and Pernici show that the multiplication and setting the variables to one can sometimes be done in the opposite order; that is, some variables can be identified to be set to one before the multiplication. This can speed the entire loop dramatically, and this optimization has been included in the Sage implementation. For details, see their paper.