The product and chain rules of differential calculus can be generalized to arbitrarily long derivatives using the general Leibniz rule and Faà di Bruno’s formula, respectively, which can each be proven with a little bit of combinatorics. These generalizations themselves can be generalized to differentiating arbitrarily long products/compositions of functions, with the associative nature of multiplication and composition allowing for an efficient algorithm for computation. Along our journey with the chain rule we enter the “BellY” of the beast of Wolfram’s functionality for computing Bell polynomials and their various generalizations.

Generalizing the Product Rule

The product rule
D
t
f(t)g(t)=
′
f
(t)g(t)+f(t)
′
g
(t)
gives the first derivative of the product of two functions:
In[]:=
D[f[t]*g[t],t]
Out[]=
g[t]
′
f
[t]+f[t]
′
g
[t]

How do we generalize to find the
n
th
derivative of a product of two functions
f
and
g
?

In[]:=
Table[D[f[t]*g[t],{t,i}],{i,0,4}]//TableForm
Out[]//TableForm=
f[t]g[t]
g[t]
′
f
[t]+f[t]
′
g
[t]
2
′
f
[t]
′
g
[t]+g[t]
′′
f
[t]+f[t]
′′
g
[t]
3
′
g
[t]
′′
f
[t]+3
′
f
[t]
′′
g
[t]+g[t]
(3)
f
[t]+f[t]
(3)
g
[t]
6
′′
f
[t]
′′
g
[t]+4
′
g
[t]
(3)
f
[t]+4
′
f
[t]
(3)
g
[t]+g[t]
(4)
f
[t]+f[t]
(4)
g
[t]
Starting with the product
f(t)g(t)=
(0)
f
(t)
(0)
g
(t)
of two functions of
t
, each time we take a derivative with respect to
t
, we are left with a sum of terms, each of the form
c
(i)
f
(t)
(j)
g
(t)
.
By the product rule and the linearity of differentiation, differentiating such a sum once more splits each such term into the expression
c
(i+1)
f
(t)
(j)
g
(t)+c
(i)
f
(t)
(j+1)
g
(t)
.
If we replace
(i)
f
(t)
with the expression
i
x
and
(j)
g
(t)
with the expression
j
y
, then we see that differentiating
c
(i)
f
(t)
(j)
g
(t)
to get
c
(i+1)
f
(t)
(j)
g
(t)+c
(i)
f
(t)
(j+1)
g
(t)
is analogous to multiplying
c
i
x
j
y
by
x+y
to get
c
i+1
x
j
y
+c
i
x
j+1
y
.
Note that
n
D
t
f(t)g(t)
comes from starting with
f(t)g(t)=
(0)
f
(t)
(0)
g
(t)
and differentiating a total of
n
times, and
n
(x+y)
comes from starting with
1=
0
x
0
y
and multiplying by
x+y
a total of
n
times.
The expansion of
n
(x+y)
comes from summing together all terms resulting from picking either
x
or
y
from each of the
n
factors of
x+y
, so after grouping together terms with the same
x
exponent and same
y
exponent, we see that each addend must be of the form
c
n-k
x
k
y
for the constant
c
counting the ways to pick
y
exactly
k
times out of
n
chances (picking
x
the other
n-k
times by default). Hence
c=
n
k

and thus we have
n
(x+y)
=
n
∑
k=0

n
k

n-k
x
k
y
,
known as the binomial theorem.
Substituting our functions
f(t)
and
g(t)
back into the binomial theorem yields
n
D
t
f(t)g(t)=
n
∑
k=0

n
k

(n-k)
f
(t)
(k)
g
(t)
,
known as the general Leibniz rule:
In[]:=
D[f[t]*g[t],{t,n}]
Out[]=
n
∑
K[1]=0
Binomial[n,K[1]]
∂
{t,n-K[1]}
g[t]
(K[1])
f
[t]

How do we further generalize to find the
n
th
derivative of a product of
m
functions?

In[]:=
Table[D[f[t]*g[t]*h[t],{t,i}],{i,0,4}]//TableForm
Out[]//TableForm=
f[t]g[t]h[t]
g[t]h[t]
′
f
[t]+f[t]h[t]
′
g
[t]+f[t]g[t]
′
h
[t]
2(g[t]
′
f
[t]+f[t]
′
g
[t])
′
h
[t]+h[t](2
′
f
[t]
′
g
[t]+g[t]
′′
f
[t]+f[t]
′′
g
[t])+f[t]g[t]
′′
h
[t]
3
′
h
[t](2
′
f
[t]
′
g
[t]+g[t]
′′
f
[t]+f[t]
′′
g
[t])+3(g[t]
′
f
[t]+f[t]
′
g
[t])
′′
h
[t]+h[t]3
′
g
[t]
′′
f
[t]+3
′
f
[t]
′′
g
[t]+g[t]
(3)
f
[t]+f[t]
(3)
g
[t]+f[t]g[t]
(3)
h
[t]
6(2
′
f
[t]
′
g
[t]+g[t]
′′
f
[t]+f[t]
′′
g
[t])
′′
h
[t]+4
′
h
[t]3
′
g
[t]
′′
f
[t]+3
′
f
[t]
′′
g
[t]+g[t]
(3)
f
[t]+f[t]
(3)
g
[t]+4(g[t]
′
f
[t]+f[t]
′
g
[t])
(3)
h
[t]+h[t]6
′′
f
[t]
′′
g
[t]+4
′
g
[t]
(3)
f
[t]+4
′
f
[t]
(3)
g
[t]+g[t]
(4)
f
[t]+f[t]
(4)
g
[t]+f[t]g[t]
(4)
h
[t]
The general Leibniz rule can be generalized with the formula
(n)
(
f
1
⋯
f
m
)
=
∑
k
1
+⋯+
k
m
=n
Multinomial[
k
1
,…,
k
m
]
m
∏
t=1
(
k
t
)
f
t
which relates to the multinomial theorem
n
(
x
1
+⋯+
x
m
)
=
∑
k
1
+⋯+
k
m
=n
Multinomial[
k
1
,…,
k
m
]
m
∏
t=1
k
t
x
t
the same way the general Leibniz rule relates to the binomial theorem:
In[]:=
D[f[t]*g[t]*h[t],{t,n}]
Out[]=
n
∑
K[1]=0
n
∑
K[2]=0
∂
{t,n-K[1]-K[2]}
h[t]Multinomial[K[1],n-K[1]-K[2],K[2]]
(K[1])
f
[t]
(K[2])
g
[t]

Generalizing the Chain Rule

The chain rule
D
t
f(g(t))=
′
f
(g(t))
′
g
(t)
gives the first derivative of the composition of two functions:
In[]:=
D[f[g[t]],t]
Out[]=
′
f
[g[t]]
′
g
[t]

How do we generalize to find the
n
th
derivative of a composition of two functions
f
and
g
?

In[]:=
Table[D[f[g[t]],{t,i}],{i,0,4}]//TableForm
Out[]//TableForm=
f[g[t]]
′
f
[g[t]]
′
g
[t]
2
′
g
[t]
′′
f
[g[t]]+
′
f
[g[t]]
′′
g
[t]
3
′
g
[t]
′′
f
[g[t]]
′′
g
[t]+
3
′
g
[t]
(3)
f
[g[t]]+
′
f
[g[t]]
(3)
g
[t]
3
′′
f
[g[t]]
2
′′
g
[t]
+6
2
′
g
[t]
′′
g
[t]
(3)
f
[g[t]]+4
′
g
[t]
′′
f
[g[t]]
(3)
g
[t]+
4
′
g
[t]
(4)
f
[g[t]]+
′
f
[g[t]]
(4)
g
[t]
Faà di Bruno’s formula generalizes the chain rule by expressing the
n
th
derivative of a composition of two functions
n
D
t
f(g(t))=
n
∑
k=1
Y
n,k

′
g
(t),
′′
g
(t),…,
(n-k+1)
g
(t)
(k)
f
(g(t))
where
Y
n,k
(
x
1
,...,
x
-k+n+1
)
is the partial Bell polynomial defined by
Y
n,k
(
x
1
,…,
x
n-k+1
)=
n
∑
m
1
=0
⋯
n
∑
m
n
=0
n!
m
1
!⋯
m
n
!
Boole[
m
1
+2
m
2
+⋯+n
m
n
n∧
m
1
+
m
2
+⋯+
m
n
k]
n
∏
s=1
m
s
(
x
s
/s!)
The partial Bell polynomial
Y
n,k
(
x
1
,...,
x
-k+n+1
)
encodes how a set of
n
objects can be partitioned into
k
nonempty subsets
by summing together terms of the form
c
j
1
x
i
1
...
j
m
x
i
m
, which indicates there are
c
ways to partition the set of
n
objects into
j
1
subsets each of size
i
1
,
j
2
subsets each of size
i
2
, etc. (so
i
1
j
1
+...+
i
m
j
m
=n
and
j
1
+...+
j
m
=k
necessarily). For instance…
In[]:=
BellY[6,2,{
x
1
,
x
2
,
x
3
,
x
4
,
x
5
}]
Out[]=
10
2
x
3
+15
x
2
x
4
+6
x
1
x
5
…indicates that when partitioning a set of six objects into two nonempty subsets, there are 10 partitions consisting of 2 subsets each of size 3, 15 partitions consisting of 1 subset of size 2 and 1 subset of size 4, and 6 partitions consisting of 1 subset of size 1 and 1 subset of size 5.
To see why Faà di Bruno’s formula works, note that
n
D
t
f(g(t))
is the sum of like terms of the form…
c*
(k)
f
(g(t))*
j
1

(1)
g
(t)
*...*
j
n

(n)
g
(t)
…so it remains to prove that the set of partitions of
{1,2,...,n}
into
k
nonempty subsets, consisting of
j
p
subsets each of size
p
for
1<=p<=n
(so
1*
j
1
+...+n*
j
n
=n
and
j
1
+...+
j
n
=k
necessarily), has cardinality
c
. See that the case
n=1
holds because
1
D
t
f(g(t))=f'(g(t))g'(t)
and the only partition of
{1}
is
{{1}}
. Set
n>=1
and assume inductively the hypothesis holds, and consider the first derivative of the above term
c*
(k)
f
(g(t))*
j
1

(1)
g
(t)
*...*
j
n

(n)
g
(t)
:
D
t
c*
(k)
f
(g(t))*
j
1

(1)
g
(t)
*...*
j
n

(n)
g
(t)
​​=c*
(k+1)
f
(g(t))*
j
1
+1

(1)
g
(t)
*...*
j
n

(n)
g
(t)
​​+c*
j
1
*
(k)
f
(g(t))*
j
1
-1

(1)
g
(t)
*
j
1
+1

(2)
g
(t)
*...*
j
n

(n)
g
(t)
​​+...​​+c*
j
n
*
(k)
f
(g(t))*
j
1

(1)
g
(t)
*...*
j
n
-1

(n)
g
(t)
*
1

(n+1)
g
(t)
The new object
n+1
must be added to the existing partitions of
{1,2,...,n}
corresponding to the term
c*
(k)
f
(g(t))*
j
1

(1)
g
(t)
*...*
j
n

(n)
g
(t)
in order to create a partition of
{1,2,...,n+1}
. It can either be added as its own subset
{n+1}
to each partition, in which case each new partition has
k+1
parts and one additional part of size 1, corresponding to the new term…
c*
(k+1)
f
(g(t))*
j
1
+1

(1)
g
(t)
*...*
j
n

(n)
g
(t)
…or it can be added to an existing subset of size
p
, in which case each new partition still has
k
parts, with one fewer of size
p
and one more of size
p+1
, corresponding to the new term…
c*
j
p
*
(k)
f
(g(t))*
j
1

(1)
g
(t)
*...*
j
p
-1

(p)
g
(t)
*
j
p
+1

(p+1)
g
(t)
*...*
j
n

(n)
g
(t)
…completing the proof.
Table of
Y
n,k
(
x
1
,...,
x
-k+n+1
)
values:
In[]:=
Table[BellY[n,k,Array[
x
#
&,n-k+1]],{n,1,7},{k,1,n}]//TableForm
Out[]//TableForm=
x
1
x
2
2
x
1
x
3
3
x
1
x
2
3
x
1
x
4
3
2
x
2
+4
x
1
x
3
6
2
x
1
x
2
4
x
1
x
5
10
x
2
x
3
+5
x
1
x
4
15
x
1
2
x
2
+10
2
x
1
x
3
10
3
x
1
x
2
5
x
1
x
6
10
2
x
3
+15
x
2
x
4
+6
x
1
x
5
15
3
x
2
+60
x
1
x
2
x
3
+15
2
x
1
x
4
45
2
x
1
2
x
2
+20
3
x
1
x
3
15
4
x
1
x
2
6
x
1
x
7
35
x
3
x
4
+21
x
2
x
5
+7
x
1
x
6
105
2
x
2
x
3
+70
x
1
2
x
3
+105
x
1
x
2
x
4
+21
2
x
1
x
5
105
x
1
3
x
2
+210
2
x
1
x
2
x
3
+35
3
x
1
x
4
105
3
x
1
2
x
2
+35
4
x
1
x
3
21
5
x
1
x
2
7
x
1
Y
n,k
(
x
1
,...,
x
-k+n+1
)
appears in abstract computations:
In[]:=
D[f[g[t]],{t,n}]
Out[]=
n
∑
K[1]=0
BellYn,K[1],Table
(K[2])
g
[t],{K[2],1,1+n-K[1]}
(K[1])
f
[g[t]]
Compare to:
Applying Faà di Bruno’s formula twice for compositions of three functions is a bit unwieldy…
…which we see is identical.

How do we modify Faà di Bruno’s formula to fit this generalization?

…so we confirm our original statement by making the following computations…

How do we define generalized partial Bell polynomials explicitly?

First way: As a special case of a generalized Bell polynomial

…which is a formula explicitly expressing the generalized partial Bell polynomials in terms of generalized Bell polynomials.
Indeed, compare…
…to the following…
…and see that they are identical.

Second way: As a special case of a Bell polynomial

…and replace it with the left side to get…
…which by Faà di Bruno’s formula equals…
…which means we have…
…and thus…
…which is an explicit formula for the generalized partial Bell polynomials in terms of the partial Bell polynomials and the generalized Bell polynomials.
Indeed, see that…
…is identical to the table of generalized partial Bell polynomials computed previously.

Recursive Algorithm for Computing the Generalized Bell Polynomial and Its Product Rule Analog

Code

Implementation

Implement Faà di Bruno's formula for a matrix having exactly two columns:
Implement the general Leibniz rule for a matrix having exactly two columns:
Apply a two-column rule to pairs of columns of the given matrix to (roughly) halve the number of columns:
Use a two-column rule repeatedly on a general matrix:

Testing

See that the output from our recursive algorithm matches that of the built-in functionality for generalized Bell polynomials:
See that the two functions return the same output, but the recursive one is significantly faster:

CITE THIS NOTEBOOK

Combinatorics, calculus, Bell polynomials, and you​
by Dan McDonald​
Wolfram Community, STAFF PICKS, September 4, 2025
​https://community.wolfram.com/groups/-/m/t/3540248