Using the Johnson-Gentleman Tree Minor Algorithm to Calculate the Permanent

Algorithm Description

The elements of A320843 can be calculated as the permanent of a (1-based) matrix \(A = a_{ij}\) where \(a_{ij} = \begin{cases} 1& \text{ if } i | j \vee j | i\\ 0& \text{ otherwise} \end{cases}\)

Considering zero-based indices in the following.

With the permanent of a square \(n \times n\) matrix \(A\) defined as

and Laplace’s expansion by minors

where \(M_{ij}\) is the submatrix of \(A\) where row \(i\) and column \(j\) are removed, a recursive approach is easy to implement, but of no practical use.

Some definitions:

Let \(I_n = \{0,\ldots,n-1\}\) be the set of row and column indices of a \(n \times n\) matrix,

Let us consider minors of \(A\) that are formed from the first rows of \(A\).

A subset \(m \subseteq I_n\) having \(|m| = k, 0 \le k \le n\) can be uniquely identified with a minor that is build from the first \(k\) rows of \(A\) and the column set \(m\). Such minor will be called a \(k\)-minor.

Furthermore, let \(p_m = \perm(m)\) the permanent of the minor \(m\).

\(m^{(k)}\) denotes an arbitrary \(k\)-minor.

In their paper, Johnson and Gentleman investigated an approach starting with trivial 1-minors: the matrix’ first row. Then the permanents of \(k-1\)-minors are combined to form \(k\)-minors. They proved that this approach, though still having a high computational complexity, is optimal when the permanent is calculated by expanding minors. (The proof is about determinants.)

Let \(\C(m) = I_n \setminus m\) be the complement set of \(m\) and \(M_k\) be the set of all \(k\)-minors.

Then for every \(k\)-minor \(m^{(k)} \in M_k\), \(k \lt n\) and every \(c \in \C(m^{(k)})\) there is a \(k+1\)-minor \(m^{(k+1)} = m^{(k)} \cup \{c\}\) that contains the summand \(a_{k, c}\, p_{m^{(k)}}\) in its expansion and all minors in its expansion are build from the extension of one of the \(k\)-minors.

This can be described by the following formula:

This is the heart of the algorithm, but it does not really describe an actual implementation. That would consist of three nested loops:

• loop over \(k, 0 \le k \lt n\)
• loop over all \(k\)-minors
• loop over the minor’s complement columns

The inner loop produces the terms \(a_{k, c_i}\, p_{m_i^{(k)}}\) that are accumulated in the slot of the corresponding \(k+1\)-minor \(m_i^{(k)} \cup \{c_i\}\).

This approach can be implemented with a number \(m, 0 \le m \lt 2^n\) representing a minor \(m\) where the number has a 1-bit at position \(j\) in its binary representation when the corresponding minor contains the column \(j\).

Starting with the single 0-minor \(\emptyset\) and \(p_\emptyset = 1\), we find the 1-minors as the elements of the first row and after \(n\) steps there is only one minor left: \(m^{(n)} = I_n\) with \(p_{I_n} = \perm(A)\).

With sufficient memory, this algorithm can be faster than Ryser’s formula, but for larger general matrices it becomes impractical due to its high memory consumption.

Two measures can be taken to make the Johnson-Gentleman algorithm usable for A320843:

Keep the list of minors minimal

A minor \(m^{(k)}\) having \(p_{m^{(k)}} = 0\) will not contribute to any larger minor. It can therefore be removed from the set of \(k\)-minors.

Arrange non-zero items

Consider a triangular matrix. Its permanent is the product of its diagonal elements. The algorithm as described above applied to an upper triangular matrix is very inefficient as it generates lots of minors that do not contribute to the final result.

On the other hand, when applied to a lower triangular matrix, there is only one \(k\)-minor for every \(k\) and the final result can be found in a few steps.

Therefore it is crucial that a matrix is transformed into a form that suits best to the algorithm.

The permanent of a matrix does not change when two rows or columns are exchanged and these may be arranged in a way that resembles a lower triangular matrix the best it can.

Sorting the columns and the rows lexicographically by zero/nonzero approaches a triangular form.

To get an impression of the preconditioning, here are the original and the conditioned matrices for \(n=20\):

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1
1 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0
1 1 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1
1 0 0 0 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1
1 1 1 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0
1 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0
1 1 0 1 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0
1 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0
1 1 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1
1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0
1 1 1 1 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0
1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0
1 1 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0
1 0 1 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0
1 1 0 1 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0
1 1 1 0 0 1 0 0 1 0 0 0 0 0 0 0 0 1 0 0
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0
1 1 0 1 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1

0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 1
0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 1
0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 0 0 1 1
0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1
0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 1 1 1 1
0 0 0 0 0 0 0 1 0 0 0 0 0 1 1 0 1 0 1 1
0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 1 0 1
0 0 0 0 0 0 1 1 0 0 0 0 0 0 1 0 1 1 0 1
0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 1 1 1 1
0 0 0 0 0 1 0 0 0 0 1 0 0 0 1 1 0 0 0 1
0 0 0 0 0 1 0 0 0 0 1 0 0 1 1 1 0 0 0 1
0 0 0 0 0 1 0 0 1 0 1 0 0 0 0 1 0 0 0 1
0 0 0 0 0 1 0 1 1 0 0 0 0 0 0 0 0 0 0 1
0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1
0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1
0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

Implementations

There are two program versions, producing identical results up to their limits.

J

Due to the nature of J, this is not too easy to comprehend. The core code calculating the permanent is:

permanent_jg =: 3 : 0
getrow =. {~ +/@,@#:@((<0 0)&{)   NB. count the bits in the first selector to identify the current row
selectors =. {.@|:@] NB. pick the selector list from the input array
values =. {:@|:      NB. pick the value list from the input array
complements =. [: I. [: -. [: |."1 ] #:~ 2 #~ #@[   NB. Find the complementary columns for each selector
nextvals =. values@>@{.@] * >@{:@] { [ f.    NB. for each complementary column pick the matrix element
                                             NB. and multiply by the minor value
nextsels =. selectors@>@{.@] + >@{:@] (34 b.) 1: f.    NB. for each complementary column calculate the
                                                       NB. corresponding selector
collectminors =. (([ , [: +/ ])/..)&,  NB. flatten the arrays of new values and new selectors,
                                       NB. collect by new minor selectors and sum over the new values
compress =. [: (1 2 $ 0)"_^:(0&=@#) 0&~:@values # ] f. NB. remove zero-valued minors,
                                                       NB. but keep at least one entity
addrow =. getrow ([ ([: compress nextsels collectminors nextvals) ] ; [ complements selectors) ] f. NB. add one row
permpow =. ]`(#@])`[    NB. in power of verb: put the right arg (the matrix) to the left,
                        NB. use the number of rows as the power and put the left arg to the right
perm =. ((<0 1) { (1 2 $ 0 1)"_ addrow^:permpow ]) f.    NB. run N steps and return the value
perm y                  NB. calculate the permanent
)

The next exercise would be avoiding zero valued minors instead of removing them afterwards.

See the full source on GitHub.

Perl

It should be easier to follow the logic in the Perl version. The core code calculating the permanent is:

# An attempt to implement the "Johnson-Gentleman tree minor algorithm":
# This is a non-recursive approach that avoids the re-examination of
# minors appearing in recursive approaches.  It does not split the task
# of calculating a determinant/permanent into smaller tasks but instead
# builds the whole result by extending from single elements.  This takes
# a lot of memory for larger matrices.
#
sub permanent ($a) {
    my $node;
    my $last = $#$a;
    my $sel;

    # Nodes are key-value pairs, where the keys are integers with bits
    # set for the selected columns forming a minor matrix and the
    # corresponding sub-permanent as values.

    # The first node is the empty minor
    $node->{0} = 1;

    # Loop over all rows
    for my $i (0 .. $last) {
        my $next;
        my $row = $a->[$i];
        # Loop over all minors of the previous node.
        for my $minor (keys %$node) {
            # Loop over all columns.  Process only nonzero elements
            # and columns that are not part of the current minor.
            # Add the minor's part to the next larger minor extended
            # by the current column.
            $row->[$_] && !(($sel = 1 << $_) & $minor) &&
            ($next->{$minor | $sel} += $row->[$_] * $node->{$minor})
            for 0 .. $last;
        }
        $node = $next;
    }
    
    # At the end, there is only one value left - the permanent of the
    # whole matrix.
    (values %$node)[0];
}

This version is more efficient as the J counterpart as there is no need to actually remove zero valued minors. They just do not appear.

See the full source on GitHub.

Results

When running these implementations with limits for virtual / real memory as 16 GB / 8 GB on my PC:

• J computes the numbers \(a(0),\ldots,a(44)\) in about 18 min
• Perl computes the numbers \(a(0),\ldots,a(49)\) in about 90 min.

And these are the results:

References

• https://oeis.org/A320843
• https://en.wikipedia.org/wiki/Permanent_(mathematics)
• https://en.wikipedia.org/wiki/Computing_the_permanent
• W. M. Gentleman, S. C. Johnson: The Evaluation of Determinants by
  Expansion by Minors and the General Problem of Substitution, Math
  Comp, Vol. 28, Nr. 126 (1974)
  https://www.ams.org/journals/mcom/1974-28-126/S0025-5718-1974-0373369-7/S0025-5718-1974-0373369-7.pdf