# Rook theoretic proofs of some identities related to Spivey’s Bell number
formula

R.B. Corcino<EMAIL_ADDRESS>R.O. Celeste<EMAIL_ADDRESS>K.J.M.
Gonzales111Author is supported by Department of Science and Technology-Science
Education Institute under its ASTHRD Program<EMAIL_ADDRESS>Mathematics and ICT Department, Cebu Normal University, 6000 Cebu City,
Philippines  Institute of Mathematics, University of the Philippines Diliman,
1101 Quezon City, Philippines

###### Abstract

We use rook placements to prove Spivey’s Bell number formula and other
identities related to it, in particular, some convolution identities involving
Stirling numbers and relations involving Bell numbers. To cover as many
special cases as possible, we work on the generalized Stirling numbers that
arise from the rook model of Goldman and Haglund. An alternative combinatorial
interpretation for the Type II generalized $q$-Stirling numbers of Remmel and
Wachs is also introduced in which the method used to obtain the earlier
identities can be adapted easily.

###### keywords:

Bell numbers , Stirling numbers , rook placement

##  1 Introduction

Let $B(n)$ denote the $n$-th Bell number and $S(n,k)$ denote the Stirling
number of the second kind. Spivey [16] established combinatorially the
following identity for Bell numbers

$B(n+m)=\sum_{k=0}^{n}\sum_{j=0}^{m}j^{n-k}\binom{n}{k}S(m,j)B(k)\,.$ (1)

Various alternative proofs and extensions of this identity have appeared in
the literature. For instance, Gould and Quaintance [6] provided a generating
function proof, which, in turn, was extended by Xu [18] in the case of Hsu and
Shuie’s [7] generalized Stirling numbers. Katriel [8] proved a $q$-analogue
using certain $q$-differential operators. Still, Belbachir and Mihoubi [2]
proved 1 using a decomposition of the Bell polynomial using a certain
polynomial basis. See also [9, Theorem 10], [10, Theorem 5.3] and [13] for
other generalizations and methods. The present authors also derived a
generalization [4, Theorem 4.5] of Spivey’s identity in the case of the
generalized $q$-Stirling numbers that arise from normal ordering.

In this short paper, we derive variations of 1 and other related identities
using the aforementioned rook model which we describe in Section 2. The
results are presented in Section 3. In Section 4, we introduce a new rook
model for the Type II generalized $q$-Stirling numbers of Remmel and Wachs
[15] which is a variation of Goldman and Haglund’s. We show how the method
used in deriving the earlier identities can be easily adapted under this
model.

##  2 The rook model

In this section, we introduce a rook model based on Goldman and Haglund’s [5].
We will later show that these two models are essentially identical and give
rise to the same rook numbers.

Let $w$ be a word consisting of the letters $U$ and $V$. If we let $U$
correspond to a unit horizontal step and $V$ a unit vertical step, then $w$
outlines a Ferrers boards (or simply, board), which we denote by $B(w)$. For
example, Figure 1 shows the board outlined by $UVUUUVVU$. Note that we allow
the rightmost columns and the bottommost rows to be empty, i.e., not
containing any cell. The initial $U$’s and final $V$’s, if any, are extraneous
when we consider rook placements, but they are natural in the context of
normal ordering (see [17]). Alternatively, we can also describe a board by the
length of its columns.

Given a board $B$, we associate to each cell of $B$ an integer called its
_pre-weight_. Fix $s\in\mathbb{R}$. A placement of $k$ rooks on a board $B$ is
a marking of $k$ cells of $B$ with “$\bullet$” such that at most one rook is
placed in each column. We assume that the rooks are placed from right to left
among $k$ chosen column. Once a rook is placed on a cell, $s-1$ is added to
the pre-weight of every cell to its left in the same row. However, if a cell
lies above a rook (in which case it said to be _cancelled_ by the rook), then
the cell is assigned the pre-weight 0. We denote by $C_{k}(B;s)$ the set of
all placement of $k$ rooks on $B$ under the rules we just described.

For the purpose of distinguishing the pre-weights of a board before and after
rooks are placed, we will call the former as _default pre-weights_. We number
the rows and columns of a board $B$ from top to bottom and from right to left,
respectively. The term pre-weight is chosen because it determines the weight
of a board. Specifically, a cell with pre-weight $t$ is assigned the weight
$q^{t}$ if it does not contain a rook, and $[t]_{q}=\frac{q^{t}-1}{q-1}$ if it
contains a rook. The weight of a rook placement $\phi$, denoted by
$\mbox{wt}(\phi)$, is defined as the product of the weight of the cells.

Denote by $J_{n}$ the board outlined by $(VU)^{n}$ whose cells have default
pre-weight 1. Note that columns of $J_{n}$ have lengths $0,1,\ldots n-1$.
Figure 3 shows a rook placement on $J_{5}$ where the values in the cells
indicate the pre-weights. This particular rook placement has weight
$q^{2s+2}[s]_{q}$.

We denote by $S_{s,q}[n,k]$ the sum of the weights of all placements of $n-k$
rooks on $J_{n}$, i.e.,

$S_{s,q}[n,k]=\sum_{\phi\in C_{n-k}(J_{n};s)}\mbox{wt}(\phi)\,.$

######  Remark 1.

If, for instance, a column has two cells where the bottom cell has pre-weight
$c_{1}$ and the top cell has pre-weight $c_{2}$, then the placement of a rook
on the bottom cell has weight $[c_{1}]_{q}$ while the placement of a rook on
the top cell has weight $q^{c_{1}}[c_{2}]_{q}$. Hence, the total weight of all
rook placements on this column is
$[c_{1}]_{q}+q^{c_{1}}[c_{2}]_{q}=[c_{1}+c_{2}]_{q}$. In general, if $c$ is
the sum of the pre-weights of the cells in a column, then the total weight of
all rook placements on the column is $[c]_{q}$.

######  Proposition 1.

The number $S_{s,q}[n,k]$ satisfies the recurrence

$S_{s,q}[n,k]=q^{s(n-1)-(s-1)(k-1)}S_{s,q}[n-1,k-1]+[s(n-1)-(s-1)k]_{q}S_{s,q}[n-1,k]\,,$

with initial conditions $S_{s,q}[n,0]=S_{s,q}[0,n]=\delta_{0,n}$, where
$\delta_{0,n}=1$ if $n=0$ and 0 if $n\neq 0$.

###### Proof.

If there is a rook on the $n$-th column, then the other columns form a rook
placement from $C_{n-1-k}(J_{n-1};s)$. Due to the placement of $n-1-k$ rooks
on the first $n-1$ columns, the $n$-th column has pre-weight
$(n-1)+(s-1)(n-1-k)=s(n-1)-(s-1)k$. Hence, the total weight contributed by all
rook placements on the $n$-th column is $[s(n-1)-(s-1)k]_{q}$.

If there is no rook on the $n$-th column, then the remaining columns form a
rook placement from $C_{n-k}(J_{n-1};s)$. Because of the placement of $n-k$
rooks on the first $n-1$ columns, the $n$-th column has pre-weight
$(n-1)+(s-1)(n-k)=s(n-1)-(s-1)(k-1)$, which implies that the $n$-th column
contributes a weight of $q^{s(n-1)-(s-1)(k-1)}$. ∎

By comparing recurrences, we see that $h^{n-k}S_{s,q}[n,k]$ equals the number
$\mathfrak{S}_{s;h}(n,k|q)$ in [11]. These numbers are coefficients of the
string $(VU)^{n}$ in its normally ordered form, i.e., an equivalent expression
where the $V$’s are to the right of $U$’s obtained using the relation
$UV=qVU+hV^{s}$. When $h=1,q=1$, the cases $s=1$ and $s=0$ produce the usual
Stirling number of the first kind and Stirling number of the second kind,
respectively. (See [11] for other special cases and their relationship with
normal ordering.)

######  Remark 2.

The original Goldman-Haglund rook model is as follows (see [5, Section 7]).
Let $B$ be a board and consider a rook placement $\phi$ on $B$ such that no
two rooks occupy the same column. Given a cell $\gamma$ in $B$, let
$v(\gamma)$ be the number of rooks strictly to the right of, and in same row
as $\gamma$. Define the weight of $\gamma$ to be

$\displaystyle\mbox{wt}(\gamma)=\begin{cases}1&\mbox{if there is a rook
below}\\\ &\hskip 21.68121pt\mbox{and in the same column as $\gamma$,}\\\
[(s-1)v(\gamma)+1]_{q}&\mbox{if $\gamma$ contains a rook,}\\\
q^{(s-1)v(\gamma)+1}&\mbox{else}\,.\end{cases}$

The weight of $\phi$, denoted by $\mbox{wt}(\phi)$ is then defined as the
product of the weight of the cells. (Note that in [5], the boards are bottom
and right justified, the rooks are placed from left to right and a rook
cancels cells that lie below it. In addition, $\alpha$ is used instead of $s$.
The definition of $\mbox{wt}(\gamma)$ above has been modified to reflect the
definitions in the current paper.)

It can be easily seen that the weight of a rook placement under the Goldman-
Haglund model and under the current model are identical. First, note that
$\mbox{wt}(\gamma)$ is precisely the pre-weight of $\gamma$. Under the latter,
a cell lying above a rook has pre-weight 0 and hence, has weight $q^{0}=1$. A
cell containing a rook and has $t$ rooks lying to its right and on the same
row will have pre-weight $t(s-1)+1$ and hence, the weight $[t(s-1)+1]_{q}$.
Finally, a cell not containing a rook and has $t$ rooks lying to its right and
on the same row will have pre-weight $t(s-1)+1$ and thus, the weight
$q^{t(s-1)}+1$. Hence, given a rook placement $\phi$, $\mbox{wt}(\phi)$ is
identical under both models.

The concept of pre-weights provides us with some degree of flexibility in
further generalizing the earlier model. For instance, instead of adding $s-1$
to the pre-weight of each cell to the left of a rook, we can fix a _rule_ $R$
which specifies, given a placement of rook in a cell, the cell (in each column
to the right of a rook) that receives the pre-weight increment. Denote by
$C^{R}_{k}(B;s)$ the set of placement of $k$ rooks on the board $B$ under this
rule. In particular, we denote by $\mathcal{R}$ the rule which specifies that
the $i$-th rook (from the right) adds a pre-weight of $s-1$ to the $i$-th cell
above the bottom cell in every column to its left. Of course, when the default
rule is used, we simply use the notation $C_{k}(B;s)$.

In addition to specifying a rule $R$, we can also assign other values for the
default pre-weights. This will be useful when we consider Type II
$q$-analogues of Stirling numbers in Section 4. For $\alpha\in\mathbb{R}$,
denote by $J^{\prime}_{n,\alpha}$ the board outlined by $(VU)^{n}V$ where the
bottom cell in each column each has default pre-weight $\alpha$. Figure 2
shows one rook placement in $C^{\mathcal{R}}_{3}(J^{\prime}_{4,\alpha};s)$,
which has weight $q^{2s+2\alpha}[\alpha]_{q}^{2}[s]_{q}$.

######  Proposition 2.

Given any rule $R$ and two boards $B_{1}$ and $B_{2}$ with the same number of
non-empty columns such that the corresponding columns have identical total
column pre-weights, we have $\sum_{\phi\in
C^{R}_{k}(B_{1};s)}\mbox{wt}(\phi)=\sum_{\phi\in
C^{R}_{k}(B_{2};s)}\mbox{wt}(\phi)$. Furthermore, given a board $B$ and two
rules $R_{1}$ and $R_{2}$, $\sum_{\phi\in
C^{R_{1}}_{k}(B;s)}\mbox{wt}(\phi)=\sum_{\phi\in
C^{R_{2}}_{k}(B;s)}\mbox{wt}(\phi)$.

###### Proof.

We can establish, using the same reasoning as Remark 1, that the sum of the
weights of all possible rook placements on a particular column is completely
determined by the column’s total pre-weight, regardless of the distribution of
the pre-weights of the individual cells. This proves the first statement. The
second statement follows from the first statement and the fact that different
rules applied on the same board result to boards with identical total column
pre-weights. ∎

##  3 Results

Define the $n$-th generalized Bell polynomial by

$B_{s,q}[n;x]=\sum_{k=0}^{n}S_{s,q}[n,k]x^{k}\,.$

The $n$-th generalized Bell number, denoted by $B_{s,q}[n]$, is defined as
$B_{s,q}[n;1]$. The $q$-binomial coefficients are given by
$\genfrac{[}{]}{0.0pt}{}{n}{k}_{q}=\frac{[n]_{q}!}{[k]_{q}![n-k]_{q}!}$, where
$[n]_{q}!=[n]_{q}[n-1]_{q}\cdots[2]_{q}[1]_{q}$. Clearly,
$\genfrac{[}{]}{0.0pt}{}{n}{k}_{q}=\genfrac{[}{]}{0.0pt}{}{n}{n-k}_{q}$. They
also satisfy the relations (see [12, Table 1 and Identity (2.2)])

$\displaystyle\genfrac{[}{]}{0.0pt}{}{n}{k}_{q}$ $\displaystyle=\sum_{0\leq
t_{1}\leq t_{2}\leq\ldots\leq t_{n-k}\leq k}q^{t_{1}+t_{2}+\cdots+t_{n-k}}$
$\displaystyle\genfrac{[}{]}{0.0pt}{}{n}{k}_{q}$
$\displaystyle=\sum_{t_{0}+t_{1}+t_{2}+\cdots+t_{k}=n-k}q^{0t_{0}+1t_{1}+2t_{2}+\cdots+kt_{k}}\,.$

In this section, we show how Spivey’s identity, specifically, its
generalization involving $S_{s,q}[n,k]$, is a consequence of a convolution
identity. We then proceed by deriving other convolution identities and the
identities for Bell numbers that arise from them.

A version of Lemma 3 and Theorem 5 already appeared in [4], but only the case
where $s\in\mathbb{N}$ and with pre-weights interpreted as subdivisions in
cells was proved in detail.

######  Lemma 3.

Let $\phi$ be a placement of $k$ rooks on the board $J_{n}$. Then, there
exists a unique (possibly empty) collection $\mathcal{C}$ of columns in
$\phi$, such that if $|\mathcal{C}|=\mu+1$, then (a) these columns have a rook
in the bottom $1,2,\ldots,\mu+1$ cells and (b) every column not in
$\mathcal{C}$ contains at least $1+t$ uncanceled cells not containing a rook,
where $t$ is the number of columns in $\mathcal{C}$ to the right of that
column.

###### Proof.

Let $\phi\in C^{R}_{k}(J^{\prime}_{n,\alpha};s)$. The elements of
$\mathcal{C}$ can be obtained iteratively as follows. Let $c_{1}$ be the first
column from the right containing a rook on the bottom cell. If there exists no
such $c_{1}$, then $\mathcal{C}=\varnothing$. Let $c_{2}$ be the next column
containing a rook in one of the bottom two cells, i.e., either on the bottom
cell or on the cell above the bottom cell. If there exists no such $c_{2}$,
then $\mathcal{C}=\\{c_{1}\\}$. We can continue this process as long as needed
until the elements of $\mathcal{C}$ are all determined. Note that this process
also shows the uniqueness of $\mathcal{C}$. ∎

######  Lemma 4.

Let $\alpha\in\mathbb{R},n,k\in\mathbb{N}$. Then, for any rule $R$,

$\displaystyle\sum_{\phi\in
C^{R}_{n-k}(J^{\prime}_{n;\alpha};s)}\mbox{wt}(\phi)$
$\displaystyle=\sum_{r=0}^{n}q^{\alpha
r}\genfrac{[}{]}{0.0pt}{}{n}{r}_{q^{s}}S_{s,q}[r,k]\prod_{i=0}^{n-r-1}[\alpha+si]_{q}$
(2) $\displaystyle\sum_{\phi\in
C^{R}_{n-k}(J^{\prime}_{n;\alpha};s)}\mbox{wt}(\phi)$
$\displaystyle=\sum_{r=0}^{n}q^{\alpha
k}S_{s,q}[n,r]\genfrac{[}{]}{0.0pt}{}{r}{k}_{q^{s-1}}\prod_{i=0}^{r-k-1}[\alpha+i(s-1)]_{q}\,.$
(3)

###### Proof.

We prove 2 first. Let $\phi\in
C^{\mathcal{R}}_{n-k}(J^{\prime}_{n;\alpha};s)$. By Lemma 3, we can write
$\phi$ uniquely as a pair $(\mathcal{C}_{\phi},\phi-\mathcal{C}_{\phi})$,
where $\mathcal{C}_{\phi}$ are the cells of $\phi$ that either lie in the
columns that satisfy the properties in the lemma with $\mu=n-r-1$ or on the
bottom $1+t$ cells of the other columns, and $\phi-\mathcal{C}_{\phi}$ are the
cells in $\phi$ not in $\mathcal{C}_{\phi}$. One sees that
$\phi-\mathcal{C}_{\phi}$ forms some rook placement in $C_{r-k}(J_{r};s)$.
This accounts for the factor $S_{s,q}[r,k]$.

Let $L_{\mu}$ be the set of all such possible $\mathcal{C}_{\phi}$, with
$\phi\in C^{\mathcal{R}}_{n-k}(J^{\prime}_{n;\alpha};s)$ with $\mu=n-r-1$ as
in Lemma 3. It suffices to compute the sum of the weights of the elements of
$L_{\mu}$. Clearly, the columns containing rooks contribute a weight of
$\prod_{i=0}^{n-r-1}[\alpha+si]_{q}$ and the bottom $r$ cells of the columns
not containing rooks contribute a weight of $q^{\alpha r}$ regardless of where
the rooks are. On the other hand, the weight contributed by the remaining
cells depends on the placement of the rooks, and varies as
$q^{st_{1}}q^{st_{2}}\cdots q^{st_{n-r}}$ for some $0\leq t_{1}\leq
t_{2}\leq\ldots\leq t_{n-r}\leq r$. Hence, these cells contribute

$\sum_{0\leq t_{1}\leq t_{2}\leq\ldots\leq t_{n-r}\leq
r}q^{st_{1}}q^{st_{2}}\cdots
q^{st_{n-r}}=\genfrac{[}{]}{0.0pt}{}{n}{r}_{q^{s}}\,.$

For 3, we consider rook placements using the default rule on the board
outlined by $(VU)^{n}V$ where the cells in the first row each have pre-weight
$\alpha$ while all the other cells have pre-weight 1. By Proposition 2, the
sum of the weight of all placements of $n-k$ rook on this board equals that of
$C^{R}_{n-k}(J^{\prime}_{n;\alpha};s)$ for any rule $R$.

A rook placement in $C_{n-k}(J^{\prime}_{n;\alpha};s)$ can then be formed as
follows. First, we place $n-r$ rooks in cells lying on rows $2,\ldots,n-1$ and
columns $2,\ldots,n$. The sum of the weights of all such placements is
$S_{s,q}[n,r]$. Then, place the remaining $r-k$ rooks on the first row. The
weight contributed by the cells in the first row containing rooks is
$\prod_{i=0}^{r-k-1}[\alpha+i(s-1)]_{q}$, and this is independent of the
choice of columns. On the other hand, the weight contributed by the cells in
the first row not containing rooks is $q^{\alpha t_{0}}q^{\alpha t_{1}}\cdots
q^{\alpha t_{r-k}}$ for some $t_{0}+t_{1}+\ldots+t_{r-k}=k$. Finally,

$\displaystyle\sum_{t_{0}+t_{1}+\ldots+t_{r-k}=k}$ $\displaystyle q^{\alpha
t_{0}}q^{\alpha t_{1}}\cdots q^{\alpha
t_{r-k}}=q^{\alpha}\genfrac{[}{]}{0.0pt}{}{r}{r-k}_{q^{s-1}}=q^{\alpha}\genfrac{[}{]}{0.0pt}{}{r}{k}_{q^{s-1}}\,.$

∎

######  Remark 3.

It is interesting to note that the individual terms on the $RHS$ of 2 and 3
are not equal for fixed $r$.

######  Theorem 5.

Let $n,m,k\in\mathbb{N}$. Then,

$\displaystyle
S_{s,q}[n+m,k]=\sum_{r=0}^{n}\sum_{j=0}^{m}S_{s,q}[m,j]\,q^{r(j(1-s)+sm)}\genfrac{[}{]}{0.0pt}{}{n}{r}_{q^{s}}S_{s,q}[r,k-j]\prod_{i=0}^{n-r-1}[j(1-s)+sm+si]_{q}\,.$
(4)

Consequently,

$B_{s,q}[n+m;x]=\sum_{r=0}^{n}\sum_{j=0}^{m}S_{s,q}[m,j]\,q^{r(j(1-s)+sm)}\genfrac{[}{]}{0.0pt}{}{n}{r}_{q^{s}}B_{s,q}[r;x]x^{j}\prod_{i=0}^{n-r-1}[j(1-s)+sm+si]_{q}\,.$
(5)

In particular,

$\displaystyle B_{s,q}[n+m]$
$\displaystyle=\sum_{r=0}^{n}\sum_{j=0}^{m}S_{s,q}[m,j]\,q^{r(j(1-s)+sm)}\genfrac{[}{]}{0.0pt}{}{n}{r}_{q^{s}}B_{s,q}[r]\prod_{i=0}^{n-r-1}[j(1-s)+sm+si]_{q}\,.$
(6)

###### Proof.

The number $S_{s,q}[n+m,k]$ equals the sum of the weights of all rook
placements in the set $C_{n+m-k}(J_{n+m;s})$. The rooks may be placed as
follows. First, place $m-j$ rooks in columns $2,\ldots,m$. The sum of the
weights of all such placements is $S_{s,q}[m,j]$. Next, place the remaining
$n+j-k$ rooks in columns $m+1,\ldots,n$. Due to the placement of the first
$m-j$ rooks, columns $m+1,\ldots,n$ form a board whose column pre-weights are
identical to those of $J^{\prime}_{n;j(1-s)+sm}$. Identity 4 follows from 2 of
Theorem 4 using the replacement $\alpha\rightarrow j(1-s)+sm,k\rightarrow
k-j$.

To obtain 5, multiply both sides of 4 by $x^{k}$ and take the sum over all
$0\leq k\leq n+m$. ∎

######  Theorem 6.

Let $n,m\in\mathbb{N}$. Then,

$\displaystyle
S_{s,q}[n+m,k]=\sum_{r=0}^{n}\sum_{j=0}^{m}S_{s,q}[m,j]q^{(j(1-s)+sm)(k-j)}\genfrac{[}{]}{0.0pt}{}{r}{k-j}_{q^{s-1}}S_{s,q}[n,r]\prod_{i=0}^{r-k+j-1}[j(1-s)+sm+i(s-1)]_{q}\,.$
(7)

Consequently,

$\displaystyle B_{s,q}[n+m;x]$
$\displaystyle=\sum_{r,j,k}S_{s,q}[m,j]q^{(j(1-s)+sm)(k-j)}\genfrac{[}{]}{0.0pt}{}{r}{k-j}_{q^{s-1}}S_{s,q}[n,r]x^{k}\prod_{i=0}^{r-k+j-1}[j(1-s)+sm+i(s-1)]_{q}\,.$
(8)

In particular,

$\displaystyle B_{s,q}[n+m]$
$\displaystyle=\sum_{r,j,k}S_{s,q}[m,j]q^{(j(1-s)+sm)(k-j)}\genfrac{[}{]}{0.0pt}{}{r}{k-j}_{q^{s-1}}S_{s,q}[n,r]\prod_{i=0}^{r-k+j-1}[j(1-s)+sm+i(s-1)]_{q}\,.$
(9)

###### Proof.

The proof is the same as that of Theorem except that we use Identity 3 of
Theorem 4. ∎

Taking cue from the identities in [12, Theorem 2.6], we derive another set of
convolution identities.

######  Theorem 7.

For $n,m,j\in\mathbb{N}$, we have

$\displaystyle S_{s,q}[n+1,m+j+1]$
$\displaystyle=\sum_{k=m}^{n}\sum_{r=0}^{n-k}S_{s,q}[k,m]\,q^{r((k-m)(s-1)+k+1)+k+(k-m)(s-1)}$
$\displaystyle\hskip
36.135pt\genfrac{[}{]}{0.0pt}{}{n-k}{r}_{q^{s}}S_{s,q}[r,j]\prod_{i=0}^{n-k-r-1}[(k-m)(s-1)+k+1+si]_{q}\,,$
(10) $\displaystyle S_{s,q}[n+1,m+j+1]$
$\displaystyle=\sum_{k=m}^{n}\sum_{r=0}^{n-k}S_{s,q}[k,m]\,q^{(j+1)(k+(k-m)(s-1))+j}$
$\displaystyle\hskip
36.135pt\genfrac{[}{]}{0.0pt}{}{r}{j}_{q^{s-1}}S_{s,q}[n-k,r]\prod_{t=0}^{r-j-1}[(k-m)(s-1)+k+1+t(s-1)]_{q}\,.$
(11)

###### Proof.

For a given rook placement in $C_{n-m-j}(J_{n+1};s)$, there exists a unique
$k$ such that there are $k-m$ rooks in columns $2,\ldots,k$, $0$ rooks in
column $k+1$ and $n-j-k$ rooks in columns $k+2,k+3,\ldots,n+1$. Indeed, if
there exists $k_{1}<k_{2}$ satisfying these properties, then there are
$k_{2}-k_{1}$ rooks in columns $k_{1}+1,k_{1}+2,\ldots,k_{2}$. This is
impossible since there is no rook in column $k_{1}+1$.

Suppose that $k-m$ rooks have been placed in columns $1,2,\ldots,k$. The
weight contributed by the $k+1$-th column is $q^{k+(k-m)(s-1)}$. Using Theorem
4 with the substitutions $\alpha\rightarrow k+1+(k-m)(s-1),n\rightarrow
n-k,k\rightarrow j,i\rightarrow r$, we obtain two different expressions for
the sum of the weights of the placements of $n-j-k$ rooks in columns
$k+2,k+3,\ldots,n+1$. ∎

######  Remark 4.

We briefly discuss some special cases.

By letting $s=0$ in 6, we get the following $q$-analogue of 1 which was
previously derived by Katriel [8]

$B_{q}[n+m]=\sum_{r=0}^{n}\sum_{j=0}^{m}S_{q}[m,j]q^{rj}\binom{n}{r}B_{q}[r][j]_{q}^{n-r}\,,$

where, $S_{0,q}[n,k]=S_{q}[n,k]$ and $B_{0,q}[n+m]=B_{q}[n+m]$ are the
$q$-Stirling number of the second kind and the $q$-Bell number, respectively.

A “dual” of Spivey’s identity was derived by Mezö [13] and is given by

$\displaystyle(n+m)!=\sum_{r=0}^{n}\sum_{j=0}^{m}c(m,j)m^{\overline{n-r}}\binom{m}{j}r!\,,$
(12)

where $a^{\overline{b}}=a(a+1)\cdots(a+b-1)$ and $c(n,k)$ is the Stirling
number of the first kind. One can check that6 reduces to 12 when $s=1,q=1$.
Using 9, we deduce another expression for $(n+m)!$, namely

$(n+m)!=\sum_{r,j,k}c(m,j)\binom{r}{k-j}c(n,r)m^{r-k+j}\,,$

or equivalently,

$(n+m)!=\sum_{r,j,k}c(m,j)\binom{r}{k}c(n,r)m^{k}\,.$

A simple combinatorial interpretation of the above identity was given by Mezö
[14]. Color the first $m$ elements red and the remaining $n$ elements blue.
Place the red elements into $j$ cycles in $c(m,j)$ ways and the blue elements
into $r$ cycles in $c(n,r)$ ways. Then, write the cycles in standard form,
i.e., such that the maximum element in each cycle is written first and the
cycles are arranged by increasing first elements. Pick $k$ cycles among the
$r$ blue cycles in $\binom{r}{k}$ ways. Insert the $k$ blue cycles among the
red elements such that if there are consecutive blue cycles after a red
element, we rewrite the blue cycles in standard form. This can be done in
$m^{k}$ ways. Finally, we remove the parentheses of the blue cycles which were
inserted and forget the colors.

Other identities for the the classical Stirling numbers and Bell numbers and
their $q$-analogues may be obtained as straight-forward corollaries of the
results in this section. For instance, Theorem 6 gives

$\displaystyle S_{q}[n+m,k]$
$\displaystyle=\sum_{r=0}^{n}\sum_{j=0}^{m}S_{q}[m,j]q^{j(k-j)}\binom{r}{k-j}_{1/q}S_{q}[n,r][j]_{q}[j-1]_{q}\cdots[k-r+1]_{q}\,,$

which appears to be new.

The proposition that follows generalizes Theorems 5, 6 and 7. As before, two
explicit forms, which we no longer state, can be obtained through an
application of Theorem 4. For the proofs, only the needed partition of the
boards are described.

######  Proposition 8.

There holds

$\displaystyle S_{s,q}[m_{1}+\cdots+m_{n},k]$
$\displaystyle=\sum_{\begin{subarray}{c}j_{1}+\cdots+j_{n}=k\\\
j_{1},\ldots,j_{n}\in\mathbb{N}\end{subarray}}\,\,\prod_{i=1}^{n}\sum\mbox{wt}(\phi)\,,$
$\displaystyle S_{s,q}[n+j-1,m_{1}+\cdots+m_{j}+j-1]$
$\displaystyle=\sum_{\begin{subarray}{c}k_{1}+\cdots+k_{j}=n\\\
k_{1},\ldots,k_{n}\in\mathbb{N}\end{subarray}}q^{\sum_{t=1}^{i-1}k_{t}+(k_{t}-m_{t})(s-1)}\prod_{i=1}^{j}\sum\mbox{wt}(\phi)\,.$

where in the first expression, the rightmost sum is taken over all rook
placements in
$C_{m_{i}-j_{i}}\left(J^{\prime}_{m_{i};\sum_{t=1}^{i-1}j_{t}(1-s)+sm_{t}};s\right)$
and in the second expression,
$C_{k_{i}-m_{i}}\left(J^{\prime}_{k_{i};\sum_{t=1}^{i-1}k_{t}+1+(k_{t}-m_{t})(s-1)};s\right)$.

###### Proof.

For the first identity: starting from the right, divide the the board
$J_{m_{1}+\cdots+m_{n}}$ into $n$ groups $G_{1},\ldots,G_{n}$ consisting of
$m_{1},\ldots,m_{n}$ adjacent columns, respectively. Then, for $1\leq i\leq
n$, place $m_{i}-j_{i}$ rooks on the columns in $G_{i}$.

For the second identity: use the fact that for every rook placement in
$C_{n-(m_{1}+\cdots+m_{j})}(J_{n+j-1})$, there exists a unique $j$-tuple
$(k_{1},\ldots,k_{j})$ satisfying $k_{1}+\cdots+k_{j}=n$ such that (a)
starting from the right, the board $J_{n+j-1}$ can be divided into $j$ groups
$G_{1},\ldots,G_{j}$ consisting of $k_{1},\ldots,k_{n}$ adjacent columns,
respectively, with one column not containing a rook in between, and (b) for
$1\leq i\leq j$, there are $k_{i}-m_{i}$ rooks on the columns in $G_{i}$. ∎

##  4 Extension to Type II $q$-analogues

Hsu and Shuie [7] introduced generalized Stirling numbers via the relations

$\displaystyle(x|\alpha)_{n}$
$\displaystyle=\sum_{k=0}^{n}S^{1}_{n,k}(\alpha,\beta,\rho)(x-\rho|\beta)_{k}$
$\displaystyle(x|\beta)_{n}$
$\displaystyle=\sum_{k=0}^{n}S^{2}_{n,k}(\alpha,\beta,\rho)(x+\rho|\alpha)_{k}$

where $(z|\alpha)_{0}=1$ and $(z|\gamma)_{n}=z(z-\gamma)\cdots(z-(n-1)\gamma)$
for every positive integer $n$. When $(\alpha,\beta,r)=(0,1,0)$,
$S^{1}_{n,k}(\alpha,\beta,\rho)$ and $S^{2}_{n,k}(\alpha,\beta,\rho)$ become
the usual Stirling number of the first kind and Stirling number of the second
kind, respectively. Since
$S^{1}_{n,k}(\alpha,\beta,\rho)=S^{2}_{n,k}(\beta,\alpha,-\rho)$, it suffices
to consider just one of these numbers, namely
$S^{1}_{n,k}(\alpha,\beta,\rho)$.

Let $[\gamma]_{p,q}=\frac{p^{\gamma}-q^{\gamma}}{p-q}$ for any
$\gamma\in\mathbb{R}$. The Type II $p,q$-analogues of these numbers were
introduced by Remmel and Wachs [15] by the following recursions ($\rho$
appears as $-r$ in [15])

$\displaystyle S^{1,p,q}_{n,k}(\alpha,\beta,\rho)$
$\displaystyle=q^{(k-1)\beta-(n-1)\alpha+\rho}S^{1,p,q}_{n-1,k-1}(\alpha,\beta,\rho)+p^{t-k\beta}[k\beta-(n-1)\alpha+\rho]S^{1,p,q}_{n-1,k}(\alpha,\beta,\rho)$
$\displaystyle S^{2,p,q}_{n,k}(\alpha,\beta,\rho)$
$\displaystyle=q^{-\rho+(k-1)\alpha-(n-1)\beta}S^{2,p,q}_{n-1,k-1}(\alpha,\beta,\rho)+p^{t+\rho-k\alpha}[k\alpha-\rho-(n-1)\beta]S^{2,p,q}_{n-1,k}(\alpha,\beta,\rho)$

where in addition
$S^{1,p,q}_{0,0}(\alpha,\beta,\rho)=S^{2,p,q}_{0,0}(\alpha,\beta,\rho)=1$ and
$S^{1,p,q}_{n,k}(\alpha,\beta,\rho)=S^{2,p,q}_{n,k}(\alpha,\beta,\rho)=0$ if
$k=0$ or $k>n$. Since
$S^{1,p,q}_{n,k}(\alpha,\beta,\rho)=S^{2,p,q}_{n,k}(\beta,\alpha,-\rho)$, we
also consider just one of these numbers, which we choose to be
$S^{1,p,q}_{n,k}(\alpha,\beta,\rho)$. When $(\alpha,\beta,\rho)=(0,1,0)$ and
$p=q=1$, $S^{1,p,q}_{n,k}(\alpha,\beta,\rho)=S(n,k)$.

Remmel and Wachs also gave combinatorial interpretations for these numbers for
a certain choice of parameters. Our goal in this section is to introduce an
alternative rook interpretation for $S^{1,1,q}_{n,k}(\alpha,\beta,\rho)$ which
adapts the method used in Section 3. To do this, we will derive the analogue
of Theorem 5 and Theorem 6. We leave the corresponding versions of the other
identities in Section 3 to the interested reader.

Let $c,d\in\mathbb{R}$ and denote by $J^{c,d}_{n}$ the board outlined by
$(VU)^{n}$ such that each bottom cell has default pre-weight $d$ and all other
cells have default pre-weight $c+d$. Let $S^{c,d}_{s,q}[n,k]$ denote the sum
of the weights of all placements of $n-k$ rooks on $J^{c,d}_{n}$.

######  Remark 5.

Since $J^{1,0}_{n}=J_{n}$, it follows that $S^{1,0}_{s,q}[n,k]=S_{s,q}[n,k]$.

######  Remark 6.

A “special case” of $J^{c,d}_{n-1}$ with $c=m,d=0$ are the $m$-jump boards
$Jb_{n,m}$ [5, Example 3] which have column heights $0,m,2m,\ldots,(n-1)m$.
One sees that these two boards have identical column pre-weights.

######  Proposition 9.

The number $S^{c,d}_{s,q}[n,k]$ satisfies the recurrence

$S^{c,d}_{s,q}[n,k]=q^{c(n-1)+d+(n-k)(s-1)}S^{c,d}_{s,q}[n-1,k-1]+[c(n-1)+d+(n-k-1)(s-1)]_{q}S^{c,d}_{s,q}[n-1,k]$

Hence,

$S^{\beta-\alpha,-\rho}_{-\beta+1,q}[n,k]=S^{1,1,q}_{n,k}(\alpha,\beta,\rho)\,.$

###### Proof.

The proof is similar to that of Proposition 1.∎

######  Theorem 10.

For $n,m,k\in\mathbb{N}$, we have

$S^{c,d}_{s,q}[n+m,k]=\sum_{r=0}^{n}\sum_{j=0}^{m}S^{c,d}_{s,q}[m,j]\,q^{r(d+mc+(m-j)(s-1))}\\\
\genfrac{[}{]}{0.0pt}{}{n}{r}_{q^{c+s-1}}S^{c,0}_{s,q}[r,k-j]\prod_{i=0}^{n-r-1}[d+mc+(m-j)(s-1)+(c+s-1)i]_{q}\,.$

This implies that

$S^{1,1,q}_{n,k}(\alpha,\beta,\rho)=\sum_{r=0}^{n}\sum_{j=0}^{m}S^{1,1,q}_{m,j}(\alpha,\beta,\rho)\,q^{r(\beta
j-\alpha
m-\rho)}\genfrac{[}{]}{0.0pt}{}{n}{r}_{q^{-\alpha}}S^{1,1,q}_{r,k-j}(\alpha,\beta,0)\prod_{i=0}^{n-r-1}[\beta
j-\rho-\alpha(m+i)]_{q}\,.$

###### Proof.

We proceed as in Theorem 5. The $m-j$ rooks in columns $2,\ldots,m$ explains
the factor $S^{c,d}_{s,q}[m,j]$. We then determine the sum of the weights of
all placements of $n+j-k$ rooks in the remaining columns.

The placement of $m-j$ rooks adds $(m-j)(s-1)$ to the total pre-weights of
columns $m+1,\ldots,n$. These columns form a board whose column pre-weights
are identical to the board outlined by $(VU)^{n}V$ such that the bottom cells
have pre-weight $d+mc+(m-j)(s-1)$ and the other cells have pre-weight $c$. For
convenience, let us denote such board by $B$.

Note that Lemma 3 is a statement regarding only the placement of rooks and
thus, it holds true regardless of the default pre-weights of the cells and the
rule used. Hence, we can use this lemma on $B$ with rule $\mathcal{R}$,
reasoning as in Theorem 5.

Let $\phi\in C^{\mathcal{R}}_{n+j-k}(B;s)$. Again, by Lemma 3, we can write
$\phi$ uniquely as $(\mathcal{C}_{\phi},\phi-\mathcal{C}_{\phi})$. Note that
$\phi-\mathcal{C}_{\phi}$ forms a rook placement in $C_{r-k+j}(J^{c,0}_{r};s)$
and not in $C_{r-k+j}(J^{c,d}_{r};s)$! This is because
$\phi-\mathcal{C}_{\phi}$ does not include any bottom cell of $\phi$.

The board $B$ is similar to $J^{\prime}_{n;\alpha}$ with
$\alpha=d+mc+(m-j)(s-1)$, the only difference is that the non-bottom cells of
$B$ have pre-weight $c$. Hence, a cell with pre-weight $c$ which receives an
additional pre-weight due to the placement of a rook following rule
$\mathcal{R}$ will have pre-weight $c+s-1$. ∎

Define
$B^{1,1,q}_{n;x}(\alpha,\beta,\rho)=\sum_{k=0}^{n}S^{1,1,q}_{n}(\alpha,\beta,\rho)$.

######  Corollary 11.

For $n,m,\in\mathbb{N}$, we have

$B^{1,1,q}_{n+m;x}(\alpha,\beta,\rho)=\sum_{r=0}^{n}\sum_{j=0}^{m}S^{1,1,q}_{m,j}(\alpha,\beta,\rho)\,q^{r(\beta
j-\alpha
m-\rho)}\genfrac{[}{]}{0.0pt}{}{n}{r}_{q^{-\alpha}}B^{1,1,q}_{r;x}(\alpha,\beta,0)\prod_{i=0}^{n-r-1}[\beta
j-\rho-\alpha(m+i)]_{q}\,.$ (13)

Letting $q=1$ produces a generalization of Spivey’s identity, which is a
variant of Xu’s result [18, Corollary 8]. Also, Mezö [13, Theorem 2] obtained
a special case of 13 for $(\alpha,\beta,\rho)=(0,1,\rho)$.

######  Theorem 12.

Let $n,m\in\mathbb{N}$. Then,

$S^{c,d}_{s,q}[n+m,k]=\sum_{r=0}^{n}\sum_{j=0}^{m}S^{c,d}_{s,q}[m,j]q^{(d+mc+(m-j)(s-1))(k-j)}\\\
\genfrac{[}{]}{0.0pt}{}{r}{k-j}_{q^{s-1}}S^{c,0}_{s,q}[n,r]\prod_{i=0}^{r-k+j-1}[d+mc+(m-j)(s-1)+i(s-1)]_{q}\,.$

This gives

$S^{1,1,q}_{n+m,k}(\alpha,\beta,\rho)=\sum_{r=0}^{n}\sum_{j=0}^{m}S^{1,1,q}_{m,j}(\alpha,\beta,\rho)q^{(\rho+\alpha
m-\beta j)(j-k)}\\\
\genfrac{[}{]}{0.0pt}{}{r}{k-j}_{q^{-\beta}}S^{1,1,q}_{n,r}(\alpha,\beta,0)\prod_{i=0}^{r-k+j-1}[\rho+\alpha
m-\beta(j+i)]_{q}\,.$ (14)

###### Proof.

Instead of using the board $J^{c,d}_{n}$, we use the board outlined by
$(VU)^{n}$ such that the bottom cells of columns $2,\ldots,m$ have pre-weight
$c+d$, the cells in row $m$ have pre-weight $c+d$, and all other cells have
pre-weight $c$. We then proceed as in the proof of Theorem 6. ∎

##  5 Some Remarks

Mezö’s dual of Spivey’s identity is based on the relation
$n!=\sum_{k=0}^{n}c(n,k)$. However, it is not true that
$[n]_{q}!=\sum_{k=0}^{n}c_{q}[n,k]$, where $c_{q}[n,k]=S_{1,q}[n,k]$. It would
be interesting to derive an expression for $[n+m]_{q}!$ similar to Mezö’s for
$q$-Stirling numbers of the first kind.

We also ask if a modification of the rook model for the Type II generalized
$q$-Stirling numbers can be made to given another combinatorial interpretation
for both the Type I and Type II $p,q$-analogues. It is also possible that the
techniques outlined in this paper may be modified to obtain the corresponding
identities for the generalized Stirling numbers that arise from other rook
models, such as those in [15] and [3].

Lastly, we note that other forms of convolution identities are given in [1,
Theorem 1] and [12, page 5]. Can these identities and their $S_{s,q}[n,k]$
versions be proved using partitions on rook boards?

## References

  * [1] T. Agoh and K. Dilcher, Generalized convolution identities for Stirling numbers of the second kind, INTEGERS 8 (2008), #A25.
  * [2] H. Belbachir and M. Mihoubi, A generalized recurrence for Bell polynomials: An alternate approach to Spivey and Gould-Quaintance formulas, Europ. J. Combin. 30 (2009), 1254-1256.
  * [3] K. Briggs, A rook theory model for the generalized $p,q$-Stirling numbers of the first and second kind, in Formal Power Series and Algebraic Combinatorics, San Diego, California (2006).
  * [4] R. Corcino, R. Celeste and K. J. Gonzales, Generalized $q$-Stirling numbers and normal ordering, _preprint_ arXiv:1407.3343.
  * [5] J. Goldman and J. Haglund, Generalized rook polynomials, J. Combin. Theory Ser. A 91 (2000), 509-530.
  * [6] H. W. Gould and J. Quaintance, Implications of Spivey’s Bell number formula, J. Integer Seq. 11 (2008), 08.3.7.
  * [7] L. Hsu and P.J. Shiue, A unified approach to generalized Stirling numbers, Adv. App. Math. 20 (1998), 366-384.
  * [8] J. Katriel, On a generalized recurrence for Bell numbers, J. Integer Seq. 11 (2008), Article 08.3.8.
  * [9] M. S. Maamra and M. Mihoubi, The $(r_{1},\ldots,r_{p})$-Bell polynomials, INTEGERS 14 (2004), #A34.
  * [10] T. Mansour, M. Schork and M. Shattuck, On a new family of generalized Stirling and Bell numbers, Electronic J. Combin. 18 (2011), #P77.
  * [11] T. Mansour, M. Schork and M. Shattuck, The generalized Stirling and Bell numbers revisited, J. Integer Seq. 15 (2012), Article 12.8.3.
  * [12] A. Médicis and P. Leroux, Generalized Stirling numbers, convolution formulae and $p,q$-analogues, Canad. J. Math 47 (1995), 474- 499\.
  * [13] I. Mezö, The dual of Spivey’s Bell number formula, J. Integer Seq. 15(2) (2012), Article 12.2.4.
  * [14] I. Mezö, personal communication.
  * [15] J. B. Remmel and M. Wachs, Rook Theory, Generalized Stirling numbers and $(p,q)$-analogues, Electronic J. Combin. 11 (2004), #R84.
  * [16] M. Spivey, A generalized recurrence for Bell numbers, J. Integer Seq. 11 (2008), Article 08.2.5.
  * [17] A. Varvak, Rook numbers and the normal ordering problem, J. Combin. Theory Ser. A 112 (2005), 292-307.
  * [18] A. Xu, Extensions of Spivey’s Bell number formula, Electronic J. Combin. 19(2) (2012), #P6.

|  |  |  |
---|---|---|---|---
|  |  |  |
|  |  |  |
Figure 1: The board outlined by $UVUUUVVU$. 0 | s | $\bullet_{s}$ | $\bullet_{\alpha}$
---|---|---|---
0 | s | α |
0 | α |  |
$\bullet_{\alpha}$ |  |  |
Figure 2: A rook placement on $J^{\prime}_{4,\alpha}$ using rule $\mathcal{R}$. 0 | 2s-1 | $\bullet_{s}$ | $\bullet_{1}$ |
---|---|---|---|---
0 | 1 | 1 |  |
0 | 1 |  |  |
$\bullet_{1}$ |  |  |  |
Figure 3: A rook placement on $J_{5}$.


# Efficient Causal Graph Discovery Using Large Language Models

Thomas Jiralerspong  Xiaoyin Chen  Yash More  Vedant Shah  Yoshua Bengio

###### Abstract

We propose a novel framework that leverages LLMs for full causal graph
discovery. While previous LLM-based methods have used a pairwise query
approach, this requires a quadratic number of queries which quickly becomes
impractical for larger causal graphs. In contrast, the proposed framework uses
a breadth-first search (BFS) approach which allows it to use only a linear
number of queries. We also show that the proposed method can easily
incorporate observational data when available, to improve performance. In
addition to being more time and data-efficient, the proposed framework
achieves state-of-the-art results on real-world causal graphs of varying
sizes. The results demonstrate the effectiveness and efficiency of the
proposed method in discovering causal relationships, showcasing its potential
for broad applicability in causal graph discovery tasks across different
domains.

LLM, large language models, causality, machine learning, ICML

##  1 Introduction

Recent advances in large language models (LLMs) have significantly improved AI
capabilities (OpenAI, 2023). LLMs have access to a vast amount of training
data, allowing them to acquire a substantial amount of knowledge, ranging from
common-sense to specific domains, such as math and science. Moreover, recent
papers have argued that complex behaviors, such as writing code, generating
long stories, and even some reasoning capabilities, can emerge from such
large-scale training (Wei et al., 2022; Rozière et al., 2023; Zhao et al.,
2023; Yao et al., 2023). Impressively, LLMs can use their knowledge to
generate plausible explanations given observations and even reason on
counterfactuals (Wang et al., 2023; Li et al., 2023; Kıcıman et al., 2023).
Given these emergent abilities to understand and learn relationships from
textual corpora, a natural question arises: Can LLMs reason about causal
relationships?

Causal discovery involves identifying the causal relationships between
different variables (Pearl, 2009; Peters et al., 2017). This process usually
results in a Directed Acyclic Graph (DAG), with edges representing the
existence and direction of causal influences. This causal graph represents the
relationships between variables which govern the data generation process. The
causal graph also allows us to determine how changes in one variable might
influence others, and serves as the foundation for further analysis related to
specific tasks, including effect inference, prediction, or attribution.
Ensuring the correctness of this causal graph is essential for the reliability
of any subsequent analysis.

Traditional causal discovery frameworks (Bello et al., 2023; Zheng et al.,
2018; Meek, 2023) rely on numerical observations and use statistical or
continuous optimization based methods to unravel confounding and causal
factors. While such rigorous numerical and abstract reasoning abilities remain
challenges for LLMs (Gendron et al., 2023; Xu et al., 2023; Imani et al.,
2023), recent works have shown that LLMs are capable of leveraging their
internal knowledge to perform commonsense and simple logical as well as causal
reasoning (Wei et al., 2022; Zhao et al., 2023; Yao et al., 2023; Kıcıman et
al., 2023). Inspired by these results, the proposed framework utilizes LLMs to
respond to causal queries solely based on metadata, such as variable names and
descriptions, without accessing the numerical observations. This is similar to
the metadata and knowledge-based reasoning that is done by human domain
experts when they construct causal graphs.

Previous works (Kıcıman et al., 2023; Choi et al., 2022; Long et al., 2023b)
have shown that LLMs are capable of inferring pairwise causal relationships
given only metadata. However, extending this pairwise framework to full graph
discovery presents additional challenges. One straightforward approach which
previous works use is to perform pairwise queries over all possible pairs of
variables. However, this method is not scalable since it requires $n\choose 2$
total queries where $n$ is the number of variables. Thus, the number of
queries required grows quadratically with the number of causal variables. As a
result, the experiments in these works are limited to very small causal
graphs.

To address the above problems, we propose a novel framework to discover full
causal graphs with the internal knowledge of LLMs that guarantees the DAG
constraint and avoids performing all pairwise queries (Figure 1). The key
insight is that we leverage the DAG property of the resulting graph and
construct the graph through a breadth-first search (BFS). Intuitively, node
expansions in BFS correspond to prompting a LLM to select the variables caused
by the current node. To efficiently construct causal graphs, we iteratively
expand nodes with a LLM in BFS order and add the discovered edges. The
proposed method achieves $O(n)$ query complexity compared to the $O(n^{2})$
query complexity of the pairwise method.

To the best of our knowledge, we present the first LLM-based causal graph
discovery method applicable to larger graphs. The proposed method is much more
efficient than previous LLM-based methods, and does not require observational
data (unlike numerical methods), while still achieving state-of-the-art or
competitive performance on three graphs of varying sizes. In addition, we show
how observational data can easily be incorporated when available to improve
the proposed method’s performance. We provide code to reproduce all
experiments at this URL.

Figure 1: Proposed framework for full graph discovery with LLMs

##  2 Related Work

Causal graph discovery has been extensively studied in many settings (Zanga &
Stella, 2023). Until recently, most proposed methods were statistical methods
that used observational data to infer causal relationships between variables
(Bello et al., 2023; Zheng et al., 2018; Meek, 2023; Spirtes et al., 1993).
Recently, many works have proposed to use LLMs for causal graph discovery
(Kıcıman et al., 2023; Choi et al., 2022; Long et al., 2023b). These LLM-based
methods do not rely on observational data. Instead, they utilize metadata
associated with the input variables, specifically variable descriptions in
natural language. This approach closely resembles how human experts construct
causal graphs using their domain knowledge. In the rest of this section, we
provide a brief overview of popular causal graph discovery algorithms (both
statistical and LLM-based).

###  2.1 Statistical Methods

####  2.1.1 Peter-Clark Algorithm (PC) (Spirtes & Glymour, 1991)

PC is a constraint-based algorithm which starts by defining a complete
undirected graph over the variables. It then performs a series of conditional
independence tests and uses these tests to remove the associated edges from
the complete graph. Finally, it orients the edges remaining in the graph using
rules that take into account acyclicity and other causal assumptions.

####  2.1.2 Greedy Equivalence Search (GES) (Meek, 2023)

GES is a score-based causal discovery algorithm, meaning that it seeks to
optimize some score that is a function of the graph. In our case, the score is
the Bayesian Information Criterion (BIC) (Schwarz, 1978). GES begins with an
empty graph. It then alternates between 2 phases: the forward phase and the
backward phase. During the forward phase, the edge that, once added, makes the
graph’s score go up the most is repeatedly added until there is no edge that
increases the score. During the backward phase, the edge that, once removed,
makes the graph’s score go up the most, is repeatedly removed until there is
no edge that increases the score. The algorithm finishes once a local maximum
is reached during the backward phase.

####  2.1.3 NOTEARS (Zheng et al., 2018)

Non-combinatorial Optimization via Trace Exponential and Augmented lagrangian
for Structure learning (NOTEARS) assumes that the data is generated linearly
according to the model $\textbf{X}=\textbf{X}W+N$ where $N$ is a noise vector
with independent entries, X is a $n\times d$ observation matrix containing $n$
i.i.d observations of $d$ variables and $W_{ij}$ can be viewed as the
probability that variable $V_{i}$ has a direct causal influence on $V_{j}$, or
equivalently, that there is a directed edge from node $i$ to node $j$ in the
causal graph between the $d$ variables.

NOTEARS thus solves the following constrained optimization problem using the
augmented Lagrangian method:

$\hat{W}=argmin_{W}\frac{1}{2n}||\textbf{X}-\textbf{X}W||^{2}_{F}+\lambda||W||_{1}$

s.t.

$exp(W\circ W)-d=0$

where $\lambda||W||_{1}$ is a regularizer and the differentiable constraint
$exp(W\circ W)-d=0$ is added to ensure that the resulting causal graph is a
directed acyclic graph (DAG) as desired.

The entries $W_{ij}$ in $W$ which are above a certain threshold are then added
to the output causal graph as directed edges from node $i$ to node $j$.

####  2.1.4 DAGMA (Bello et al., 2023)

Directed Acyclic Graphs via M-matrices for Acyclicity (DAGMA) solves the same
optimization problem as NOTEARS but using a different acyclicity constraint
based on the log-determinant function.

###  2.2 LLM-based methods

####  2.2.1 Pairwise Method

Kıcıman et al. (2023); Choi et al. (2022); Long et al. (2023b) use pairwise
queries to infer the causal relationship between 2 variables at a time. Given
variables $A$ and $B$, there are three possible causal relationships,
$A\rightarrow B$, $A\leftarrow B$, and no relationship. This query can be
formulated as a question in natural language with the following prompt:

    <A>: <Description A>
    <B>: <Description B>

    Given the above information, which of
    the following is the most likely:
    A. Changing <A> causes a change in <B>
    B. Changing <B> causes a change in <A>
    C. There is no causal relationship
    between <A> and <B>

To extend to full graph discovery, Kıcıman et al. (2023) propose to perform
pairwise queries over all possible pairs. However, this method is not scalable
since it requires $n\choose 2$ $\in O(n^{2})$ total queries, where $n$ is the
number of variables, so they are only able to perform experiments on a small
causal graph with 12 variables.

##  3 Methods

###  3.1 Main Method

We now propose a framework for efficient full causal graph discovery using a
LLM. Instead of querying the LLM for each pair of nodes, we instruct the LLM
to identify all potential effects of a single variable in one completion. This
approach reduces the amount of queries required from $O(n^{2})$ to $O(n)$
since visiting each node once is sufficient to determine all edges originating
from that node.
With this new prompting method, we can treat each query as a node expansion in
a BFS algorithm and gradually construct the causal graph by traversing it
using BFS. To perform the BFS, we also need to establish the order in which
nodes are explored. In this regard, we draw inspiration from the DAG property
of causal graphs and utilize their topological orders for BFS. Specifically,
the proposed framework consists of three stages: 1) an Initialization stage
wherein the LLM is prompted to identify variables that are not causally
affected by any other variables in the graph, 2) an Expansion stage wherein
the LLM is instructed to find variables caused by the current node, and 3) an
Insertion stage where the variables proposed by the LLM are added to the BFS
queue and the edges proposed by the LLM are checked to make sure they do not
cause cycles before being inserted. We start with the initialization stage,
then alternate between the expansion and insertion stages until all nodes have
been visited. The pseudocode for the proposed algorithm is shown in Algorithm
1.
The initialization stage determines variables that are not caused by any other
variables. These selected variables are added to the the BFS queue to be used
later as starting points for the BFS. We use the following initialization
prompt to provide the LLM with descriptions of all variables and ask the model
to select independent variables. Note that the first paragraph of this prompt
is task-specific. Here we use the prompt for experiments on the Neuropathic
Pain Dataset (see section 4.1) as an example:

    You are a helpful assistant to a
    neuropathic pain diagnosis expert.
    The following factors are key
    variables related to neuropathic pain
    diagnosis which have various causal
    effects on each other. Our goal is
    to construct a causal graph between
    these variables.

    <A>: Description of variable A
    <B>: Description of variable B
    ...

    Now you are going to use the data to
    construct a causal graph. You will
    start with identifying the variable(s)
    that are unaffected by any other
    variables.
    Think step by step. Then, provide your
    final answer (variable names only)
    within the tags <Answer>...</Answer>.

The expansion stage expands each node in the search tree by prompting the LLM
to select all variables that are caused by the current node. Given a queue of
nodes that have not been visited yet (initialized by the initialization
stage), we visit the nodes in a first-in-first-out order, i.e., the BFS order.
We use the following expansion prompt to provide the LLM with the current
graph structure (i.e. the causal relationships already predicted) and let the
model predict a list of affected variables:

    Given <Independent Variables> is(are)
    not affected by any other variable and
    the following causal relationships:

    A causes B, C, D
    C causes D, E

    Select the variables that are caused by
    <Currently Visited Node>.
    Think step by step. Then, provide your
    final answer (variable names only)
    within the tags <Answer>...</Answer>.

In the insertion stage, for each variable predicted by the model, we first add
it to the queue of nodes to be visited. Next, for each edge from the current
node to its expanded children, we perform a cycle check before adding this
edge to the predicted causal graph. This check is necessary because, although
the LLM provides the current graph structure, we observe that it often
generates edges that violate the DAG constraint. To check for cycles, we
perform a depth-first search in the predicted graph with the edge to be added.
If adding the edge creates a cycle, we do not add it, thus guaranteeing that
the resulting graph is a DAG.
We alternate between the expansion and insertion stages until all nodes have
been visited. The number of queries required for the proposed algorithm is
$O(n)$, as BFS only visits each node once. This is in contrast to the pairwise
method (Sec. 2.2), which requires $O(n^{2})$ queries.

###  3.2 Adding Statistics to the Prompt

We also experiment with adding relevant statistics to the prompt in the form
of the Pearson correlation coefficients between the variables involved in the
prompt. Even though correlation does not imply causation, we hypothesize that
adding such statistics to the prompt helps the LLM by giving it statistical
relationships which it can use when its internal knowledge is not enough to
determine the causal relationship between 2 variables.

For the proposed method, this is done by adding the following text to the
expansion stage prompt after the list of causal relationships:

    Additionally, the Pearson Correlation
    Coefficients between <Currently Visited
    Node> and the other variables are as
    follows:
    <A>: <Pearson Correlation Coefficient
         between <Currently Visited Node>
         and A>
    <B>: <Pearson Correlation Coefficient
         between <Currently Visited Node>
         and B>
    ...

For the pairwise method (Sec. 2.2), the following text is added to the prompt
after the name and description of the variables:

    Additionally, the Pearson Correlation
    Coefficient between <A> and <B> is
    <Correlation Coefficient>.

Lastly, it is important to mention that the search process takes place within
a single multi-turn chat environment. This enables the LLM to base each prompt
on the preceding chat history.

Algorithm 1 BFS with LLMs

LM $p_{\theta}$, descriptions of variables $X$, initial variable selector
$I()$, expansion generator $E()$, cycle checker $CheckCycle()$

$G\leftarrow\\{\\}$ $\triangleright$ Create an empty graph to store the
result.

$frontier,visited\leftarrow I(p_{\theta},X)$ $\triangleright$ With
initialization prompt.

while $frontier$ is not empty do

$toVisit\leftarrow frontier\left[0\right]$

$frontier.remove(toVisit)$

$visited.add(toVisit)$

for $node$ in $E(p_{\theta},G)$ do $\triangleright$ Expand with expansion
prompt.

if not $CheckCycle(G,toVisit,node)$ then

$\triangleright$ Check if adding $toVisit\rightarrow node$ will create cycle.

$G.add((toVisit,node))$

end if

if $node$ not in $frontier\cup visited$ then

$frontier.add(node)$

end if

end for

end while

return $G$

##  4 Experiments and Results

###  4.1 Experimental Setup

We use the GPT-4 API (OpenAI, 2023) for all of our experiments. We evaluate
the proposed approach on three causal graphs of varying sizes:

  1. 1.

Asia (Lauritzen & Spiegelhalter, 1988). The Asia causal graph is a small
causal graph that models lung cancer in patients having potentially recently
visited Asia. The ground-truth graph has 8 nodes and 8 edges. Some examples of
variables in the graph are Visit to Asia, Smoking, Bronchitis, and Lung
Cancer.

  2. 2.

Child (Spiegelhalter et al., 1993). The Child causal graph is a medium-sized
causal graph which models congenital heart disease in newborn babies. The
ground-truth graph has 20 nodes and 25 edges. Some examples of variables in
the graph are Birth Asphyxia, Lung Flow, and Chest XRay.

  3. 3.

Neuropathic Pain (Tu et al., 2019) The Neuropathic Pain causal graph is a very
large causal graph consisting of 221 nodes and 770 edges. It contains
relationships between different nerves and the associated symptoms that
patients express. It mainly consists of symptom diagnoses, pattern diagnoses,
and pathophysiological diagnoses. The symptom diagnoses describe the
discomfort of patients. The pattern diagnoses identify nerve roots that lead
to symptom patterns. The pathophysiological diagnoses identify
discoligamentous injuries.

For the Neuropathic Pain causal graph, we use the same prompts as described in
Section 3.1. We use similar prompts with appropriate task specific
modifications for the other two datasets.

We get the Asia and Child causal graphs and data from the bnlearn package
(Scutari, 2010), and their descriptions from Long et al. (2023a). We get the
Neuropathic Pain causal graph data and descriptions from Tu et al. (2019).

We compare the proposed method, as well as the proposed method plus
observational statistics from 1000 and 10000 samples, against the numerical
methods of GES, PC, NOTEARS, and DAGMA (all discussed in section 2.2), and
against pairwise queries, and pairwise queries plus observational statistics
from 1000 and 10000 samples. We use the implementations of GES and PC from the
causal-learn package (Zheng et al., 2023), the implementation of NOTEARS from
Zheng et al. (2018), and the implementation of DAGMA from Bello et al. (2023).
For the numerical methods, we conduct experiments in the low (100 samples),
medium (1000 samples), and high (10000 samples) data regimes.

For NOTEARS and DAGMA, we try $\lambda$ values in $\\{0,0.005,0.1,0.2\\}$ and
report the results for the best performing $\lambda$ value.

For the LLM methods, we use the newest gpt-4-0125-preview checkpoint. We try
sampling temperatures in $\\{0,0.5,0.7,1.0\\}$ and report the results for the
best temperature setting.

###  4.2 Metrics

Causal graph discovery can be viewed as a classification task where each pair
of variables (A,B) must be classified as either A causes B, B causes A, or
Neither. We therefore report standard classification metrics such as accuracy,
precision, recall, and F score.

Following Kıcıman et al. (2023), we also report the Normalized Hamming
Distance (NHD). In a graph with $m$ nodes, the NHD is given by
$\sum_{i,j=1}^{m}\frac{1}{m^{2}}\mathbbm{1}_{G_{ij}\neq G^{\prime}_{ij}}$,
i.e. the number of edges that are present in one graph and not the other,
divided by the total number of all possible edges. One problem with the NHD is
that since most graphs are relatively sparse, trivially outputting 0 predicted
edges achieves a relatively good NHD. For example, for the Child graph, which
has 20 nodes and 25 edges, if we output 0 predicted edges then we achieve a
NHD of $\frac{1}{(20^{2})}25=0.0625$, which is better than the NHD achieved by
all the methods. Intuitively, predicting an edge is much riskier when using
this metric, so algorithms that predict more edges are punished. To remedy
this and take into account the number of predicted edges, again following
Kıcıman et al. (2023), we compute the ratio between the NHD and the baseline
NHD of an algorithm that outputs the same number of edges but with all of them
being incorrect. The lower the ratio, the better the method does compared to
the worst baseline that outputs the same number of edges. We report the NHD,
number of predicted edges (NPE), baseline NHD, and ratio between NHD and
baseline NHD for all algorithms.

As our main metrics for comparison, we use the ratio between NHD and baseline
NHD as well as the F score since these metrics are the least biased by the
number of predicted edges.

---
Method | Acc. $\uparrow$ | Prec. $\uparrow$ | Recall $\uparrow$ | F Score $\uparrow$ | NPE | NHD $\downarrow$ | Ref. NHD | NHD Ratio $\downarrow$
Numerical Methods |  |  |  |  |  |  |  |
100 Samples |  |  |  |  |  |  |  |
GES | 0.23 | 0.38 | 0.38 | 0.38 | 8 | 0.28 | 0.44 | 0.63
PC | 0.33 | 0.75 | 0.38 | 0.5 | 4 | 0.24 | 0.48 | 0.5
NOTEARS ($\lambda=0.02$) | 0.5 | 1 | 0.5 | 0.67 | 4 | 0.11 | 0.33 | 0.33
DAGMA ($\lambda=0.02$) | 0.22 | 0.67 | 0.25 | 0.36 | 3 | 0.28 | 0.44 | 0.64
1000 Samples |  |  |  |  |  |  |  |
GES | 0.67 | 0.86 | 0.75 | 0.8 | 7 | 0.06 | 0.31 | 0.2
PC | 0.5 | 0.71 | 0.63 | 0.67 | 7 | 0.1 | 0.31 | 0.33
NOTEARS ($\lambda=0.005$) | 0.44 | 0.8 | 0.5 | 0.62 | 5 | 0.10 | 0.27 | 0.38
DAGMA ($\lambda=0.01$) | 0.5 | 1 | 0.5 | 0.67 | 4 | 0.11 | 0.33 | 0.33
10000 Samples |  |  |  |  |  |  |  |
GES | 0.7 | 0.78 | 0.88 | 0.82 | 9 | 0.06 | 0.35 | 0.18
PC | 0.55 | 0.67 | 0.75 | 0.71 | 9 | 0.1 | 0.35 | 0.29
NOTEARS ($\lambda=0.005$) | 0.44 | 0.8 | 0.5 | 0.62 | 5 | 0.14 | 0.36 | 0.38
DAGMA ($\lambda=0.005$) | 0.5 | 1 | 0.5 | 0.67 | 4 | 0.11 | 0.33 | 0.33
LLM Methods |  |  |  |  |  |  |  |
Pairwise | 0.33 | 0.35 | 0.88 | 0.5 | 20 | 0.22 | 0.44 | 0.5
Pairwise + 1000 samples | 0.38 | 0.43 | 0.75 | 0.54 | 14 | 0.20 | 0.45 | 0.45
Pairwise + 10000 samples | 0.47 | 0.47 | 1 | 0.64 | 17 | 0.14 | 0.39 | 0.36
Ours | 0.88 | 1 | 0.88 | 0.93 | 7 | 0.016 | 0.23 | 0.067
Ours + 1000 samples | 0.8 | 0.8 | 1 | 0.89 | 10 | 0.031 | 0.28 | 0.11
Ours + 10000 samples | 0.8 | 0.8 | 1 | 0.89 | 10 | 0.031 | 0.28 | 0.11
Table 1: Results on the Asia causal graph (8 nodes, 8 edges). The proposed
method with no observational statistics outperforms all other methods with a F
score of 0.93 and a NHD ratio of 0.067.
---
Method | Acc. $\uparrow$ | Prec. $\uparrow$ | Recall $\uparrow$ | F Score $\uparrow$ | NPE | NHD $\downarrow$ | Ref. NHD | NHD Ratio $\downarrow$
Numerical Methods |  |  |  |  |  |  |  |
100 Samples |  |  |  |  |  |  |  |
GES | 0.24 | 0.41 | 0.36 | 0.38 | 22 | 0.11 | 0.18 | 0.62
PC | 0.19 | 0.37 | 0.28 | 0.32 | 19 | 0.12 | 0.17 | 0.68
NOTEARS ($\lambda=0.01$) | 0.14 | 0.26 | 0.24 | 0.25 | 23 | 0.1 | 0.13 | 0.75
DAGMA ($\lambda=0.005$) | 0.14 | 0.2 | 0.32 | 0.24 | 41 | 0.13 | 0.17 | 0.76
1000 Samples |  |  |  |  |  |  |  |
GES | 0.31 | 0.41 | 0.47 | 0.44 | 34 | 0.08 | 0.15 | 0.53
PC | 0.29 | 0.38 | 0.56 | 0.45 | 37 | 0.085 | 0.155 | 0.55
NOTEARS ($\lambda=0.005$) | 0.2 | 0.41 | 0.28 | 0.33 | 17 | 0.086 | 0.13 | 0.67
DAGMA ($\lambda=0.005$) | 0.24 | 0.43 | 0.36 | 0.39 | 21 | 0.07 | 0.12 | 0.61
10000 Samples |  |  |  |  |  |  |  |
GES | 0.42 | 0.45 | 0.84 | 0.58 | 47 | 0.08 | 0.18 | 0.42
PC | 0.26 | 0.32 | 0.6 | 0.42 | 47 | 0.101 | 0.18 | 0.58
NOTEARS ($\lambda=0.005$) | 0.19 | 0.37 | 0.28 | 0.32 | 19 | 0.093 | 0.14 | 0.68
DAGMA ($\lambda=0$) | 0.22 | 0.4 | 0.32 | 0.36 | 20 | 0.08 | 0.12 | 0.64
LLM Methods |  |  |  |  |  |  |  |
Pairwise | 0.13 | 0.16 | 0.4 | 0.23 | 61 | 0.165 | 0.215 | 0.76
Pairwise + 1000 samples | 0.16 | 0.19 | 0.52 | 0.28 | 25 | 0.17 | 0.24 | 0.72
Pairwise + 10000 samples | 0.14 | 0.16 | 0.56 | 0.25 | 86 | 0.21 | 0.28 | 0.75
Ours | 0.3 | 0.44 | 0.48 | 0.46 | 27 | 0.07 | 0.13 | 0.54
Ours + 1000 samples | 0.4 | 0.58 | 0.56 | 0.57 | 24 | 0.065 | 0.15 | 0.43
Ours + 10000 samples | 0.46 | 0.59 | 0.68 | 0.63 | 25 | 0.055 | 0.15 | 0.37
Table 2: Results on the Child causal graph (20 nodes, 25 edges). The proposed
method with observational statistics from 10000 samples outperforms all other
methods with a F score of 0.63 and a NHD ratio of 0.37.
---
Method | Acc. $\uparrow$ | Prec. $\uparrow$ | Recall $\uparrow$ | F Score $\uparrow$ | NPE | NHD $\downarrow$ | Ref. NHD | NHD Ratio $\downarrow$
Numerical Methods |  |  |  |  |  |  |  |
100 Samples |  |  |  |  |  |  |  |
GES | N/A | N/A | N/A | N/A | N/A | N/A | N/A | N/A
PC | 0.02 | 0.08 | 0.03 | 0.04 | 311 | 0.022 | 0.023 | 0.96
NOTEARS ($\lambda=0.01$) | 0.029 | 0.32 | 0.03 | 0.06 | 76 | 0.1 | 0.11 | 0.94
DAGMA ($\lambda=0.01$) | 0.03 | 0.3 | 0.03 | 0.05 | 77 | 0.10 | 0.11 | 0.95
1000 Samples |  |  |  |  |  |  |  |
GES | N/A | N/A | N/A | N/A | N/A | N/A | N/A | N/A
PC | 0.04 | 0.09 | 0.07 | 0.08 | 559 | 0.025 | 0.027 | 0.92
NOTEARS ($\lambda=0.01$) | 0.02 | 0.47 | 0.02 | 0.04 | 36 | 0.34 | 0.35 | 0.96
DAGMA ($\lambda=0$) | 0.03 | 0.098 | 0.040 | 0.057 | 316 | 0.040 | 0.042 | 0.94
10000 Samples |  |  |  |  |  |  |  |
GES | N/A | N/A | N/A | N/A | N/A | N/A | N/A | N/A
PC | N/A | N/A | N/A | N/A | N/A | N/A | N/A | N/A
NOTEARS ($\lambda=0.005$) | 0.03 | 0.45 | 0.031 | 0.058 | 53 | 0.19 | 0.20 | 0.94
DAGMA ($\lambda=0$) | 0.033 | 0.14 | 0.040 | 0.063 | 214 | 0.059 | 0.063 | 0.94
LLM Methods |  |  |  |  |  |  |  |
Pairwise | N/A | N/A | N/A | N/A | N/A | N/A | N/A | N/A
Ours | 0.217 | 0.583 | 0.251 | 0.351 | 331 | 0.014 | 0.022 | 0.643
Table 3: Results on the Neuropathic Pain causal graph (221 nodes, 770 edges).
We omit the results for GES and pairwise queries because they are intractable
to use on a graph of this size. All methods except the proposed method fail
catastrophically. The proposed method performs reasonably well with a F score
of 0.351 and a NHD ratio of 0.643.

###  4.3 Results

Tables 1, 2 and 3 show the results of our experiments on the Asia causal
graph, the Child causal graph, and the Neuropathic Pain causal graph
respectively. We bold the best result in the F Score column and NHD Ratio
column since those are the main metrics we are looking at for comparison.

On the Asia causal graph, the proposed method without observational statistics
performs the best of all the methods by far, achieving an F score of 0.93 and
a NHD ratio of 0.067. Of the numerical methods, in the low data regime,
NOTEARS performs the best, while GES performs the best in the medium and high
data regimes. Interestingly, adding observational statistics to the prompt
does not improve the performance of the proposed method for this graph,
potentially because it is small enough for the LLM to be able to use its
common sense knowledge to discover the causal relationships between all of its
variables, so the addition of observational statistics is unhelpful and only
adds unnecessary noise to the prompt.

On the Child causal graph, the proposed method with observational statistics
from 10000 samples performs the best out of all the methods, with an F score
of 0.63 and a NHD Ratio of 0.37. The proposed method with observational
statistics from 1000 samples performs almost exactly as well as the best
baseline (GES with 10000 samples) while using 10x less data. Finally, the
proposed method without observational statistics is beaten only by GES with
10000 samples and performs as well as the next best baseline (GES with 1000
samples) while using no observational data. These results also showcase how
incorporating observational statistics into the LLM’s prompt can significantly
improve performance.

On the Neuropathic Pain causal graph, running pairwise queries would require
$221\choose 2$ $=24310$ queries, which is impractical. Similarly, GES requires
considering all 24310 possible edges at each step of its search, which is
intractable. We thus omit the results for both GES and pairwise queries on the
Neuropathic Pain graph. In addition, PC runs out of memory when run on the
Neuropathic Pain graph with 10000 samples, so we also omit those results from
the table.

On this much larger graph, we see that even the methods which are able to run
fail catastrophically, barely outperforming the baseline which predicts the
same number of edges but all wrong (as seen by the NHD Ratios which are all
above 0.9). In contrast, the proposed method achieves reasonable performance
with an F Score of 0.351 and a NHD Ratio of 0.643, thus emerging as the only
method which is able to handle larger graphs such as this one.

Since the proposed method without observational data already outperforms the
other methods significantly on the Neuropathic Pain causal graph and the
number of queries required for this large graph is quite significant even
considering the efficiency of The proposed method, we choose not to run the
proposed method with observational statistics on this graph due to budget
constraints.

##  5 Limitations

The proposed method has the same limitations as any method that relies on
LLMs.

First, since LLMs rely on the knowledge from their training data to synthesize
the causal relationships, the proposed approach only works on real-world data.
Second, due to their black-box nature, the generations made by the LLM lack
attribution or grounding, limiting transparency. Finally, the performance of
the proposed method is dependent on the quality of the LLM used.

##  6 Conclusion

In this work, we introduce a novel method for causal graph discovery which
leverages a breadth-first search approach in conjunction with Large Language
Models. In contrast to statistical approaches, the proposed method eliminates
the need for observational data and extra training time, ensuring quick and
easy applicability across diverse problem domains. When observational data is
available, we show that it can be incorporated into the proposed method to
further improve performance. Compared to existing LLM-based methods, the
proposed approach achieves superior efficiency by scaling the number of
queries linearly with the number of variables, rather than quadratically.
Furthermore, the proposed method achieves state-of-the-art performance on
three causal graphs of varying sizes, and is the only method tested to achieve
reasonable performance on the large Neuropathic Pain causal graph. The
proposed method thus emerges as a more efficient, flexible, and effective
solution to causal graph discovery.

##  7 Future Work

One avenue that future work should aim to explore is further combining the
common sense knowledge of LLMs with information from observational data by
employing hybrid statistical/LLM approaches. Potential strategies include
using the LLM-generated causal graph as a prior for statistical methods
(similar to Choi et al. (2022)), employing statistical methods to validate
causal relationships suggested by the LLM against observational data, and
incorporating different observational statistics into the LLM prompt.

Another avenue to explore is using advanced prompting strategies like Tree of
Thoughts (Yao et al., 2023), which are known to significantly improve
performance in complex multi-step reasoning tasks like causal graph discovery.

A final avenue to explore is trying the proposed method with different LLMs,
in order to provide insights into how the size and quality of the LLM
influences the proposed method’s performance.

##  8 Impact Statement

Causal graph discovery has many potential applications in medicine
(understanding Alzheimer’s pathophysiology), biology (constructing gene
regulatory networks), and psychology (understanding the causal relationships
between psychiatric disorders) (Zanga & Stella, 2023) which, if successful,
could have far-reaching positive impacts on society. We hope that our work
takes a step towards achieving those positive impacts by advancing the field
of causal graph discovery.

## References

  * Bello et al. (2023) Bello, K., Aragam, B., and Ravikumar, P.  Dagma: Learning dags via m-matrices and a log-determinant acyclicity characterization, 2023.
  * Choi et al. (2022) Choi, K., Cundy, C., Srivastava, S., and Ermon, S.  Lmpriors: Pre-trained language models as task-specific priors.  _arXiv preprint arXiv: 2210.12530_ , 2022.
  * Gendron et al. (2023) Gendron, G., Bao, Q., Witbrock, M., and Dobbie, G.  Large language models are not strong abstract reasoners.  _arXiv preprint arXiv: 2305.19555_ , 2023.
  * Imani et al. (2023) Imani, S., Du, L., and Shrivastava, H.  Mathprompter: Mathematical reasoning using large language models.  _Annual Meeting of the Association for Computational Linguistics_ , 2023.  doi: 10.48550/arXiv.2303.05398.
  * Kıcıman et al. (2023) Kıcıman, E., Ness, R., Sharma, A., and Tan, C.  Causal reasoning and large language models: Opening a new frontier for causality.  _arXiv preprint arXiv: 2305.00050_ , 2023.
  * Lauritzen & Spiegelhalter (1988) Lauritzen, S. L. and Spiegelhalter, D. J.  Local computations with probabilities on graphical structures and their application to expert systems.  _Journal of the Royal Statistical Society. Series B (Methodological)_ , 50(2):157–224, 1988.  ISSN 00359246.  URL http://www.jstor.org/stable/2345762.
  * Li et al. (2023) Li, Y., Xu, M., Miao, X., Zhou, S., and Qian, T.  Large language models as counterfactual generator: Strengths and weaknesses.  _arXiv preprint arXiv: 2305.14791_ , 2023.
  * Long et al. (2023a) Long, S., Piché, A., Zantedeschi, V., Schuster, T., and Drouin, A.  Causal discovery with language models as imperfect experts, 2023a.
  * Long et al. (2023b) Long, S., Schuster, T., Piché, A., of Family Medicine, D., University, M., Mila, de Montreal, U., and Research, S.  Can large language models build causal graphs?  _arXiv preprint arXiv: 2303.05279_ , 2023b.
  * Meek (2023) Meek, C.  Graphical Models: Selecting causal and statistical models.  4 2023.  doi: 10.1184/R1/22696393.v1.  URL https://kilthub.cmu.edu/articles/thesis/Graphical_Models_Selecting_causal_and_statistical_models/22696393.
  * OpenAI (2023) OpenAI.  Gpt-4 technical report.  _arXiv preprint arXiv: 2303.08774_ , 2023.
  * Pearl (2009) Pearl, J.  _Causality_.  Cambridge University Press, Cambridge, UK, 2 edition, 2009.  ISBN 978-0-521-89560-6.  doi: 10.1017/CBO9780511803161.
  * Peters et al. (2017) Peters, J., Janzing, D., and Schlkopf, B.  _Elements of Causal Inference: Foundations and Learning Algorithms_.  The MIT Press, 2017.  ISBN 0262037319.
  * Rozière et al. (2023) Rozière, B., Gehring, J., Gloeckle, F., Sootla, S., Gat, I., Tan, X. E., Adi, Y., Liu, J., Remez, T., Rapin, J., Kozhevnikov, A., Evtimov, I., Bitton, J., Bhatt, M., Ferrer, C. C., Grattafiori, A., Xiong, W., Défossez, A., Copet, J., Azhar, F., Touvron, H., Martin, L., Usunier, N., Scialom, T., and Synnaeve, G.  Code llama: Open foundation models for code.  _arXiv preprint arXiv: 2308.12950_ , 2023.
  * Schwarz (1978) Schwarz, G.  Estimating the dimension of a model.  _The Annals of Statistics_ , 6(2):461–464, 1978.  ISSN 00905364.  URL http://www.jstor.org/stable/2958889.
  * Scutari (2010) Scutari, M.  Learning bayesian networks with the bnlearn R package.  _Journal of Statistical Software_ , 35(3):1–22, 2010.  doi: 10.18637/jss.v035.i03.
  * Spiegelhalter et al. (1993) Spiegelhalter, D. J., Dawid, A. P., Lauritzen, S. L., and Cowell, R. G.  Bayesian analysis in expert systems.  _Statistical Science_ , 8(3):219–247, 1993.  ISSN 08834237.  URL http://www.jstor.org/stable/2245959.
  * Spirtes & Glymour (1991) Spirtes, P. and Glymour, C.  An algorithm for fast recovery of sparse causal graphs.  _Social Science Computer Review_ , 9(1):62–72, 1991.  doi: 10.1177/089443939100900106.  URL https://doi.org/10.1177/089443939100900106.
  * Spirtes et al. (1993) Spirtes, P., Glymour, C., and Scheines, R.  _Causation, Prediction, and Search_ , volume 81.  01 1993.  ISBN 978-1-4612-7650-0.  doi: 10.1007/978-1-4612-2748-9.
  * Tu et al. (2019) Tu, R., Zhang, K., Bertilson, B., Kjellström, H., and Zhang, C.  Neuropathic pain diagnosis simulator for causal discovery algorithm evaluation.  _Neural Information Processing Systems_ , 2019.
  * Wang et al. (2023) Wang, P., Chan, A., Ilievski, F., Chen, M., and Ren, X.  PINTO: Faithful language reasoning using prompt-generated rationales.  In _The Eleventh International Conference on Learning Representations_ , 2023.  URL https://openreview.net/forum?id=WBXbRs63oVu.
  * Wei et al. (2022) Wei, J., Wang, X., Schuurmans, D., Bosma, M., Ichter, B., Xia, F., Chi, E., Le, Q., and Zhou, D.  Chain-of-thought prompting elicits reasoning in large language models.  _arXiv preprint arXiv: 2201.11903_ , 2022.
  * Xu et al. (2023) Xu, Y., Li, W., Vaezipoor, P., Sanner, S., and Khalil, E. B.  Llms and the abstraction and reasoning corpus: Successes, failures, and the importance of object-based representations.  _arXiv preprint arXiv: 2305.18354_ , 2023.
  * Yao et al. (2023) Yao, S., Yu, D., Zhao, J., Shafran, I., Griffiths, T. L., Cao, Y., and Narasimhan, K.  Tree of thoughts: Deliberate problem solving with large language models, 2023.
  * Zanga & Stella (2023) Zanga, A. and Stella, F.  A survey on causal discovery: Theory and practice, 2023.
  * Zhao et al. (2023) Zhao, Z., Lee, W. S., and Hsu, D.  Large language models as commonsense knowledge for large-scale task planning.  _arXiv preprint arXiv: 2305.14078_ , 2023.
  * Zheng et al. (2018) Zheng, X., Aragam, B., Ravikumar, P., and Xing, E. P.  Dags with no tears: Continuous optimization for structure learning, 2018.
  * Zheng et al. (2023) Zheng, Y., Huang, B., Chen, W., Ramsey, J., Gong, M., Cai, R., Shimizu, S., Spirtes, P., and Zhang, K.  Causal-learn: Causal discovery in python.  _arXiv preprint arXiv:2307.16405_ , 2023.