Compute Permanent of a Matrix with Ryser's Algorithm

Introduction

This article introduces the permanent of a matrix (as a function mapping matrices to real numbers) and in particular a straightforward implementation of Ryser's formula to compute the exact permanent in O(n*2^n).

Although permanents are defined similarly to the determinant, the latter being computable in T(n³) (see my CSML project on CodeProject), no polynomial algorithm for permanent computation is known so far.

Furthermore, the computation of a permanent of a 0-1-matrix A is #P-complete, and a polynomial algorithm for permanents would particularly imply P = NP. However, approximation is well possible: "When the entries of A are nonnegative, the permanent can be computed approximately in probabilistic polynomial time, up to an error of eM, where M is the value of the permanent and e > 0 is arbitrary." (See Wiki article.)

Definition

Let

be a real n by n matrix, then the permanent is defined by

(Leibniz's formula), where

For instance, consider the n by n 1-matrix U (all entries are ones), then we receive per(U) = n! -- why? There are n! permutations of {1,...,n}, and each permutated product in Leibniz's formula is equal to one.

Complexity of Leibniz's formula: (n-1)*(n!) flops, which is super-exponential.

Sadly, there is no geometric interpretation of what a permanent is (determinants can be considered an n-dimensional volume); its use is mainly in combinatorics and graph theory: The permanent is a count for the number of perfect matchings in a bipartite graph (ibd).

Moreover, there is a recursive algorithm similar to the Laplacian development of the determinant. You can develop by the i-th row

or by the j-th column

where the i-j-minor of A is defined by

This algorithm is still super-exponential, but it is well suited to calculate permanents of very small matrices by hand.

Ryser's Formula

The fastest known exact algorithm for the computation of the permanent is due to Ryser:

where

is the set of all n by k submatrices of A, and

is the i-th row sum of X. Ideally, Ryser's algorithm needs (n-1)(2^n - 1) flops; the implementation presented here runs in O((n^2)(2^n)).

The Code

At first, we need to transform a decimal to a binary array with an additional element (the last one) saving the number of ones in the array:

inline int* dec2binarr(long n, int dim)
{
    // note: res[dim] will save the sum res[0]+...+res[dim-1]
    int* res = (int*)calloc(dim + 1, sizeof(int));   
    int pos = dim - 1;

    // note: this will crash if dim < log_2(n)...
    while (n > 0)
    {
        res[pos] = n % 2;
        res[dim] += res[pos];
        n = n / 2; // integer division        
        pos--;
    }

    return res;
}

Here comes the actual algorithm:

long per(int* A, int n) // expects n by n matrix encoded as vector
{
    long sum = 0;
    long rowsumprod, rowsum;
    int* chi = new int[n + 1];    
    double C = (double)pow((double)2, n); 

    // loop all 2^n submatrices of A
    for (int k = 1; k < C; k++)
    {
        rowsumprod = 1;
        chi = dec2binarr(k, n); // characteristic vector        

        // loop columns of submatrix #k
        for (int m = 0; m < n; m++)
        {
            rowsum = 0;

            // loop rows and compute rowsum
            for (int p = 0; p < n; p++)
                rowsum += chi[p] * A[m * n + p];
        
            // update product of rowsums
            rowsumprod *= rowsum;    
        
            // (optional -- use for sparse matrices)
            // if (rowsumprod == 0) break;    
        }        
        
        sum += (long)pow((double)-1, n - chi[n]) * rowsumprod;        
    }    

    return sum;
} 

Notes

1. The submatrices are chosen by use of a characteristic vector chi (only the columns are considered, where chi[p] == 1). To retrieve the t from Ryser's formula, we need to save the number n-t, as is done in chi[n]. Then we get t = n - chi[n].
2. The matrix parameter A is expected to be a one-dimensional integer array -- should the matrix be encoded row-wise or column-wise? -- It doesn't matter. The permanent is invariant under row-switching and column-switching, and it is .
3. Further enhancements: If any rowsum equals zero, the entire rowsumprod becomes zero, and thus the m-loop can be broken. Since if-statements are relatively expensive compared to integer operations, this might save time only for sparse matrices (where most entries are zeros).
4. If anyone finds a polynomial algorithm for permanents, he will get rich and famous (at least in the computer science world).

History

• 13^th November, 2007: Initial post