Untying The Gordian Knot via Experimental Mathematics

TL;DR

This paper advocates for an experimental mathematics approach in enumerative combinatorics, demonstrating that computational data and heuristic guessing can lead to rigorous solutions for complex problems.

Contribution

It introduces a methodology where computational experiments guide the discovery and proof of combinatorial results, emphasizing the value of heuristic and experimental techniques.

Findings

Computer-generated data can suggest correct solutions.

Heuristic 'guess and verify' methods can be rigorously justified.

The approach simplifies solving complex enumerative problems.

Abstract

In this methodological article on experimental-yet-rigorous enumerative combinatorics, we use two instructive case studies, to show that often, just like Alexander the Great before us, the simple, "cheating" solution to a hard problem is the best. So before you spend days (and possibly years) trying to answer a mathematical question by analyzing and trying to 'understand' its structure, let your computer generate enough data, and then let it guess the answer. Often its guess can be proved by a quick 'hand-waving' (yet fully rigorous) 'meta-argument'. Since our purpose is to illustrate a methodology, we include many details, as well as Maple source-code.

Equations

$$V=\{v_{ij}|1\le i\le k,1\le j\le n\},$$

$$S=\{\{v_{ij},v_{i^{\prime }{j}^{\prime }}\}||i-i^{\prime }|+|j-j^{\prime }|=1\}.$$

$${\mathcal{B}}_3=\{\{\{1,2,3\}\},\{\{1\},\{2,3\},\{\{1,2\},\{3\}\},\{\{1,3\},\{2\}\},\{\{1\},\{2\},\{3\}\}\}.$$

$${\mathcal{B}}_3=\{123,1/23,12/3,13/2,1/2/3\}.$$

$$A_k(i,j)=|\{A\subseteq E_{n+1}|A\diamond P_j=P_i\}|.$$

$${\mathcal{B}}_2=\{12,1/2\},E_{n+1}=\{\{v_{1,n},v_{1,n+1}\},\{v_{2,n},v_{2,n+1}\},\{v_{1,n+1},v_{2,n+1}\}\}.$$

$$A=\left[\begin{array}{cc}3& 1\\ 2& 1\end{array}\right].$$

$$v_n=\left[\begin{array}{c}{T}_1(n)\\ T_2(n)\end{array}\right],$$

$$v_n=A{v}_{n-1}.$$

$${\chi }_{\lambda }(A)={\lambda }^2-4\lambda +1.$$

$$A^2-4A+1=0.$$

$$T_1(n)=4T_1(n-1)-T_1(n-2),$$

$$\frac{x}{1-4x+x^2}.$$

$$A_3=\left[\begin{array}{ccccc}8& 3& 3& 4& 1\\ 4& 3& 2& 2& 1\\ 4& 2& 3& 2& 1\\ 1& 0& 0& 1& 0\\ 3& 2& 2& 2& 1\end{array}\right].$$

$$L=\left[\begin{array}{cccc}{a}_{12}+\cdots +a_{1n}& -a_{12}& \dots & -a_{1,n}\\ -a_{21}& a_{21}+\cdots +a_{2n}& \dots & -a_{2,n}\\ ⋮& ⋮& \ddots & ⋮\\ -a_{n1}& -a_{n2}& \dots & a_{n1}+\cdots +a_{n,n-1}\end{array}\right].$$

$$\left|\begin{array}{cccc}{a}_{12}+\cdots +a_{1n}& -a_{12}& \dots & -a_{1,n-1}\\ -a_{21}& a_{21}+\cdots +a_{2n}& \dots & -a_{2,n-1}\\ ⋮& ⋮& \ddots & ⋮\\ -a_{n-1,1}& -a_{n-1,2}& \dots & a_{n-1,1}+\cdots +a_{n-1,n}\end{array}\right|.$$

$$L[n]=\sum _{i=1}^{d}R[i]L[n-i].$$

$$\mathtt{CtoR}([[1,1],[1,1]],\mathtt{t})=\frac{1}{-t^2-t+1}.$$

$$F_k(t)=\sum _{n=0}^{\infty }s(k,n)t^n.$$

$$F_1(t)=\frac{t}{1-t}.$$

$$F_2=\frac{t}{t^2-4t+1}.$$

$$F_3=\frac{-t^3+t}{t^4-15t^3+32t^2-15t+1}.$$

$$F_4=\frac{t^7-49t^5+112t^4-49t^3+t}{t^8-56t^7+672t^6-2632t^5+4094t^4-2632t^3+672t^2-56t+1}.$$

$$F_5=\frac{N_5}{D_5}$$

$$N_5=-t^{15}+1440t^{13}-26752t^{12}+185889t^{11}-574750t^{10}+708928t^9-708928t^7$$

$$+574750t^6-185889t^5+26752t^4-1440t^3+t,$$

$$D_5=t^{16}-209t^{15}+11936t^{14}-274208t^{13}+3112032t^{12}-19456019t^{11}+70651107t^{10}-152325888t^9$$

$$+196664896t^8-152325888t^7+70651107t^6-19456019t^5+3112032t^4-274208t^3+11936t^2-209t+1.$$

$$F_6=\frac{N_6}{D_6}$$

$$N_6=t^{31}-33359t^{29}+3642600t^{28}-173371343t^{27}+4540320720t^{26}-70164186331t^{25}$$

Full text

Peter Paule

This article is dedicated to Peter Paule, one of the great pioneers of experimental mathematics and symbolic computation. In particular, it is greatly inspired by his masterpiece, co-authored with Manuel Kauers, ‘The Concrete Tetrahedron’ [KP], where a whole chapter is dedicated to our favorite ansatz, the $C-$finite ansatz.

Introduction

Once upon a time there was a knot that no one could untangle, it was so complicated. Then came Alexander the Great and, in one second, cut it with his sword.

Analogously, many mathematical problems are very hard, and the current party line is that in order for it be considered solved, the solution, or answer, should be given a logical, rigorous, deductive proof.

Suppose that you want to answer the following question:

Find a closed-form formula, as an expression in $n$, for the real part of the $n$-th complex root of the Riemann zeta function, $\zeta(s)$ .

Let’s call this quantity $a(n)$. Then you compute these real numbers, and find out that $a(n)=\frac{1}{2}$ for $n\leq 1000$. Later you are told by Andrew Odlyzko that $a(n)=\frac{1}{2}$ for all $1\leq n\leq 10^{10}$. Can you conclude that $a(n)=\frac{1}{2}$ for all $n$? We would, but, at this time of writing, there is no way to deduce it rigorously, so it remains an open problem. It is very possible that one day it will turn out that $a(n)$ (the real part of the $n$-th complex root of $\zeta(s)$) belongs to a certain ansatz, and that checking it for the first $N_{0}$ cases implies its truth in general, but this remains to be seen.

There are also frameworks, e.g. Pisot sequences (see [ESZ], [Z2]), where the *inductive * approach fails miserably.

On the other hand, in order to (rigorously) prove that $1^{3}+2^{3}+3^{3}+\dots+n^{3}=(n(n+1)/2)^{2}$, for every positive integer $n$, it suffices to check it for the five special cases $0\leq n\leq 4$, since both sides are polynomials of degree $4$, hence the difference is a polynomial of degree $\leq 4$, given by five ‘degrees of freedom’.

This is an example of what is called the ‘$N_{0}$ principle’. In the case of a polynomial identity (like this one), $N_{0}$ is simply the degree plus one.

But our favorite ansatz is the $C$-finite ansatz. A sequence of numbers $\{a(n)\}$ ($0\leq n<\infty$) is $C$-finite if it satisfies a linear recurrence equation with constant coefficients. For example the Fibonacci sequence that satisfies $F(n)-F(n-1)-F(n-2)=0$ for $n\geq 2$.

The $C$-finite ansatz is beautifully described in Chapter 4 of the masterpiece ‘The Concrete Tetrahedron’ ([KP]), by Manuel Kauers and Peter Paule, and discussed at length in [Z1].

Here the ‘$N_{0}$ principle’ also holds (see [Z3]), i.e. by looking at the ‘big picture’ one can determine a priori, a positive integer, often not that large, such that checking that $a(n)=b(n)$ for $1\leq n\leq N_{0}$ implies that $a(n)=b(n)$ for all $n>0$.

A sequence $\{a(n)\}_{n=0}^{\infty}$ is $C$-finite if and only if its (ordinary) generating function $f(t):=\sum_{n=0}^{\infty}a(n)\,t^{n}$ is a rational function of $t$, i.e. $f(t)=P(t)/Q(t)$ for some polynomials $P(t)$ and $Q(t)$. For example, famously, the generating function of the Fibonacci sequence is $t/(1-t-t^{2})$.

Phrased in terms of generating functions, the $C$-finite ansatz is the subject of chapter 4 of yet another masterpiece, Richard Stanley’s ‘Enumerative Combinatorics’ (volume 1) ([S]). There it is shown, using the ‘transfer matrix method’ (that originated in physics), that in many combinatorial situations, where there are finitely many states, one is guaranteed, a priori, that the generating function is rational.

Alas, finding this transfer matrix, at each specific case, is not easy! The human has to first figure out the set of states, and then using human ingenuity, figure out how they interact.

A better way is to automate it. Let the computer do the research, and using ‘symbolic dynamical programming’, the computer, automatically, finds the set of states, and constructs, all by itself (without any human pre-processing) the set of states and the transfer matrix. But this may not be so efficient for two reasons. First, at the very end, one has to invert a matrix with symbolic entries, hence compute symbolic determinants, that is time-consuming. Second, setting up the ‘infra-structure’ and writing a program that would enable the computer to do ‘machine-learning’ can be very daunting.

In this article, we will describe two case studies where, by ‘general nonsense’, we know that the generating functions are rational, and it is easy to bound the degree of the denominator (alias the order of the recurrence satisfied by the sequence). Hence a simple-minded, empirical, approach of computing the first few terms and then ‘fitting’ a recurrence (equivalently rational function) is possible.

The first case-study concerns counting spanning trees in families of grid-graphs, studied by Paul Raff ([R]), and F.J. Faase ([F]). In their research, the human first analyzes the intricate combinatorics, manually sets up the transfer matrix, and only at the end lets a computer-algebra system evaluate the symbolic determinant.

Our key observation, that enabled us to ‘cut the Gordian knot’ is that the terms of the studied sequences are expressible as numerical determinants. Since computing numerical determinants is so fast, it is easy to compute sufficiently many terms, and then fit the data into a rational function. Since we easily have an upper bound for the degree of the denominator of the rational function, everything is rigorous.

The second case-study is computing generating functions for sequences of determinants of ‘almost diagonal Toeplitz matrices’. Here, in addition to the ‘naive’ approach of cranking enough data and then fitting it into a rational function, we also describe the ‘symbolic dynamical programming method’, that surprisingly is faster for the range of examples that we considered. But we believe that for sufficiently large cases, the naive approach will eventually be more efficient, since the ‘deductive’ approach works equally well for the analogous problem of finding the sequence of permanents of these almost diagonal Toeplitz matrices, for which the naive approach will soon be intractable.

This article may be viewed as a tutorial, hence we include lots of implementation details, and Maple code. We hope that it will inspire readers (and their computers!) to apply it in other situations

Accompanying Maple Packages

This article is accompanied by three Maple packages, GFMatrix.txt, JointConductance.txt, and SpanningTrees.txt, all available from the url

http://sites.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/gordian.html

In that page there are also links to numerous sample input and output files.

The human approach to enumerating spanning trees of grid graphs

In order to illustrate the advantage of “keeping it simple”, we will review the human approach to the enumeration task that we will later redo using the ‘Gordian knot’ way. While the human approach is definitely interesting for its own sake, it is rather painful.

Our goal is to enumerate the number of spanning trees in certain families of graphs, notably grid graphs and their generalizations. Let’s examine Paul Raff’s interesting approach described in his paper Spanning Trees in Grid Graph [R]. Raff’s approach was inspired by the pioneering work of F. J. Faase ([F]).

The goal is to find generating functions that enumerate spanning trees in grid graphs and the product of an arbitrary graph and a path or a cycle.

Grid graphs have two parameters, let’s call them $k$ and $n$. For a $k\times n$ grid graph, let’s think of $k$ as fixed while $n$ is the discrete input variable of interest.

Definition The $k\times n$ grid graph $G_{k}(n)$ is the following graph given in terms of its vertex set $V$ and edge set $E$:

[TABLE]

[TABLE]

The main idea in the human approach is to consider the collection of set-partitions of $[k]=\{1,2,\dots,k\}$ and figure out the transition when we extend a $k\times n$ grid graph to a $k\times(n+1)$ one.

Let $\mathcal{B}_{k}$ be the collection of all set-partitions of $[k]$. $B_{k}=|\mathcal{B}_{k}|$ is called the $k$-th Bell number. Famously, the exponential generating function of $B_{k}$, namely $\sum_{k=0}^{\infty}\frac{B_{k}}{k!}\,t^{k}$, equals $e^{e^{t}-1}$.

A lexicographic ordering on $\mathcal{B}_{k}$ is defined as follows:

Definition Given two partitions $P_{1}$ and $P_{2}$ of $[k]$, for $i\in[k]$, let $X_{i}$ be the block of $P_{1}$ containing $i$ and $Y_{i}$ be the block of $P_{2}$ containing $i$. Let $j$ be the minimum number such that $X_{i}\neq Y_{i}$. Then $P_{1}<P_{2}$ iff

1. $|P_{1}|<|P_{2}|$ or

2. $|P_{1}|=|P_{2}|$ and $X_{j}\prec Y_{j}$ where $\prec$ denotes the normal lexicographic order.

For example, here is the ordering for $k=3$:

[TABLE]

For simplicity, we can rewrite it as follows:

[TABLE]

Definition Given a spanning forest $F$ of $G_{k}(n)$, the partition induced by $F$ is obtained from the equivalence relation

$i\sim j\Longleftrightarrow v_{n,i},v_{n,j}$ are in the same component of $F$.

For example, the partition induced by any spanning tree of $G_{k}(n)$ is $123\dots k$ because by definition, in a spanning tree, all $v_{n,i},1\leq i\leq k$ are in the same component. For the other extreme, where every component only consists of one vertex, the corresponding set-partition is $1/2/3/\dots/k-1/k$ because no two $v_{n,i},v_{n,j}$ are in the same component for $1\leq i<j\leq k$.

Definition Given a spanning forest $F$ of $G_{k}(n)$ and a set-partition $P$ of $[k]$, we say that $F$ is consistent with $P$ if:

1. The number of trees in $F$ is precisely $|P|$.

2. $P$ is the partition induced by $F$.

Let $E_{n}$ be the set of edges $E(G_{k}(n))\backslash E(G_{k}(n-1))$, then $E_{n}$ has $2k-1$ members.

Given a forest $F$ of $G_{k}(n-1)$ and some subset $X\subseteq E_{n}$, we can combine them to get a forest of $G_{k}(n)$ as follows. We just need to know how many subsets of $E_{n}$ can transfer a forest consistent with some partition to a forest consistent with another partition. This leads to the following definition:

Definition Given two partitions $P_{1}$ and $P_{2}$ in $\mathcal{B}_{k}$, a subset $X\subseteq E_{n}$ transfers from $P_{1}$ to $P_{2}$ if a forest consistent with $P_{1}$ becomes a forest consistent with $P_{2}$ after the addition of $X$. In this case, we write $X\diamond P_{1}=P_{2}$.

With the above definitions, it is natural to define a $B_{k}\times B_{k}$ transfer matrix $A_{k}$ by the following:

[TABLE]

Let’s look at the $k=2$ case as an example. We have

[TABLE]

For simplicity, let’s call the edges in $E_{n+1}$ $e_{1},e_{2},e_{3}$. Then to transfer the set-partition $P_{1}=12$ to itself, we have the following three ways: $\{e_{1},e_{2}\},\{e_{1},e_{3}\},\{e_{2},e_{3}\}$. In order to transfer the partition $P_{2}=1/2$ into $P_{1}$, we only have one way, namely: $\{e_{1},e_{2},e_{3}\}$. Similarly, there are two ways to transfer $P_{1}$ to $P_{2}$ and one way to transfer $P_{2}$ to itself Hence the transfer matrix is the following $2\times 2$ matrix:

[TABLE]

Let $T_{1}(n),T_{2}(n)$ be the number of forests of $G_{k}(n)$ which are consistent with the partitions $P_{1}$ and $P_{2}$, respectively. Let

[TABLE]

then

[TABLE]

The characteristic polynomial of $A$ is

[TABLE]

By the Cayley-Hamilton Theorem, $A$ satisfies

[TABLE]

Hence the recurrence relation for $T_{1}(n)$ is

[TABLE]

the sequence is $\{1,4,15,56,209,780,2911,10864,40545,151316,\dots\}$ (OEIS A001353) and the generating function is

[TABLE]

Similarly, for the $k=3$ case, the transfer matrix

[TABLE]

The transfer matrix method can be generalized to general graphs of the form $G\times P_{n}$, especially cylinder graphs.

As one can see, we had to think very hard. First we had to establish a ‘canonical’ ordering over set-partitions, then define the consistence between partitions and forests, then look for the transfer matrix and finally worry about initial conditions.

Rather than think so hard, let’s compute sufficiently many terms of the enumeration sequence, and try to guess a linear recurrence equation with constant coefficients, that would be provable a posteriori just because we know that there exists a transfer matrix without worrying about finding it explicitly. But how do we generate sufficiently many terms? Luckily, we can use the celebrated Matrix Tree Theorem.

The Matrix Tree Theorem

Matrix Tree Theorem If $A=(a_{ij})$ is the adjacency matrix of an arbitrary graph $G$, then the number of spanning trees is equal to the determinant of any co-factor of the Laplacian matrix $L$ of $G$, where

[TABLE]

For instance, taking the $(n,n)$ co-factor, we have that the number of spanning trees of $G$ equals

[TABLE]

Since computing determinants for numeric matrices is very fast, we can find the generating functions for the number of spanning trees in grid graphs and more generalized graphs by experimental methods, using the C-finite ansatz.

The GuessRec Maple procedure

Our engine is the Maple procedure GuessRec(L) that resides in the Maple packages accompanying this article. We used the ‘vanilla’, straightforward, linear algebra approach for guessing, using undetermined coefficients. A more efficient way is via the celebrated Berlekamp-Massey algorithm ([Wi]). Since the guessing part is not the bottle-neck of our approach ( it is rather the data-generation part), we preferred to keep it simple.

Naturally, we need to collect enough data. The input is the data (given as a list) and the output is a conjectured recurrence relation derived from that data.

Procedure GuessRec(L) inputs a list, L, and attempts to output a linear recurrence equation with constant coefficients satisfied by the list. It is based on procedure GuessRec1(L,d) that looks for such a recurrence of order $d$.

The output of GuessRec1(L,d) consists of the the list of initial $d$ values (‘initial conditions’) and the recurrence equation represented as a list. For instance, if the input is $L=[1,1,1,1,1,1]$ and $d=1$, then the output will be $[[1],[1]]$; if the input is $L=[1,4,15,56,209,780,2911,10864,40545,151316]$ as the $k=2$ case for grid graphs and $d=2$, then the output will be $[[1,4],[4,-1]]$. This means that our sequence satisfies the recurrence $a(n)=4a(n-1)-a(n-2)$, subject to the initial conditions $a(0)=1,a(1)=4$.

Here is the Maple code:

GuessRec1:=proc(L,d) local eq,var,a,i,n:

if nops(L)<=2*d+2 then

 print(‘The list must be of size >=‘, 2*d+3 ):

 RETURN(FAIL):

fi:

var:=$\{$seq(a[i],i=1..d)$\}$:

eq:=$\{$seq(L[n]-add(a[i]*L[n-i],i=1..d),n=d+1..nops(L))$\}$:

var:=solve(eq,var):

if var=NULL then

 RETURN(FAIL):

else

 RETURN([[op(1..d,L)],[seq(subs(var,a[i]),i=1..d)]]):

fi:

end:

The idea is that having a long enough list $L$ $(|L|>2d+2)$ of data, we use the data after the $d$-th one to discover whether there exists a linear recurrence relation, the first $d$ data points being the initial condition. With the unknowns $a_{1},a_{2},\dots,a_{d}$, we have a linear systems of no less than $d+3$ equations. If there is a solution, it is extremely likely that the recurrence relation holds in general. The first list of length $d$ in the output constitutes the list of initial conditions while the second list, $R$, codes the linear recurrence, where $[R[1],\dots R[d]]$ stands for the following recurrence:

[TABLE]

Here is the Maple procedure GuessRec(L):

GuessRec:=proc(L) local gu,d:

for d from 1 to trunc(nops(L)/2)-2 do

 gu:=GuessRec1(L,d):

 if gu<>FAIL then

  RETURN(gu):

fi:

od:

FAIL:

end:

This procedure inputs a sequence $L$ and tries to guess a recurrence equation with constant coefficients satisfying it. It returns the initial values and the recurrence equation as a pair of lists. Since the length of $L$ is limited, the maximum degree of the recurrence cannot be more than $\lfloor|L|/2-2\rfloor$. With this procedure, we just need to input $L=[1,4,15,56,209,780,2911,10864,40545,151316]$ to get the recurrence (and initial conditions) $[[1,4],[4,-1]]$.

Once the recurrence relation, let’s call it S, is discovered, procedure CtoR(S,t) finds the generating function for the sequence. Here is the Maple code:

CtoR:=proc(S,t) local D1,i,N1,L1,f,f1,L:

if not (type(S,list) and nops(S)=2 and type(S[1],list) and type(S[2],list)

 and nops(S[1])=nops(S[2]) and type(t, symbol) ) then

 print(‘Bad input‘):

 RETURN(FAIL):

fi:

D1:=1-add(S[2][i]*t**i,i=1..nops(S[2])):

N1:=add(S[1][i]t*(i-1),i=1..nops(S[1])):

L1:=expand(D1*N1):

L1:=add(coeff(L1,t,i)*t**i,i=0..nops(S[1])-1):

f:=L1/D1:

L:=degree(D1,t)+10:

f1:=taylor(f,t=0,L+1):

if expand([seq(coeff(f1,t,i),i=0..L)])<>expand(SeqFromRec(S,L+1)) then

print([seq(coeff(f1,t,i),i=0..L)],SeqFromRec(S,L+1)):

 RETURN(FAIL):

else

 RETURN(f):

fi:

end:

Procedure SeqFromRec used above (see the package) simply generates many terms using the recurrence.

Procedure CtoR(S,t) outputs the rational function in $t$, whose coefficients are the members of the C-finite sequence $S$. For example:

[TABLE]

Briefly, the idea is that the denominator of the rational function can be easily determined by the recurrence relation and we use the initial condition to find the starting terms of the generating function, then multiply it by the denominator, yielding the numerator.

Application of GuessRec for enumerating spanning trees of grid graphs and $G\times P_{n}$

With the powerful procedures GuessRec and CtoR, we are able to find generating functions for the number of spanning trees of generalized graphs of the form $G\times P_{n}$. We will illustrate the application of GuessRec to finding the generating function for the number of spanning trees in grid graphs.

First, using procedure GridMN(k,n), we get the $k\times n$ grid graph.

Then, procedure SpFn uses the Matrix Tree Theorem to evaluate the determinant of the co-factor of the Laplacian matrix of the grid graph which is the number of spanning trees in this particular graph. For a fixed $k$, we need to generate a sufficiently long list of data for the number of spanning trees in $G_{k}(n),n\in[l(k),u(k)]$. The lower bound $l(k)$ can’t be too small since the first several terms are the initial condition; the upper bound $u(k)$ can’t be too small as well since we need sufficient data to obtain the recurrence relation. Notice that there is a symmetry for the recurrence relation, and to take advantage of this fact, we modified GuessRec to get the more efficient GuessSymRec (requiring less data). Once the recurrence relation, and the initial conditions, are given, applying CtoR(S,t) will give the desirable generating function, that, of course, is a rational function of $t$. All the above is incorporated in procedure GFGridKN(k,t) which inputs a positive integer $k$ and a symbol $t$, and outputs the generating function whose coefficient of $t^{n}$ is the number of spanning trees in $G_{k}(n)$, i.e. if we let $s(k,n)$ be the number of spanning trees in $G_{k}(n)$, the generating function

[TABLE]

We now list the generating functions $F_{k}(t)$ for $1\leq k\leq 7$: Except for $k=7$, these were already found by Raff[R] and Faase[F], but it is reassuring that, using our new approach, we got the same output. The case $k=7$ seems to be new.

Theorem 1 The generating function for the number of spanning trees in $G_{1}(n)$ is:

[TABLE]

Theorem 2 The generating function for the number of spanning trees in $G_{2}(n)$ is:

[TABLE]

Theorem 3 The generating function for the number of spanning trees in $G_{3}(n)$ is:

[TABLE]

Theorem 4 The generating function for the number of spanning trees in $G_{4}(n)$ is:

[TABLE]

For $5\leq k\leq 7$, since the formulas are too long, we present their numerators and denominators separately.

Theorem 5 The generating function for the number of spanning trees in $G_{5}(n)$ is:

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Theorem 6 The generating function for the number of spanning trees in $G_{6}(n)$ is:

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Theorem 7 The generating function for the number of spanning trees in $G_{7}(n)$ is:

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Note that, surprisingly, the degree of the denominator of $F_{7}(t)$ is $48$ rather than the expected $64$ since the first six generating functions’ denominator have degree $2^{k-1}$, $1\leq k\leq 6$. With a larger computer, one should be able to compute $F_{k}$ for larger $k$, using this experimental approach.

Generally, for an arbitrary graph $G$, we consider the number of spanning trees in $G\times P_{n}$. With the same methodology, a list of data can be obtained empirically with which a generating function follows.

Joint Resistance

The original motivation for the Matrix Tree Theorem, first discovered by Kirchhoff (of Kirchhoff’s laws fame) came from the desire to efficiently compute joint resistances in an electrical network.

Suppose one is interested in the joint resistance in an electric network in the form of a grid graph between two diagonal vertices $[1,1]$ and $[k,n]$. We assume that each edge has resistance $1$ Ohm. To obtain it, all we need is, in addition for the number of spanning trees (that’s the numerator), the number of spanning forests $SF_{k}(n)$ of the graph $G_{k}(n)$ that have exactly two components, each component containing exactly one of the members of the pair $\{[1,1],[k,n]\}$ (this is the denominator). The joint resistance is just the ratio.

In principle, we can apply the same method to obtain the generating function $S_{k}$. Empirically, we found that the denominator of $S_{k}$ is always the square of the denominator of $F_{k}$ times another polynomial $C_{k}$. Once the denominator is known, we can find the numerator in the same way as above. So our focus is to find $C_{k}$.

The procedure DenomSFKN(k,t) in the Maple package JointConductance.txt, calculates $C_{k}$. For $2\leq k\leq 4$, we have

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Remark By looking at the output of our Maple package, we conjectured that $R(k,n)$, the resistance between vertex $[1,1]$ and vertex $[k,n]$ in the $k\times n$ grid graph, $G_{k}(n)$, where each edge is a resistor of $1$ Ohm, is asymptotically $n/k$, for any fixed $k$, as $n\rightarrow\infty$. We proved it rigorously for $k\leq 6$, and we wondered whether there is a human-generated “electric proof”. Naturally we emailed Peter Doyle, the co-author of the delightful masterpiece [DS], who quickly came up with the following argument.

Making the horizontal resistors into almost resistance-less gold wires gives the lower bound $R(k,n)\geq(n-1)/k$ since it is a parallel circuit of $k$ resistors of $n-1$ Ohms. For an upper bound of the same order, put 1 Ampere in at [1,1] and out at $[k,n]$, routing $1/k$ Ampere up each of the $k$ verticals. The energy dissipation is $k(n-1)/k^{2}+C(k)=(n-1)/k+C(k)$, where the constant $C(k)$ is the energy dissipated along the top and bottom resistors. Specifically, $C(k)=2(1-1/k)^{2}+(1-2/k)^{2}+\dots+(1/k)^{2})$. So $(n-1)/k\leq R(k,n)\leq(n-1)/k+C(k)$.

We thank Peter Doyle for his kind permission to reproduce this electrifying argument.

The statistic of the number of vertical edges in spanning trees of grid graphs

Often in enumerative combinatorics, the class of interest has natural ‘statistics’, like height, weight, and IQ for humans. Recall that the naive counting is

[TABLE]

getting a number. Define:

[TABLE]

where $f:=A\rightarrow\mathbb{Z}$ is the statistic of interest. To go from the weighted enumeration (a certain Laurent polynomial) to straight enumeration, one simply plugs-in $x=1$, i.e. $|A|_{1}=|A|$.

The scaled random variable is defined as follows. Let $E(f)$ and $Var(f)$ be the expectation and variance, respectively, of the statistic $f$ defined on $A$, and define the scaled random variable, for $a\in A$, by

[TABLE]

In this section, we are interested in the statistic ‘number of vertical edges’, defined on spanning trees of grid graphs. For given $k$ and $n$, let, as above, $G_{k}(n)$ denote the $k\times n$ grid-graph. Let $\mathcal{F}_{k,n}$ be its set of spanning trees. If the weight is 1, then $\sum_{f\in\mathcal{F}_{k,n}}1=|\mathcal{F}_{k,n}|$ is the naive counting. Now let’s define a natural statistic

$ver(T)$ = the number of vertical edges in the spanning tree $T$,

and the weight $w(T)=v^{ver(T)}$, then the weighted counting follows:

[TABLE]

where $\mathcal{F}_{k,n}$ is the set of spanning trees of $G_{k}(n)$.

We define the bivariate generating function

[TABLE]

More generally, with our Maple package GFMatrix.txt, and procedure VerGF, we are able to obtain the bivariate generating function for an arbitrary graph of the form $G\times P_{n}$. The procedure VerGF takes inputs $G$ (an arbitrary graph), $N$ (an integer determining how many data we use to find the recurrence relation) and two symbols $v$ and $t$.

The main tool for computing VerGF is still the Matrix Tree Theorem and GuessRec. But we need to modify the Laplacian matrix for the graph. Instead of letting $a_{ij}=-1$ for $i\neq j$ and $\{i,j\}\in E(G\times P_{n})$, we should consider whether the edge $\{i,j\}$ is a vertical edge. If so, we let $a_{i,j}=-v,a_{j,i}=-v$. The diagonal elements which are $(-1)\times$ (the sum of the rest entries on the same row) should change accordingly. The following theorems are for grid graphs when $2\leq k\leq 4$ while $k=1$ is a trivial case because there are no vertical edges.

Theorem 8 The bivariate generating function for the weighted counting according to the number of vertical edges of spanning trees in $G_{2}(n)$ is:

[TABLE]

Theorem 9 The bivariate generating function for the weighted counting according to the number of vertical edges vertical edges of spanning trees in $G_{3}(n)$ is:

[TABLE]

Theorem 10 The bivariate generating function for the weighted counting according to the number of vertical edges of spanning trees in $G_{4}(n)$ is:

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

With the Maple package BiVariateMoms.txt and its Story procedure from

http://sites.math.rutgers.edu/~zeilberg/tokhniot/BiVariateMoms.txt,

the expectation, variance and higher moments can be easily analyzed. We calculated up to the 4th moment for $G_{2}(n)$. For $k=3,4$, you can find the output files from

http://sites.math.rutgers.edu/~yao/OutputStatisticVerticalk=3.txt

http://sites.math.rutgers.edu/~yao/OutputStatisticVerticalk=4.txt

Theorem 11 The moments of the statistic: the number of vertical edges in the spanning trees of $G_{2}(n)$ are as follows:

Let $b$ be the largest positive root of the polynomial equation

[TABLE]

whose floating-point approximation is 3.732050808, then the size of the $n$-th family (i.e. straight enumeration) is very close to

[TABLE]

The average of the statistics is, asymptotically

[TABLE]

The variance of the statistics is, asymptotically

[TABLE]

The skewness of the statistics is, asymptotically

[TABLE]

The kurtosis of the statistics is, asymptotically

[TABLE]

Application of the C-finite Ansatz to computing generating functions of determinants (and permanents) of almost-diagonal Toeplitz matrices

So far, we have seen applications of the $C$-finite ansatz methodology for automatically computing generating functions for enumerating spanning trees/forests for certain infinite families of graphs.

The second case study is completely different, and in a sense more general, since the former framework may be subsumed in this new context.

Definition Diagonal matrices $A$ are square matrices in which the entries outside the main diagonal are [math], i.e. $a_{ij}=0$ if $i\neq j$.

Definition An almost-diagonal Toeplitz matrix $A$ is a square matrices in which $a_{i,j}=0$ if $j-i\geq k_{1}$ or $i-j\geq k_{2}$ for some fixed positive integers $k_{1},k_{2}$ and $\forall i_{1},j_{1},i_{2},j_{2}$, if $i_{1}-j_{1}=i_{2}-j_{2}$, then $a_{i_{1}j_{1}}=a_{i_{2}j_{2}}$.

For simplicity, we use the notation $L=$[n, [the first $k_{1}$ entries in the first row], [the first $k_{2}$ entries in the first column]] to denote the $n\times n$ matrix with these specifications. Note that this notation already contains all information we need to reconstruct this matrix. For example, [6, [1,2,3], [1,4]] is the matrix

[TABLE]

The following is the Maple procedure DiagMatrixL (in our Maple package GFMatrix.txt), which inputs such a list $L$ and outputs the corresponding matrix.

DiagMatrixL:=proc(L) local n, r1, c1,p,q,S,M,i:

n:=L[1]:

r1:=L[2]:

c1:=L[3]:

p:=nops(r1)-1:

q:=nops(c1)-1:

if r1[1] <> c1[1] then

 return fail:

fi:

S:=[0$(n-1-q), seq(c1[q-i+1],i=0..q-1), op(r1), 0$(n-1-p)]:

M:=[0$n]:

for i from 1 to n do

 M[i]:=[op(max(0,n-1-q)+q+2-i..max(0,n-1-q)+q+1+n-i,S)]:

od:

return M:

end:

For this matrix, $k_{1}=3$ and $k_{2}=2$. Let $k_{1},k_{2}$ be fixed and $M_{1},M_{2}$ be two lists of numbers or symbols of length $k_{1}$ and $k_{2}$ respectively, $A_{k}$ is the almost-diagonal Toeplitz matrix represented by the list $L_{k}=[k,M_{1},M_{2}]$. Note that the first elements in the lists $M_{1}$ and $M_{2}$ must be identical.

Having fixed two lists $M_{1}$ of length $k_{1}$ and $M_{2}$ of length $k_{2}$, (where $M_{1}[1]=M_{2}[1]$), it is of interest to derive automatically, the generating function (that is always a rational function for reasons that will soon become clear), $\sum_{k=0}^{\infty}a_{k}\,t^{k}$, where $a_{k}$ denotes the determinant of the $k\times k$ almost-diagonal Toeplitz matrix whose first row starts with $M_{1}$, and first column starts with $M_{2}$. Analogously, it is also of interest to do the analogous problem when the determinant is replaced by the permanent.

Here is the Maple procedure GFfamilyDet which takes inputs (i) $A$: a name of a Maple procedure that inputs an integer $n$ and outputs an $n\times n$ matrix according to some rule, e.g., the almost-diagonal Toeplitz matrices, (ii) a variable name $t$, (iii) two integers $m$ and $n$ which are the lower and upper bounds of the sequence of determinants we consider. It outputs a rational function in $t$, say $R(t)$, which is the generating function of the sequence.

GFfamilyDet:=proc(A,t,m,n) local i,rec,GF,B,gu,Denom,L,Numer:

L:=[seq(det(A(i)),i=1..n)]:

rec:=GuessRec([op(m..n,L)])[2]:

gu:=solve(B-1-add(t**i*rec[i]*B,i=1..nops(rec)), B):

Denom:=denom(subs(gu,B)):

Numer:=Denom*(1+add(L[i]*t**i, i=1..n)):

Numer:=add(coeff(Numer,t,i)*t**i, i=0..degree(Denom,t)):

Numer/Denom:

end:

Similarly we have procedure GFfamilyPer for the permanent. Let’s look at an example. The following is a sample procedure which considers the family of almost diagonal Toeplitz matrices which the first row $[2,3]$ and the first column $[2,4,5]$.

SampleB:=proc(n) local L,M:

L:=[n, [2,3], [2,4,5]]:

M:=DiagMatrixL(L):

end:

Then GFfamilyDet(SampleB, t, 10, 50) will return the generating function

[TABLE]

It turns out, that for this problem, the more ‘conceptual’ approach of setting up a transfer matrix also works well. But don’t worry, the computer can do the ‘research’ all by itself, with only a minimum amount of human pre-processing.

We will now describe this more conceptual approach, that may be called symbolic dynamical programming, where the computer sets up, automatically, a finite-state scheme, by dynamically discovering the set of states, and automatically figures out the transfer matrix.

The Transfer Matrix method for almost-diagonal Toeplitz matrices

Recall from Linear Algebra 101, the

Cofactor Expansion Let $|A|$ denote the determinant of an $n\times n$ matrix $A$, then

[TABLE]

where $M_{ij}$ is the $(i,j)$-minor.

We’d like to consider the Cofactor Expansion for almost-diagonal Toeplitz matrices along the first row. For simplicity, we assume while $a_{i,j}=0$ if $j-i\geq k_{1}$ or $i-j\geq k_{2}$ for some fixed positive integers $k_{1},k_{2}$, and if $-k_{2}<j_{1}-i_{1}<j_{2}-i_{2}<k_{1}$, then $a_{i_{1}j_{1}}\neq a_{i_{2}j_{2}}$. Under this assumption, for any minors we obtain through recursive Cofactor Expansion along the first row, the dimension, the first row and the first column should provide enough information to reconstruct the matrix.

For an almost-diagonal Toeplitz matrix represented by $L=$[Dimension, [the first $k_{1}$ entries in the first row], [the first $k_{2}$ entries in the first column]], any minor can be represented by [Dimension, [entries in the first row up to the last nonzero entry], [entries in the first column up to the last nonzero entry]].

Our goal in this section is the same as the last one, to get a generating function for the determinant or permanent of almost-diagonal Toeplitz matrices $A_{k}$ with dimension $k$. Once we have those almost-diagonal Toeplitz matrices, the first step is to do a one-step expansion as follows:

ExpandMatrixL:=proc(L,L1)

local n,R,C,dim,R1,C1,i,r,S,candidate,newrow,newcol,gu,mu,temp,p,q,j:

n:=L[1]:

R:=L[2]:

C:=L[3]:

p:=nops(R)-1:

q:=nops(C)-1:

dim:=L1[1]:

R1:=L1[2]:

C1:=L1[3]:

if R1=[] or C1=[] then

 return :

elif R[1]<>C[1] or R1[1]<>C1[1] or dim>n then

 return fail:

else

S:={}:

gu:=[0$(n-1-q), seq(C[q-i+1],i=0..q-1), op(R), 0$(n-1-p)]:

candidate:=[0$nops(R1),R1[-1]]:

for i from 1 to nops(R1) do

 mu:=R1[i]:

for j from n-q to nops(gu) do

  if gu[j]=mu then

   candidate[i]:=gu[j-1]:

  fi:

 od:

od:

for i from n-q to nops(gu) do

 if gu[i] = R1[2] then

  temp:=i:

  break:

 fi:

od:

for i from 1 to nops(R1) do

 if i = 1 then

  mu:=[R1[i]*(-1)**(i+1), [dim-1,[op(i+1..nops(candidate), candidate)], [seq(gu[temp-i],i=1..temp-n+q)]]]:

  S:=S union mu:

 else

  mu:=[R1[i]*(-1)**(i+1), [dim-1, [op(1..i-1, candidate), op(i+1..nops(candidate), candidate)], [op(2..nops(C1), C1)]]]:

  S:=S union mu:

 fi:

od:

 return S:

fi:

end:

The ExpandMatrixL procedure inputs a data structure $L=$ [Dimension, first_row=[ ], first_col=[ ]] as the matrix we start and the other data structure $L1$ as the current minor we have, expands $L1$ along its first row and outputs a list of [multiplicity, data structure].

We would like to generate all the ”children” of an almost-diagonal Toeplitz matrix regardless of the dimension, i.e., two lists $L$ represent the same child as long as their first_rows and first_columns are the same, respectively. The set of ”children” is the scheme of the almost diagonal Toeplitz matrices in this case.

The following is the Maple procedure ChildrenMatrixL which inputs a data structure $L$ and outputs the set of its ”children” under Cofactor Expansion along the first row:

ChildrenMatrixL:=proc(L) local S,t,T,dim,U,u,s:

dim:=L[1]:

S:={[op(2..3,L)]}:

T:={seq([op(2..3,t[2])],t in ExpandMatrixL(L,L))}:

while T minus S <> {} do

 U:=T minus S:

 S:=S union T:

 T:={}:

 for u in U do

 T:=T union {seq([op(2..3,t[2])],t in ExpandMatrixL(L,[dim,op(u)]))}:

od:

od:

for s in S do

 if s[1]=[] or s[2]=[] then

  S:=S minus {s}:

 fi:

od:

S:

end:

After we have the scheme $S$, by the Cofactor Expansion of any element in the scheme, a system of algebraic equations follows. For children in $S$, it’s convenient to let the almost-diagonal Toeplitz matrix be the first one $C_{1}$ and for the rest, any arbitrary ordering will do. For example, if after Cofactor Expansion for $C_{1}$, $c_{2}$ ”copies” of $C_{2}$ and $c_{3}$ ”copies” of $C_{3}$ are obtained, then the equation will be

[TABLE]

However, if the above equation is for $C_{i},i\neq 1$, i.e. $C_{i}$ is not the almost-diagonal Toeplitz matrix itself, then the equation will be slightly different:

[TABLE]

Here $t$ is a symbol as we assume the generating function is a rational function of $t$.

Here is the Maple code that implements how we get the generating function for the determinant of a family of almost-diagonal Toeplitz matrices by solving a system of algebraic equations:

GFMatrixL:=proc(L,t) local S,dim,var,eq,n,A,i,result,gu,mu:

dim:=L[1]:

S:=ChildrenMatrixL(L):

S:=[[op(2..3,L)], op(S minus {[op(2..3,L)]})]:

n:=nops(S):

var:={seq(A[i],i=1..n)}:

eq:={}:

for i from 1 to 1 do

 result:=ExpandMatrixL(L,[dim,op(S[i])]):

 for gu in result do

  if gu[2][2]=[] or gu[2][3]=[] then

   result:=result minus {gu}:

  fi:

 od:

 eq:=eq union {A[i] - 1 - add(gu[1]tA[CountRank(S, [op(2..3, gu[2])])], gu in result)}:

od:

for i from 2 to n do

 result:=ExpandMatrixL(L,[dim,op(S[i])]):

 for gu in result do

 if gu[2][2]=[] or gu[2][3]=[] then

  result:=result minus gu:

 fi:

od:

eq:=eq union {A[i] - add(gu[1]tA[CountRank(S, [op(2..3, gu[2])])], gu in result)}:

od:

gu:=solve(eq, var)[1]:

subs(gu, A[1]):

end:

GFMatrixL([20, [2, 3], [2, 4, 5]], t) returns

[TABLE]

Compared to empirical approach, the ‘symbolic dynamical programming’ method is faster and more efficient for the moderate-size examples that we tried out. However, as the lists will grow larger, it is likely that the former method will win out, since with this non-guessing approach, it is equally fast to get generating functions for determinants and permanents, and as we all know, permanents are hard.

The advantage of the present method is that it is more appealing to humans, and does not require any ‘meta-level’ act of faith. However, both methods are very versatile and are great experimental approaches for enumerative combinatorics problems. We hope that our readers will find other applications.

Summary

Rather than trying to tackle each enumeration problem, one at a time, using ad hoc human ingenuity each time, building up an intricate transfer matrix, and only using the computer at the end as a symbolic calculator, it is a much better use of our beloved silicon servants (soon to become our masters!) to replace ‘thinking’ by ‘meta-thinking’, i.e. develop experimental mathematics methods that can handle many different types of problems. In the two case studies discussed here, everything was made rigorous, but if one can make semi-rigorous and even non-rigorous discoveries, as long as they are interesting, one should not be hung up on rigorous proofs. In other words, if you can find a rigorous justification (like in these two case studies) that’s nice, but if you can’t, that’s also nice!

Acknowledgment

Many thanks are due to a very careful referee that pointed out many minor, but annoying, errors. Also thanks to Peter Doyle for permission to include his elegant electric argument.