###### Abstract

We provide a novel representation of the total n-th derivative of the
multivariate composite function $f\circ g$, i.e. a generalized Faà di Bruno’s
formula. To this end, we make use of properties of the Kronecker product and
the n-th derivative of the left-composite $f$, which allow the use of a
multivariate form of partial Bell polynomials to represent the generalized Faà
di Bruno’s formula. We further show that standard recurrence relations that
hold for the univariate partial Bell polynomial also hold for the multivariate
partial Bell polynomial under a simple transformation. We apply this
generalization of Faà di Bruno’s formula to the computation of multivariate
moments of the normal distribution.

##  1 Introduction

In this paper we set out to construct an easy to handle representation of
higher total derivatives of composite functions which might be vector-valued
and multivariate. We generalize the well-known Faà Di Bruno’s formula in a way
which is identical in structure to the univariate case. Specifically, we
provide a parsimonious yet valuable representation for the n-th total
derivative using a multivariate generalization of partial Bell polynomials:

$d^{n}f(g(x))=\left(\sum_{k=1}^{n}f_{g^{k}}\textbf{B}^{g}_{n,k}\right)(dx)^{\otimes
n},$ (1.1)

where $\textbf{B}^{g}_{n,k}$ is identical in structure to the univariate
partial Bell polynomial over the respective matrices
$\left(g_{x},g_{x^{2}},\ldots,g_{x^{n-k+1}}\right)$ containing all partial
derivatives of $g$ w.r.t $x$ at all orders $l=1,2,...,n-k+1$ of sizes
$n_{y}\times n_{x}^{l}$, with scalar multiplication substituted by the
Kronecker product $\otimes$ and with $f_{g^{k}}$ being a $n_{f}\times
n_{y}^{k}$ matrix collecting all k-th order partial derivatives of $f$ w.r.t
$g$.111Representing all k-th order partial derivatives in one vector or matrix
instead of using multi-indices or tensors is crucial for this representation.
For example, such a representation has been advanced in Holmquist (1996) and
Magnus (2010).

Faà Di Bruno’s formula for the univariate case is a well-known result. An
overview of its history can be found in Johnson (2002). Multiple papers
generalize Faà Di Bruno’s formula for vector-valued and multivariate composite
$f\circ g$ to obtain higher order partial derivatives. Examples for
multivariate cases can be found in Gzyl (1986), Constantine and Savits (1996),
Hardy (2006) and Leipnik and Pearce (2007). However, they do not employ Bell
polynomials to express the combinatorial structure of the higher order
composite derivatives. Noschese and Ricci (2003) and Schumann (2019) make use
of Bell polynomials to derive higher order partial derivatives and Di Nardo,
Guarino, and Senato (2011) applies classical umbra calculus. A generalization
of Faà Di Bruno’s formula to obtain rather the n-th total derivative are
provided by Mishkov (2000) and for the matrix of n-th order derivatives by
Chacón and Duong (2021). However, both representations to do not exploit the
powerful properties of Bell polynomials.

In this paper, we offer a generalized representation of the n-th total
derivative for a vector-valued and multivariate composite $f\circ g$ that
takes advantage of both the properties of Bell polynomials as well as the
Kronecker product. Firstly, we thereby systematically collect all higher order
partial derivatives in a novel multivariate Bell Polynomial and make it thus
directly applicable to compute e.g. Taylor series expansions. Secondly, the
novel multivariate partial Bell polynomials inherit the common recurrence
relationships familiar from the univariate case which leads us directly to
apply them to a multivariate Faà Di Bruno’s formula. Thirdly, we apply this
novel representation to the computation of arbitrary moments of the
multivariate normal distribution.

We start by recalling the univariate partial Bell polynomial and present a
proof of the univariate version of Faà di Bruno’s formula using induction in
Section 2. This section serves as an expository and pedagogical guideline,
since these results are already well-known in the univariate case, but our
proof for the multivariate version works along similar lines. In Section 3, we
introduce the notation for matrix-valued derivatives and present our
multivariate version of the partial Bell polynomial. In this section, we also
lay out the issues when attempting to use the multivariate partial Bell
polynomial and how these can be resolved. In Section 4, we show, on the basis
of commutation properties of the Kronecker product, that the familiar
recurrence relations for partial Bell polynomials also hold in the
multivariate case. We finally prove the generalized version of Faà Di Bruno’s
formula in Section 5. In Section 6, we apply the generalization to the
computation of moments of the multivariate normal distribution. We add some
concluding remarks in Section 7.

Notation and Definitions: Where it is not defined differently, functions
$f(y)$ and $g(x)$ are vectors of size $n_{f}\times 1$ and $n_{y}\times 1$,
respectively, given by $y=g\left(x\right)$ with
$f:\mathbb{R}^{n_{y}}\rightarrow\mathbb{R}^{n_{f}}$ and
$g:\mathbb{R}^{n_{x}}\rightarrow\mathbb{R}^{n_{y}}$222Obviously $f,g$ can be
defined on more narrow domains as long as their are sufficiently often
differentiable on these.. Both $f$ and $g$ are at least n times
differentiable. Further let $x=\left(x_{1},\ldots,x_{n_{x}}\right)^{\prime}$
and $y=\left(y_{1},\ldots,y_{n_{y}}\right)^{\prime}$ be column vectors.
$\mathcal{M}_{r\times c}\left(\mathbb{R}\right)$ represents the set of all
matrices over the real numbers with $r$ rows and $c$ columns. $I_{n}$ denotes
the identity matrix of size $n\times n$. $Sym(m)$ is the symmetric group over
the set of elements given by $\left\\{1,\ldots m\right\\}$ and $\sigma\in
Sym(m)$ represents a permutation defined over these $m$ elements.

##  2 Univariate Faà Di Bruno

To arrive at some intuition and guidance for the proof of the multivariate Faà
Di Bruno’s formula, we start by providing a proof of the univariate case. We
define the partial Bell polynomial333For reference see for example Comtet
(1974). over the function $g$ as follows

######  Definition 1 (Univariate Partial Bell Polynomial).

For some $n,k\in\mathbb{N}$ with $n\geq k$ and a function
$g:\mathbb{R}\rightarrow\mathbb{R}$ which is at least $n-k+1$ times
differentiable on some domain $D\subset\mathbb{R}$, we can define the partial
Bell polynomial by

$\displaystyle\begin{array}[]{l}B_{n,k}(g_{x}(x),g_{x^{2}}(x),\cdots,g_{x^{(n-k+1)}}(x))\\\\[6.45831pt]
=\sum{n!\over j_{1}!j_{2}!\cdots j_{n-k+1}!}\left({g_{x}(x)\over
1!}\right)^{j_{1}}\left({g_{x^{2}}(x)\over
2!}\right)^{j_{2}}\cdots\left({g_{x^{(n-k+1)}}(x)\over(n-k+1)!}\right)^{j_{n-k+1}},\end{array}$
(2.4)

where $\left(g_{x},g_{x^{2}},\ldots,g_{x^{n-k+1}}\right)$ denote the
univariate partial derivatives of $g$ w.r.t $x$ at all orders
$i=1,2,...,n-k+1$, and where the sum is taken over all sequences $j_{1}$,
$j_{2}$, $j_{3}$,…., $j_{n-k+1}$ of non-negative integers such that these two
conditions are satisfied:

$\displaystyle\begin{array}[]{cc}j_{1}+j_{2}+\cdots+j_{n-k+1}&=k\\\\[6.45831pt]
j_{1}+2j_{2}+3j_{3}+\cdots+(n-k+1)j_{n-k+1}&=n.\end{array}$ (2.7)

To ease notation, we set

$B^{g}_{n,k}:=B_{n,k}(g_{x}(x),g_{x^{2}}(x),\cdots,g_{x^{(n-k+1)}}(x)).$ (2.8)

Further it holds that $B^{g}_{n,k}=0$ for $k>n$, $n=0$, and $k=0$.

Before we proceed we require the following Lemma:

######  Lemma 1 (Differentiation Recurrence of Univariate Partial Bell
Polynomials).

For any scalars $n,k$ and functions $g$ as defined above, it holds that

$B^{g}_{n+1,k}=g_{x}(x)B^{g}_{n,k-1}+\frac{d}{dx}B^{g}_{n,k}.$ (2.9)

###### Proof.

The result follows from two well known recurrence relations for the partial
Bell polynomial444See for example Charalambides (2018).

$\displaystyle B_{n,k}^{g}$ $\displaystyle=\sum_{i=1}^{n-k+1}{n-1\choose
i-1}g_{x^{i}}B_{n-i,k-1}^{g}$ (2.10) $\displaystyle\frac{d}{dx}B_{n,k}^{g}$
$\displaystyle=\sum_{i=1}^{n-k+1}{n\choose i}g_{x^{i+1}}B_{n-i,k-1}^{g}$
(2.11)

Starting from the right-hand-side (RHS) of (2.9) and using (2.11) we can see
that

$\displaystyle\quad g_{x}B^{g}_{n,k-1}+\frac{d}{dx}B^{g}_{n,k}$ (2.12)
$\displaystyle=$ $\displaystyle\quad
g_{x}B^{g}_{n,k-1}+\sum_{i=1}^{n-k+1}{n\choose i}g_{x^{i+1}}B_{n-i,k-1}^{g}$
(2.13) $\displaystyle=$ $\displaystyle\quad
g_{x}B^{g}_{n,k-1}+\sum_{i=2}^{n-k+2}{n\choose i-1}g_{x^{i}}B_{n-i+1,k-1}^{g}$
(2.14) $\displaystyle=$ $\displaystyle\quad
g_{x}B^{g}_{n,k-1}-g_{x}B_{n,k-1}^{g}+\sum_{i=1}^{n-k+2}{n\choose
i-1}g_{x^{i}}B_{n-i+1,k-1}^{g}$ (2.15) $\displaystyle=$
$\displaystyle\quad\sum_{i=1}^{n-k+2}{n\choose
i-1}g_{x^{i}}B_{n-i+1,k-1}^{g}.$ (2.16)

By setting $n:=n+1$ in equation (2.10) and comparing to (2.16), the proof is
complete. ∎

Using Lemma 1 it is now straightforward to prove the univariate and scalar-
valued case of Faà Di Bruno by induction:

######  Theorem 1 (Univariate Faà Di Bruno’s Formula Using Bell Polynomials).

For some $n,k\in\mathbb{N}$ with $n\geq k$ and function
$g:\mathbb{R}\rightarrow\mathbb{R}$ which is at least $n-k+1$ differentiable
on some domain $D_{g}\subset\mathbb{R}$ and a function
$f:\mathbb{R}\rightarrow\mathbb{R}$ which is at least $n$ times differentiable
on some domain $D_{f}\subset\mathbb{R}$, the $n$-th order composite derivative
of $f\circ g$ is given by the Faà Di Bruno’s formula:

$\displaystyle\frac{d^{n}}{dx^{n}}f(g(x))=\sum_{k=1}^{n}f_{g^{k}}(g(x))B^{g}_{n,k}.$
(2.17)

###### Proof.

We prove the univariate case by induction. Start of induction for $n=1$ gives

$\displaystyle\frac{d}{dx}f(g(x))$
$\displaystyle=\sum_{k=1}^{1}f_{g^{k}}(g(x))B_{n,k}^{g}$ (2.18)
$\displaystyle=f_{g}\left(g(x)\right)B_{1,1}^{g}$ (2.19)
$\displaystyle=f_{g}\left(g(x)\right)g_{x}(x),$ (2.20)

which follows directly from the chain rule. Now suppose (2.17) holds for some
arbitrary $n\in\mathbb{N}$. Then we can write

$\displaystyle\frac{d^{n+1}}{dx^{n+1}}f(g(x))$
$\displaystyle=\frac{d}{dx}\left[\frac{d^{n}}{dx^{n}}f(g(x))\right]$ (2.21)
$\displaystyle=\frac{d}{dx}\left[\sum_{k=1}^{n}f_{g^{k}}(g(x))B^{g}_{n,k}\right]$
(2.22)
$\displaystyle=\sum_{k=1}^{n}f_{g^{k+1}}\left(g(x)\right)g_{x}(x)B^{g}_{n,k}+f_{g^{k}}\left(g(x)\right)\frac{d}{dx}B^{g}_{n,k}.$
(2.23)

We can factor for $f_{g^{k}}\left(g(x)\right)$ after re-parameterizing the
sums and noting that
$B^{g}_{n,0}\left(\cdot\right)=B^{g}_{n,n+1}\left(\cdot\right)=0$.

$\displaystyle=\sum_{k=1}^{n+1}f_{g^{k}}\left(g(x)\right)g_{x}(x)B^{g}_{n,k-1}+\sum_{k=1}^{n+1}f_{g^{k}}\left(g(x)\right)\frac{d}{dx}B^{g}_{n,k}$
(2.24)
$\displaystyle=\sum_{k=1}^{n+1}f_{g^{k}}\left(g(x)\right)\left(g_{x}(x)B^{g}_{n,k-1}+\frac{d}{dx}B^{g}_{n,k}\right).$
(2.25)

Because

$B^{g}_{n+1,k}=g_{x}(x)B^{g}_{n,k-1}+\frac{d}{dx}B^{g}_{n,k},$ (2.26)

as given by Lemma 1, we get

$\displaystyle\frac{d^{n+1}}{dx^{n+1}}f(g(x))=\sum_{k=1}^{n+1}f_{g^{k}}\left(g(x)\right)B^{g}_{n+1,k}.$
(2.27)

∎

##  3 Multivariate Bell Polynomials: Definition and the Issue

In this section, we define a multivariate and matrix-valued version of the
partial Bell polynomial. To that purpose, we also quickly present the notation
for vector- and matrix-valued representations of derivatives and differentials
used in this paper. Let $F$ be a matrix of size $s\times t$. We define a
differential operator ${\partial\over\partial x^{\prime}}$ w.r.t the transpose
of $x$ by

${\partial\over\partial x^{\prime}}=\left({\partial\over\partial
x_{1}},\ldots,{\partial\over\partial x_{n_{x}}}\right).$ (3.28)

First-order derivatives of matrices and vectors are then defined by

$F_{x}={\partial F\over\partial x^{\prime}}:=F\otimes{\partial\over\partial
x^{\prime}}.$ (3.29)

Analogously, higher-order derivatives are defined by

$F_{x^{k}}={\partial F\over\partial(x^{\prime})^{\otimes
k}}:=F\otimes\left({\partial\over\partial x^{\prime}}\right)^{\otimes k}.$
(3.30)

This notation implies some structure on the ordering of partial derivative
within ${\partial F\over\partial x^{\prime}}$.555See Tracy and Sultan (1993),
Sultan and Tracy (1996), and Rogers (1980) who also make use of this notation.
More explicitly this structure corresponds to the replacement of each
component $F_{i,j}$ of $F$ by its respective row vector of derivatives, i.e.

${\partial F_{i,j}\over\partial x^{\prime}}=\left({\partial
F_{i,j}\over\partial x_{1}},\ldots,{\partial F_{i,j}\over\partial
x_{d}}\right).$ (3.31)

When F is a vector, denoted by $f$ in the following, we can write the
corresponding first-order differential using

$df=f_{x}dx$ (3.32)

and the higher-order differentials by

$d^{k}f=f_{x^{k}}\left(dx\right)^{\otimes k},$ (3.33)

where the Kronecker product power notation refers to $\left(dx\right)^{\otimes
k}=\underbrace{dx\otimes\cdots\otimes dx}_{k\;\text{times}}$.

We are now equipped to define our generalization of the partial Bell
polynomial.

######  Definition 2 (Multivariate and Matrix-Valued Partial Bell Polynomial).

Let $g_{x^{l}}(x)$ denote the $n_{y}\times n_{x}^{l}$ matrix of the $l$th-
order partial derivatives of $g(\cdot)$, and where $\textbf{B}^{g}_{n,k}$ is a
$n_{y}^{k}\times n_{x}^{n}$ matrix given by

$\displaystyle\begin{array}[]{l}\textbf{B}^{g}_{n,k}:=\textbf{B}_{n,k}(g_{x}(x),g_{x^{2}}(x),\cdots,g_{x^{(n-k+1)}}(x))=\\\\[6.45831pt]
\sum{n!\over j_{1}!j_{2}!\cdots j_{n-k+1}!}\left({g_{x}(x)\over
1!}\right)^{\otimes j_{1}}\otimes\left({g_{x^{2}}(x)\over 2!}\right)^{\otimes
j_{2}}\otimes\cdots\otimes\left({g_{x^{(n-k+1)}}(x)\over(n-k+1)!}\right)^{\otimes
j_{n-k+1}},\end{array}$ (3.36)

where the sum is taken over all sequences $j_{1}$, $j_{2}$, $j_{3}$,….,
$j_{n-k+1}$ of non-negative integers such that these two conditions are
satisfied:

$\displaystyle\begin{array}[]{cc}j_{1}+j_{2}+\cdots+j_{n-k+1}&=k\\\\[6.45831pt]
j_{1}+2j_{2}+3j_{3}+\cdots+(n-k+1)j_{n-k+1}&=n.\end{array}$ (3.39)

This definition is identical in structure to the univariate case with the
exception that scalar multiplication is substituted by the Kronecker product.
The bold writing of $\textbf{B}^{g}_{n,k}$ emphasizes that this Bell
polynomial is a matrix object and not a scalar.

#### Non-Recurrence of Multivariate Partial Bell Polynomials - Example

It is important to stress that the recurrence relations of univariate Bell
Polynomials do not simply carry over to the multivariate and matrix-valued
case as the Kronecker products do not commute. To see this, suppose a
multivariate and matrix-valued recurrence relation was given by

$\textbf{B}_{n,k}^{g}=\sum_{i=1}^{n-k+1}{n-1\choose
i-1}\left(\textbf{B}_{n-i,k-1}^{g}\otimes g_{x^{i}}\right),$ (3.40)

as in the univariate case. Consider for given $g$ the following partial Bell
polynomials:

$\displaystyle\textbf{B}^{g}_{1,1}$ $\displaystyle=g_{x}$ (3.41)
$\displaystyle\textbf{B}^{g}_{2,1}$ $\displaystyle=g_{x^{2}}$ (3.42)
$\displaystyle\textbf{B}^{g}_{3,2}$ $\displaystyle=3\left(g_{x}\otimes
g_{x^{2}}\right).$ (3.43)

If the recurrence was to hold true, it would require that

$\textbf{B}^{g}_{3,2}=\textbf{B}^{g}_{2,1}\otimes
g_{x}+2\left(\textbf{B}^{g}_{1,1}\otimes g_{x^{2}}\right).$ (3.44)

Simply plugging in the given partial Bell polynomials then leads to

$3\left(g_{x}\otimes g_{x^{2}}\right)\neq g_{x^{2}}\otimes
g_{x}+2\left(g_{x}\otimes g_{x^{2}}\right).$ (3.45)

The left-hand-side does not equal the right-hand-side because the Kronecker
product does not commute.

#### Reestablishing Recurrence of Multivariate Partial Bell Polynomials -
Example cont.

However, it is well known that Kronecker products can be commuted with the
help of so-called commutation matrices. Therefore, given the existence of
commutation matrices $L$ and $R$, we can write

$g_{x^{2}}\otimes g_{x}+2\left(g_{x}\otimes
g_{x^{2}}\right)=L\left(g_{x}\otimes g_{x^{2}}\right)R+2\left(g_{x}\otimes
g_{x^{2}}\right).$ (3.46)

We will introduce commutation matrices in more detail in the next section.
Suppose now, that the function $g$ is the inner function of the composite
function $f\circ g$, and that $dx$ is a differential of the arguments of $g$.
Then it turns out, we will be able to write

$\displaystyle f_{g^{2}}\left(\textbf{B}^{g}_{2,1}\otimes
g_{x}+2\left(\textbf{B}^{g}_{1,1}\otimes g_{x^{2}}\right)\right)(dx)^{\otimes
3}$ $\displaystyle=f_{g^{2}}\left(g_{x^{2}}\otimes g_{x}+2\left(g_{x}\otimes
g_{x^{2}}\right)\right)(dx)^{\otimes 3}$ (3.47)
$\displaystyle=f_{g^{2}}\left(L\left(g_{x}\otimes
g_{x^{2}}\right)R+2\left(g_{x}\otimes g_{x^{2}}\right)\right)(dx)^{\otimes 3}$
(3.48) $\displaystyle=f_{g^{2}}L\left(g_{x}\otimes
g_{x^{2}}\right)R(dx)^{\otimes 3}+f_{g^{2}}2\left(g_{x}\otimes
g_{x^{2}}\right)(dx)^{\otimes 3}$ (3.49)
$\displaystyle=f_{g^{2}}3\left(g_{x}\otimes g_{x^{2}}\right)(dx)^{\otimes 3},$
(3.50) $\displaystyle=f_{g^{2}}\textbf{B}^{g}_{3,2}(dx)^{\otimes 3}$ (3.51)

because $f_{g^{2}}L=f_{g^{2}}$ and $R(dx)^{\otimes 3}=(dx)^{\otimes 3}$.
$f_{g^{2}}L=f_{g^{2}}$ is true because the commutation matrix $L$
coincidentally only commutes the last two colored factors in
$f_{g^{2}}=f\otimes{\color[rgb]{0,0,1}\partial\over\partial
x^{\prime}}\otimes{\color[rgb]{1,0,0}\partial\over\partial x^{\prime}}$, which
does not change the overall derivative. Similarly, the commutation matrix $R$
will commute the colored objects in ${\color[rgb]{0,0,1}dx\otimes
dx}\otimes{\color[rgb]{1,0,0}dx}$, which also doesn’t change the differential.

In the following we will show the generalization of this insight for $n,k\geq
0$ in $\textbf{B}^{g}_{n,k}$. To that purpose, we start by generalizing the
equality of (3.49) to (3.50), after which we show how we are always able to
write multivariate partial Bell polynomial recurrences as in (3.47) and
(3.51).

#### Outline

In the following we outline the steps to prove that we are able to generalize
recurrence relations for multivariate partial Bell Polynomials. Taking the
example of $\textbf{B}_{3,2}$ as a guideline, in the first step, Step 1., we
show that for all matrices $D_{k},X_{n}$ with

$\displaystyle D_{k}$ $\displaystyle:=D^{\otimes k}$ (3.52) $\displaystyle
X_{n}$ $\displaystyle:=X^{\otimes n},$ (3.53)

where $D$ is an arbitrary real-valued $1\times n_{y}$ vector, and $X$ an
arbitrary real-valued $n_{x}\times 1$ vector, we can write

$D_{k}\left(g_{x^{2}}\otimes g_{x}\right)X_{n}=D_{k}\left(g_{x}\otimes
g_{x^{2}}\right)X_{n}.$ (3.54)

However, we will prove this more generally for so-called partial base
polynomial of the form

$g_{x}^{\otimes j_{1}}\otimes g_{x^{2}}^{\otimes j_{2}}\otimes\ldots\otimes
g_{x^{n-k+1}}^{\otimes j_{n-k+1}}.$ (3.55)

We define partial Bell base polynomials and present relevant results for
commutation matrices in the next section in detail. In principle, this will
show that the ordering of the factors within a base polynomial does not
matter. In the second step, Step 2., we then generalize the recurrence
relation

$D_{k}\textbf{B}^{g}_{3,2}X_{n+1}=D_{k}\left(\textbf{B}^{g}_{2,1}\otimes
g_{x}+2\left(\textbf{B}^{g}_{1,1}\otimes g_{x^{2}}\right)\right)X_{n+1}.$
(3.56)

and show that all recurrence relationships in Lemma 1 hold in the multivariate
case. In the general case, this amounts to showing that

$D_{k}\textbf{B}_{n,k}^{g}X_{n+1}=D_{k}\left[\sum_{i=1}^{n-k+1}{n-1\choose
i-1}\left(\textbf{B}_{n-i,k-1}^{g}\otimes g_{x^{i}}\right)\right]X_{n+1}$
(3.57)

holds. It turns out that we can rewrite any univariate partial Bell polynomial
recurrence into a multivariate partial Bell polynomial recurrence in that
from.

By choosing $D_{k}=f_{g^{k}}$ and $X_{n}=\left(dx\right)^{\otimes n}$, Faà Di
Bruno’s formula for the multivariate case follows using a similar induction
proof as in the univariate case.

##  4 Recurrence of Multivariate Bell Polynomials

In this section we present well-known results concerning the commutation
matrices for Kronecker products, and prove the generalization of multivariate
partial Bell polynomial recurrence relations.

#### Commutation Matrices

In general, two matrices connected by a Kronecker product do not commute.
However one can construct so-called commutation matrices, so that for a given
matrix $A_{1}$ of size $r_{1}\times c_{1}$ and matrix $A_{2}$ of size
$r_{2}\times c_{2}$ we get,

$K_{r_{2}r_{1}}\left(A_{1}\otimes A_{2}\right)K_{c_{1}c_{2}}=A_{2}\otimes
A_{1}.$ (4.58)

$K_{r_{2}r_{1}}$ and $K_{c_{1}c_{2}}$ are commutation matrices of sizes
$r_{2}r_{1}\times r_{2}r_{1}$ and $c_{1}c_{2}\times c_{1}c_{2}$, respectively.
It will also be useful to note some properties of commutation matrices:

  1. 1.

$K_{mn}=K_{nm}$

  2. 2.

$K_{mn}K_{mn}^{\prime}=K_{mn}K_{nm}=I_{mn}$

  3. 3.

$K_{1n}=K_{n1}=I_{n}$.

For a more detailed introduction and an explicit representation of commutation
matrices see Magnus and Neudecker (1979) and Magnus and Neudecker (1999).

One can easily generalize the commutation result for matrices to arbitrarily
many matrices for all cyclical permutations.777See Magnus and Neudecker
(1979), Theorem 3.1 (xi). However, it is also possible to generalize the
commutation result to $m$ Kronecker factor matrices for arbitrary permutations
$\sigma$. Early versions of that result can be found in Holmquist (1985). We
use a recent result in D’Angeli and Donno (2017), as it is closer in notation
to our case, and also delivers some useful corollaries we will make us of.

######  Theorem 2 (Commutation of Kronecker Products, Theorem 3.12. in
D’Angeli and Donno (2017)).

For each $i=1,\ldots,m$, let
$A_{i}=\left\\{a_{i,kl}\right\\}\in\mathcal{M}_{r_{i}\times
c_{i}}\left(\mathbb{R}\right)$ be a matrix with row-index $k=1,\ldots,r_{i}$
and column-index $l=1,\ldots,c_{i}$. Further let $\sigma\in
Sym\left(m\right)$, then there exist so-called shuffle matrices $L$ and $R$
such that

$L\left(A_{1}\otimes A_{2}\otimes\ldots\otimes
A_{m}\right)R=A_{\sigma^{-1}\left(1\right)}\otimes
A_{\sigma^{-1}\left(2\right)}\otimes\ldots\otimes
A_{\sigma^{-1}\left(m\right)},$ (4.59)

where the shuffle matrices are given by

$\displaystyle L=$ $\displaystyle\sum_{\begin{subarray}{c}i_{j}=1,\ldots
r_{j}\\\
j=1,\ldots,k\end{subarray}}E^{i_{\sigma^{-1}\left(1\right)},i_{1}}_{r_{\sigma^{-1}\left(1\right)}\times
r_{1}}\otimes
E^{i_{\sigma^{-1}\left(2\right)},i_{2}}_{r_{\sigma^{-1}\left(2\right)}\times
r_{2}}\otimes\ldots\otimes
E^{i_{\sigma^{-1}\left(m\right)},i_{m}}_{r_{\sigma^{-1}\left(m\right)}\times
r_{m}}$ (4.60) $\displaystyle R=$
$\displaystyle\sum_{\begin{subarray}{c}j_{h}=1,\ldots c_{h}\\\
h=1,\ldots,k\end{subarray}}E^{j_{1},j_{\sigma^{-1}\left(1\right)}}_{c_{1}\times
c_{\sigma^{-1}\left(1\right)}}\otimes
E^{j_{2},j_{\sigma^{-1}\left(2\right)}}_{c_{2}\times
c_{\sigma^{-1}\left(2\right)}}\otimes\ldots\otimes
E^{j_{m},j_{\sigma^{-1}\left(m\right)}}_{c_{m}\times
c_{\sigma^{-1}\left(m\right)}}$ (4.61)

and $E^{i,j}_{m\times n}$ is a elementary matrices of size $m\times n$ and a 1
at position $\left(i,j\right)$, and zeros otherwise.

Similarly useful will be corollary 4.2 in D’Angeli and Donno (2017), which we
will present here in slightly different form.

######  Corollary 1 (Corollary 4.2. in D’Angeli and Donno (2017)).

Define

$K^{\sigma}_{r}:=\sum_{\begin{subarray}{c}i_{j}=1,\ldots r_{j}\\\
j=1,\ldots,k\end{subarray}}E^{i_{\sigma^{-1}\left(1\right)},i_{1}}_{r_{\sigma^{-1}\left(1\right)}\times
r_{1}}\otimes
E^{i_{\sigma^{-1}\left(2\right)},i_{2}}_{r_{\sigma^{-1}\left(2\right)}\times
r_{2}}\otimes\ldots\otimes
E^{i_{\sigma^{-1}\left(m\right)},i_{m}}_{r_{\sigma^{-1}\left(m\right)}\times
r_{m}},$ (4.62)

where $r=\left(r_{1},\ldots,r_{m}\right)$ and $\sigma\in Sym(m)$. Then it
follows that

$K^{\sigma^{-1}}_{r_{\sigma^{-1}}}K^{\sigma}_{r}=I_{\prod_{i}r_{i}},$ (4.63)

with
$r_{\sigma^{-1}}=\left(r_{\sigma^{-1}(1)},\ldots,r_{\sigma^{-1}(m)}\right)$.

$K^{\sigma}_{r}$ is identical to $P^{\sigma}_{r}$ in D’Angeli and Donno
(2017). We use $K$ instead of $P$ to emphasize the similarity to commutation
matrices. Lastly, note that matrices $L$ and $R$ in Theorem 2 are then given
by $L=P^{\sigma}_{r}$ and $R=P^{\sigma^{-1}}_{c_{\sigma^{-1}}}$.

We can now proceed to prepare the proof of step 1 using the following lemma.

######  Lemma 2.

For a permutation $\sigma\in Sym(k)$ and for matrices

$\displaystyle D_{k}$ $\displaystyle=D^{\otimes k}$ (4.64) $\displaystyle
X_{n}$ $\displaystyle=X^{\otimes n},$ (4.65)

where $D\in\mathcal{M}_{1\times n_{y}}\left(\mathbb{R}\right)$ and
$X\in\mathcal{M}_{n_{x}\times 1}\left(\mathbb{R}\right)$, there exist shuffle
matrices $K^{\sigma^{-1}}_{c_{\sigma^{-1}}}$ and $K^{\sigma}_{r}$ such that

$D_{k}K^{\sigma^{-1}}_{c_{\sigma^{-1}}}=D_{k}$ (4.67)

with

$K^{\sigma^{-1}}_{c_{\sigma^{-1}}}=\sum_{\begin{subarray}{c}j_{h}=1,\ldots
n_{y}\\\
h=1,\ldots,k\end{subarray}}E^{j_{1},j_{{\sigma}^{-1}\left(1\right)}}_{n_{y}\times
n_{y}}\otimes
E^{j_{2},j_{{\sigma}^{-1}\left(2\right)}}_{c}\otimes\ldots\otimes
E^{j_{k},j_{{\sigma}^{-1}\left(k\right)}}_{n_{y}\times n_{y}}$ (4.68)

and

$K^{\sigma}_{r}X_{n}=X_{n}$ (4.69)

with

$K^{\sigma}_{r}=\sum_{\begin{subarray}{c}i_{j}=1,\ldots r_{j}\\\
j=1,\ldots,k\end{subarray}}E^{i_{\sigma^{-1}\left(1\right)},i_{1}}_{r_{{\sigma^{-1}\left(1\right)}}\times
r_{1}}\otimes
E^{i_{\sigma^{-1}\left(2\right)},i_{2}}_{r_{{\sigma^{-1}\left(2\right)}}\times
r_{2}}\otimes\ldots\otimes
E^{i_{\sigma^{-1}\left(k\right)},i_{k}}_{r_{{\sigma^{-1}\left(k\right)}}\times
r_{k}},$ (4.70)

###### Proof.

Notice that $D_{k}=D\otimes\ldots\otimes D$ is invariant under any permutation
${\sigma}$, as all factors are identical. We can then write down a shuffle
matrix $K^{\sigma^{-1}}_{c_{\sigma^{-1}}}$ following Theorem 2 by

$K^{\sigma^{-1}}_{c_{\sigma^{-1}}}=\sum_{\begin{subarray}{c}j_{h}=1,\ldots
n_{y}\\\
h=1,\ldots,k\end{subarray}}E^{j_{1},j_{{\sigma}^{-1}\left(1\right)}}_{n_{y}\times
n_{y}}\otimes E^{j_{2},j_{{\sigma}^{-1}\left(2\right)}}_{n_{y}\times
n_{y}}\otimes\ldots\otimes
E^{j_{k},j_{{\sigma}^{-1}\left(k\right)}}_{n_{y}\times n_{y}},$ (4.71)

since $c_{j}=n_{y}$ for all $j$. Because $D_{k}$ is a row vector, the left-
handed shuffle matrix is just 1 and then by construction it holds that
$D_{k}K^{\sigma^{-1}}_{c_{\sigma^{-1}}}=D_{k}$.

A similar argument can be made for $X_{n}$. The only difference here is that
we don’t want to shuffle all individual $n$ matrices but only a $k$ matrix
decomposition out of all $n$ matrices. Given a permutation $\sigma$ the
shuffle is then given by

$K^{\sigma}_{r}=\sum_{\begin{subarray}{c}i_{j}=1,\ldots r_{j}\\\
j=1,\ldots,k\end{subarray}}E^{i_{\sigma^{-1}\left(1\right)},i_{1}}_{r_{{\sigma^{-1}\left(1\right)}}\times
r_{1}}\otimes
E^{i_{\sigma^{-1}\left(2\right)},i_{2}}_{r_{{\sigma^{-1}\left(2\right)}}\times
r_{2}}\otimes\ldots\otimes
E^{i_{\sigma^{-1}\left(k\right)},i_{k}}_{r_{{\sigma^{-1}\left(k\right)}}\times
r_{k}},$ (4.72)

where rows $r_{j}$ depend on the choice of the $k$ matrix decomposition of
$n$. Since $X_{n}$ is a column vector the right-handed shuffle matrix is 1 and
since all matrices are identical, $X_{n}$ is invariant to any shuffle matrix,
so it follows that $X_{n}=K^{\sigma}_{r}X_{n}$.

∎

We proceed by defining the notion of a partial Bell base polynomial. This will
be helpful, as these represent the smallest meaningful part of a Bell
polynomial to which we can apply commutation matrices.

######  Definition 3 (Partial Bell Base Polynomial).

We a define a partial Bell ‘base polynomial’ using

$\tilde{\textbf{B}}^{g}_{n,k,j}:=g_{x}^{\otimes j_{1}}\otimes
g_{x^{2}}^{\otimes j_{2}}\otimes\ldots\otimes g_{x^{n-k+1}}^{\otimes
j_{n-k+1}},$ (4.73)

where $j=\left(j_{1},\ldots,j_{n-k+1}\right)$.

Note that $\tilde{\textbf{B}}^{g}_{n,k,j}$ is a Kronecker product of $k$
matrices because of

$\displaystyle\tilde{\textbf{B}}^{g}_{n,k,j}:=$
$\displaystyle\underbrace{g_{x}\otimes\ldots\otimes
g_{x}}_{j_{1}\;\text{times}}\otimes\ldots\otimes\underbrace{g_{x^{n-k+1}}\otimes\ldots\otimes
g_{x^{n-k+1}}}_{j_{n-k+1}\;\text{times}}$ (4.74) $\displaystyle:=$
$\displaystyle g^{\left(1\right)}\otimes\ldots\otimes
g^{\left(j_{1}\right)}\otimes\ldots\otimes
g^{\left(k-j_{n-k+1}+1\right)}\otimes\ldots\otimes g^{\left(k\right)}$ (4.75)

using the fact that $\sum_{k}j_{l}=k$.

######  Example 1.

With $n=3,k=2,j=\left(1,1\right)$, we get

$\tilde{\textbf{B}}^{g}_{3,2,\left(1,1\right)}=g_{x}\otimes g_{x^{2}}.$ (4.76)

Strictly speaking, we could omit the subindices for $n$ and $k$ in the base
polynomial notation. However, it turns out to be useful to be explicitly
keeping track of which partial Bell polynomial a specific base polynomial
belongs to, without computing that using values in $j$. For example:

$g_{x}\otimes\tilde{\textbf{B}}^{g}_{3,2,\left(1,1\right)}=\tilde{\textbf{B}}^{g}_{4,3,\left(2,1\right)},$
(4.77)

because $j_{1}+j_{2}=3$ and $j_{1}+2j_{2}=4$.

######  Definition 4 (Permutations of Partial Bell Base Polynomials).

Given a permutation $\sigma\in Sym(k)$ over these $k$ matrix factors , then
$\tilde{\textbf{B}}^{g,\sigma}_{n,k,j}$ describes the Bell base polynomial
under a permutation $\sigma$

$\tilde{\textbf{B}}^{g,\sigma}_{n,k,j}:=g^{\sigma\left(1\right)}\otimes\ldots\otimes
g^{\sigma\left(j_{1}\right)}\otimes\ldots\otimes
g^{\sigma\left(k-j_{n-k+1}+1\right)}\otimes\ldots\otimes
g^{\sigma\left(k\right)}$ (4.78)

Based on the previous discussion, it should be immediately clear that in
general
$\tilde{\textbf{B}}^{g}_{n,k,j}\neq\tilde{\textbf{B}}^{g,\sigma}_{n,k,j}$.

We can justify calling $\tilde{\textbf{B}}^{g}_{n,k,j}$ ‘base polynomials’
because partial Bell polynomials can be represented as linear combinations

$\textbf{B}_{n,k}^{g}=\sum_{j}\alpha_{j}\tilde{\textbf{B}}^{g}_{n,k,j},$
(4.79)

where $j=\left(j_{1},\ldots,j_{n-k+1}\right)$ and
$\alpha_{j}=\frac{n!}{j_{1}!j_{2}!\ldots
j_{n-k+1}!\;2^{j_{2}}3^{j_{3}}(n-k+1)^{j_{n-k+1}}}$. The multi-index $j$
iterates over all vectors which satisfy $\sum_{l}j_{l}=k$ and
$\sum_{l}lj_{l}=n$.

Based in this definition and using the previous lemma we are now able to show
that under a simple transformation Bell base polynomial are invariant to
permutations of their ordering.

######  Theorem 3 (Proof of Step 1, Order Invariance of Partial Bell Base
Polynomials).

For Bell base polynomials

$\tilde{\textbf{B}}^{g}_{n,k,j}:=g_{x}^{\otimes j_{1}}\otimes
g_{x^{2}}^{\otimes j_{2}}\otimes\ldots\otimes g_{x^{n-k+1}}^{\otimes
j_{n-k+1}},$ (4.80)

where $j=\left(j_{1},\ldots,j_{n-k+1}\right)$ with $\sum_{l}j_{l}=k$ and
$\sum_{l}lj_{l}=n$, we can shuffle the component matrices arbitrarily such
that

$D_{k}\tilde{\textbf{B}}^{g}_{n,k,j}X_{n}=D_{k}\tilde{\textbf{B}}^{g,\sigma}_{n,k,j}X_{n},$
(4.81)

where $\tilde{\textbf{B}}^{g,\sigma}_{n,k,j}$ is a reshuffling of individual
matrices in $\tilde{\textbf{B}}^{g}_{n,k,j}$ under some permutation $\sigma$
over $\left\\{1,\ldots,k\right\\}$.

###### Proof.

Take the Bell base polynomial of the form

$\tilde{\textbf{B}}^{g}_{n,k,j}:=g_{x}^{\otimes j_{1}}\otimes
g_{x^{2}}^{\otimes j_{2}}\otimes\ldots\otimes g_{x^{n-k+1}}^{\otimes
j_{n-k+1}},$ (4.82)

where $j=\left(j_{1},\ldots,j_{n-k+1}\right)$ and $\sum_{l}j_{l}=k$ and
$\sum_{l}lj_{l}=n$. $\tilde{\textbf{B}}^{g}_{n,k,j}$ now consists of exactly
$k$ individual matrices of the form $g_{x^{l}}$ for $l=1,\ldots,n-k+1$. Each
product component has $n_{y}$ rows and $n_{x}^{l}$ columns. Given some
permutation $\sigma$ over $\left\\{1,\ldots,k\right\\}$, Theorem 2 delivers
shuffle matrices $K^{\sigma}_{r}$ and $K^{\sigma^{-1}}_{c}$:

$K^{\sigma}_{r}=\sum_{\begin{subarray}{c}i_{j}=1,\ldots n_{y}\\\
j=1,\ldots,k\end{subarray}}E^{i_{\sigma^{-1}\left(1\right)},i_{1}}_{n_{y}\times
n_{y}}\otimes E^{i_{\sigma^{-1}\left(2\right)},i_{2}}_{n_{y}\times
n_{y}}\otimes\ldots\otimes
E^{i_{\sigma^{-1}\left(k\right)},i_{k}}_{n_{y}\times n_{y}}$ (4.83)

with

$\left(r_{1},r_{2},\ldots,r_{k}\right)=\left(\underbrace{n_{y},\ldots,n_{y}}_{k\;\text{times}}\right).$
(4.84)

$K^{\sigma^{-1}}_{c_{\sigma^{-1}}}=\sum_{\begin{subarray}{c}j_{h}=1,\ldots
c_{h}\\\
h=1,\ldots,k\end{subarray}}E^{j_{1},j_{\sigma^{-1}\left(1\right)}}_{c_{1}\times
c_{\sigma^{-1}\left(1\right)}}\otimes
E^{j_{2},j_{\sigma^{-1}\left(2\right)}}_{c_{2}\times
c_{\sigma^{-1}\left(2\right)}}\otimes\ldots\otimes
E^{j_{k},j_{\sigma^{-1}\left(k\right)}}_{c_{k}\times
c_{\sigma^{-1}\left(k\right)}},$ (4.85)

with

$\left(c_{1},c_{2},\ldots,c_{k}\right)=\left(\underbrace{n_{x},\ldots,n_{x}}_{j_{1}\;\text{times}},\ldots,\underbrace{n_{x}^{l},\ldots,n_{x}^{l}}_{j_{l}\;\text{times}},\ldots,\underbrace{n_{x}^{n-k+1},\ldots,n_{x}^{n-k+1}}_{j_{n-k+1}\;\text{times}}\right).$
(4.86)

The resulting Kronecker product of matrices using a permutation $\sigma$ is
then given by

$\tilde{\textbf{B}}^{g,\sigma}_{n,k,j}=K^{\sigma}_{r}\tilde{\textbf{B}}^{g}_{n,k,j}K^{\sigma^{-1}}_{c}.$
(4.87)

Following Lemma 2 we can choose columns of the shuffle matrix of $D_{k}$, i.e.
$K^{\sigma^{-1}}_{c_{\sigma^{-1}}}$ such that $c=r$. This is possible, because
in this case the number of rows of all base polynomial factors, $g_{x^{l}}$ is
$n_{y}$ as is the number of columns of $D_{k}$. Using the same permutation
$\sigma$ we therefore get using Corrolary 1

$K^{\sigma^{-1}}_{r_{\sigma^{-1}}}K^{\sigma}_{r}=I_{n_{y}^{k}}.$ (4.88)

Now for $X_{n}=X\otimes\cdots\otimes X$ choose a decomposition of $k$ matrices
such that their rows line up with columns
$\left(c_{1},c_{2},\ldots,c_{k}\right)$ of the $k$ factors of
$\tilde{\textbf{B}}^{g}_{n,k,j}$. This is possible since each individual $X$
is a column vector of size $n_{x}\times 1$, such that we can write for a given
choice of indices $j$:

$X_{n}=\underbrace{\left(X\otimes\ldots\otimes
X\right)}_{j_{1}\;\text{times}}\otimes\ldots\otimes\underbrace{\left(X^{\otimes
l}\otimes\ldots\otimes X^{\otimes
l}\right)}_{j_{l}\;\text{times}}\otimes\ldots\otimes\underbrace{\left(X^{\otimes
n-k+1}\otimes\ldots\otimes X^{\otimes
n-k+1}\right)}_{j_{n-k+1}\;\text{times}}.$ (4.89)

Because by definition $\sum_{l}lj_{l}=n$ and $\sum_{l}j_{l}=k$ and because
each individual $X^{\otimes l}$ has $n_{x}^{l}$ rows, its rows are given by

$\left({r}_{1},{r}_{2},\ldots,{r}_{k}\right)=\left(\underbrace{n_{x},\ldots,n_{x}}_{j_{1}\;\text{times}},\ldots,\underbrace{n_{x}^{l},\ldots,n_{x}^{l}}_{j_{l}\;\text{times}},\ldots,\underbrace{n_{x}^{n-k+1},\ldots,n_{x}^{n-k+1}}_{j_{n-k+1}\;\text{times}}\right).$
(4.90)

Now choose a shuffle matrix for the same permutation $\sigma$ as before given
by Lemma 2 with $K^{\sigma}_{c}$ where we choose $c=r$, where $r$ is given by
definition (4.90). By Corrolary 1 it follows that
$K^{\sigma^{-1}}_{c_{\sigma^{-1}}}K^{\sigma}_{c}=I_{n_{x}^{n}}$. Hence we get

$\displaystyle D_{k}\tilde{\textbf{B}}^{g}_{n,k,j}X_{n}$
$\displaystyle=(D_{k}K^{\sigma^{-1}}_{r_{\sigma^{-1}}})(K^{\sigma}_{r}\tilde{\textbf{B}}^{g}_{n,k,j}K^{\sigma^{-1}}_{c_{\sigma^{-1}}})(K^{\sigma}_{c}X_{n}).$
(4.91)
$\displaystyle=D_{k}\left(K^{\sigma^{-1}}_{r_{\sigma^{-1}}}K^{\sigma}_{r}\right)\tilde{\textbf{B}}^{g,\sigma}_{n,k,j}\left(K^{\sigma^{-1}}_{c_{\sigma^{-1}}}K^{\sigma}_{c}\right)X_{n}$
(4.92) $\displaystyle=D_{k}\tilde{\textbf{B}}^{g,\sigma}_{n,k,j}X_{n},$ (4.93)

which concludes the proof. ∎

Having shown the invariance of Bell base polynomial under pre- and post-
multiplication of $D_{k}$ and $X_{n}$, we can proceed to prove that partial
Bell recurrence relations can similarly be written in that way.

######  Proposition 1 (Proof of Step 2, Invariance of Bell Recurrence
Relations).

The following partial Bell polynomial recurrences hold

  1. i)

$D_{k}\textbf{B}_{n,k}^{g}X_{n}=D_{k}\left[\sum_{i=1}^{n-k+1}{n-1\choose
i-1}\left(\textbf{B}_{n-i,k-1}^{g}\otimes g_{x^{i}}\right)\right]X_{n}$

  2. ii)

$D_{k}\frac{d}{dx}\textbf{B}_{n,k}^{g}X_{n+1}=D_{k}\left[\sum_{i=1}^{n-k+1}{n\choose
i}\left(\textbf{B}_{n-i,k-1}^{g}\otimes g_{x^{i+1}}\right)\right]X_{n+1}$

  3. iii)

$D_{k}\textbf{B}^{g}_{n+1,k}X_{n+1}=D_{k}\left[\left(\textbf{B}^{g}_{n,k-1}\otimes
g_{x}\right)+\frac{\partial\textbf{B}^{g}_{n,k}}{\partial
x^{\prime}}\right]X_{n+1}$

###### Proof.

Here we will prove recurrence relation i) explicitly. Recurrence ii) then
follows using the identical arguments and additionally using the extended
Kronecker product rule in Appendix 7.158. The third recurrence then follows
from Lemma 1 and recurrences i) and ii), as the proof of Lemma 1 does not
require that the Kronecker product commutes. The general idea is to show that
using the previously established results on commutation of Kronecker products
in Theorem 3, these recurrences hold in the multivariate case exactly because
they also hold in the univariate case. To start, recall that the partial Bell
polynomial can be written as

$\textbf{B}_{n,k}^{g}=\sum_{j}\alpha_{j}\tilde{\textbf{B}}^{g}_{n,k,j},$
(4.94)

where $j=\left(j_{1},\ldots,j_{n-k+1}\right)$ and
$\alpha_{j}=\frac{n!}{j_{1}!j_{2}!\ldots
j_{n-k+1}!\;2^{j_{2}}3^{j_{3}}(n-k+1)^{j_{n-k+1}}}$. Multi-index $j$ iterates
over all vectors which satisfy that $\sum_{l}j_{l}=k$ and $\sum_{l}lj_{l}=n$.
It follows that we can write for all $i\in\left\\{1,2,\ldots,n-i-k+2\right\\}$

$\displaystyle\textbf{B}_{n-i,k-1}^{g}\otimes g_{x^{i}}$
$\displaystyle=\left(\sum_{j^{(i)}}\alpha_{j^{(i)}}\tilde{\textbf{B}}^{g}_{n-i,k-1,j^{(i)}}\right)\otimes
g_{x^{i}}$ (4.95)
$\displaystyle=\sum_{j^{(i)}}\alpha_{j^{(i)}}\left(\tilde{\textbf{B}}^{g}_{n-i,k-1,j^{(i)}}\otimes
g_{x^{i}}\right)$ (4.96)
$\displaystyle=\sum_{j^{(i)}}\alpha_{j^{(i)}}\left(g_{x}^{\otimes
j^{(i)}_{1}}\otimes g_{x^{2}}^{\otimes j^{(i)}_{2}}\otimes\ldots\otimes
g_{x^{n-i-k+2}}^{\otimes j^{(i)}_{n-i-k+2}}\otimes g_{x^{i}}\right),$ (4.97)

with $j^{(i)}=\left(j^{(i)}_{1},\ldots,j^{(i)}_{n-i-k+2}\right)$ and
$\sum_{l}j^{(i)}=k-1$ and $\sum_{l}lj^{(i)}=n-i$ . For any $i$, Theorem 2
provides matrices $K_{r_{i}}^{\sigma}$ and
$K_{c_{i,\sigma^{-1}}}^{\sigma^{-1}}$ such that

$\displaystyle g_{x}^{\otimes j^{(i)}_{1}}\otimes\cdots\otimes
g_{x^{n-i-k+2}}(x)^{\otimes j^{(i)}_{n-i-k+2}}\otimes g_{x^{i}}$ (4.98)
$\displaystyle=$ $\displaystyle K_{r_{i}}^{\sigma}\left(g_{x}^{\otimes
j^{(i)}_{1}}\otimes\cdots\otimes
g_{x^{i}}^{\otimes(j^{(i)}_{i}+1)}\otimes\cdots\otimes
g_{x^{n-i-k+2}}(x)^{\otimes
j^{(i)}_{n-i-k+2}}\right)K_{c_{i,\sigma^{-1}}}^{\sigma^{-1}},$ (4.99)

where $r_{i}$ and $c_{i}$ are the respective vectors of rows and columns.
Because $D_{k}K_{c_{i,\sigma^{-1}}}^{\sigma^{-1}}=D_{k}$ and
$K_{r_{i}}^{\sigma}X_{n}=X_{n}$ for all $i$, we can apply Theorem 3, which
then delivers for the RHS of recurrence i) the following:

$\displaystyle D_{k}\left[\sum_{i=1}^{n-k+1}{n-1\choose
i-1}\left(B_{n-i,k-1}^{g}\otimes g_{x^{i}}\right)\right]X_{n}$ (4.100)
$\displaystyle=$ $\displaystyle D_{k}\left[\sum_{i=1}^{n-k+1}{n-1\choose
i-1}\sum_{j^{(i)}}\alpha_{j^{(i)}}\left(g_{x}^{\otimes
j^{(i)}_{1}}\otimes\cdots
g_{x^{i}}^{\otimes(j^{(i)}_{i}+1)}\otimes\cdots\otimes
g_{x^{n-i-k+2}}(x)^{\otimes j^{(i)}_{n-i-k+2}}\right)\right]X_{n}.$ (4.101)

Now, notice that for each $i$, we can write

$\displaystyle\tilde{\textbf{B}}^{g}_{n,k,\tilde{j}^{(i)}}$
$\displaystyle:=g_{x}^{\otimes j_{1}}\otimes\cdots
g_{x^{i}}^{\otimes(j_{i}+1)}\otimes\cdots\otimes g_{x^{n-i-k+2}}(x)^{\otimes
j_{n-i-k+2}}$ (4.103)
$\displaystyle=g_{x}^{\otimes\tilde{j}^{\left(i\right)}_{1}}\otimes\cdots
g_{x^{i}}^{\otimes\tilde{j}^{\left(i\right)}_{i}}\otimes\cdots\otimes
g_{x^{n-i-k+2}}(x)^{\otimes\tilde{j}^{\left(i\right)}_{n-i-k+2}}$ (4.104)

with $\tilde{j}^{\left(i\right)}$ given by

$\tilde{j}^{\left(i\right)}=\left(\tilde{j}^{\left(i\right)}_{1},\ldots,\tilde{j}^{\left(i\right)}_{i},\ldots\tilde{j}^{\left(i\right)}_{n-k+1+(1-i)},\textbf{0}\right)=\left(j^{(i)}_{1},\ldots,j^{(i)}_{i}+1,\ldots
j^{(i)}_{n-k+1+(1-i)},\textbf{0}\right),$ (4.105)

where bold 0 is a row vector of zeros such that $\tilde{j}^{\left(i\right)}$
is of size $1\times(n-k+1)$. It easy to check that for multi-indices
$\tilde{j}^{\left(i\right)}$ it now must hold that

$\sum_{l}\tilde{j}^{\left(i\right)}_{l}=k\quad\text{and}\quad\sum_{l}l\tilde{j}^{\left(i\right)}_{l}=n,$
(4.106)

justifying writing the base polynomial as
$\tilde{\textbf{B}}^{g}_{n,k,\tilde{j}^{(i)}}$.

Hence, if follows

$\displaystyle D_{k}\left[\sum_{i=1}^{n-k+1}{n-1\choose
i-1}\left(g_{x^{i}}\otimes\textbf{B}_{n-i,k-1}^{g}\right)\right]X_{n}$ (4.107)
$\displaystyle=$ $\displaystyle D_{k}\left[\sum_{i=1}^{n-k+1}{n-1\choose
i-1}\sum_{j^{(i)}}\alpha_{j^{(i)}}\tilde{\textbf{B}}^{g}_{n,k,\tilde{j}^{(i)}}\right]X_{n}$
(4.108) $\displaystyle=$ $\displaystyle
D_{k}\left[\sum_{i=1}^{n-k+1}\sum_{j^{(i)}}{n-1\choose
i-1}\alpha_{j^{(i)}}\tilde{\textbf{B}}^{g}_{n,k,\tilde{j}^{(i)}}\right]X_{n}.$
(4.109)

Next, the left-hand-side (LHS) of i) can be written as

$D_{k}\textbf{B}_{n,k}^{g}X_{n}=D_{k}\sum_{j}\alpha_{j}\tilde{\textbf{B}}^{g}_{n,k,{j}}X_{n},$
(4.110)

where $j$ is such that

$\sum_{l}j_{l}=k\quad\text{and}\quad\sum_{l}lj_{l}=n.$ (4.111)

Therefore, by combining equations (4.110) and (4.109) the following must hold
for i) to be true

$D_{k}\left[\sum_{j}\alpha_{j}\tilde{\textbf{B}}^{g}_{n,k,{j}}\right]X_{n}=D_{k}\left[\sum_{i=1}^{n-k+1}\sum_{j^{(i)}}{n-1\choose
i-1}\alpha_{j^{(i)}}\tilde{\textbf{B}}^{g}_{n,k,\tilde{j}^{(i)}}\right]X_{n}.$
(4.112)

Since $\tilde{\textbf{B}}^{g}_{n,k,{j}}$ and
$\tilde{\textbf{B}}^{g}_{n,k,\tilde{j}^{(i)}}$ have the same matrix size, it
would suffice to show that the rearrangement of the double-sum on the RHS
leads to the sum on the LHS.

However, we know that exactly this must be the case, since in the univariate
case it holds that

$\sum_{j}\alpha_{j}\tilde{{B}}^{g}_{n,k,{j}}=\sum_{i=1}^{n-k+1}\sum_{j^{(i)}}{n-1\choose
i-1}\alpha_{j^{(i)}}\tilde{{B}}^{g}_{n,k,\tilde{j}^{(i)}},$ (4.113)

where the non-bold Bell base polynomials $\tilde{{B}}^{g}_{n,k,{j}}$ and
$\tilde{{B}}^{g}_{n,k,\tilde{j}^{(i)}}$ represent their univariate
counterparts. The multivariate recurrence relation i) follows. ∎

##  5 Generalized Multivariate Faà Di Bruno’s Formula

We now prove our multivariate and matrix-valued version of Faà Di Bruno’s
formula by using the invariance of recurrence relations shown in Proposition
1. To do that, we must first identify matrices $D_{k}$ and $X_{n}$ for use in
Faà Di Bruno.

######  Lemma 3.

For composite functions $\left(f\circ
g\right):\mathbb{R}^{n_{x}}\rightarrow\mathbb{R}^{n_{f}}$, by setting
$D_{k}=f_{g^{k}}$ and $X_{n+1}=\left(dx\right)^{\otimes n+1}$ it follows

$f_{g^{k}}\textbf{B}_{n+1,k}^{g}\left(dx\right)^{\otimes
n+1}=f_{g^{k}}\left[\left(g_{x}\otimes\textbf{B}^{g}_{n,k-1}\right)+\frac{\partial\textbf{B}^{g}_{n,k}}{\partial
x^{\prime}}\right]\left(dx\right)^{\otimes n+1}$ (5.114)

###### Proof.

The statement follows from Proposition 1 by setting $D_{k}=f_{g^{k}}$ and
$X_{n+1}=\left(dx\right)^{\otimes n+1}$. Choosing $X=dx$, which is a
$n_{x}\times 1$ column-vector, satisfies the condition in Lemma 2. For
$D_{k}=f_{g^{k}}$ we can write the derivative as

$f_{g^{k}}=f\otimes\underbrace{{\partial\over\partial
g^{\prime}}\otimes\cdots\otimes{\partial\over\partial
g^{\prime}}}_{k\quad\text{times}}.$ (5.115)

Setting $D={\partial\over\partial g^{\prime}}$, we get $f\otimes D^{\otimes
k}$. We know that any permutation of base polynomials in
$\textbf{B}_{n+1,k}^{g}$ will be invariant with respect to $D^{\otimes k}$.
However the same holds with respect to $f\otimes D^{\otimes k}$. This follows
from Lemma 2 with

$\displaystyle f\otimes D^{\otimes k}$ $\displaystyle=f\otimes D^{\otimes
k}K^{\sigma^{-1}}_{c}$ (5.116) $\displaystyle=\left(f\cdot
1\right)\otimes\left(D^{\otimes k}K^{\sigma^{-1}}_{c}\right)$ (5.117)
$\displaystyle=\left(f\otimes D^{\otimes k}\right)\left(1\otimes
K^{\sigma^{-1}}_{c}\right)$ (5.118) $\displaystyle=\left(f\otimes D^{\otimes
k}\right)K^{\sigma^{-1}}_{c}.$ (5.119)

∎

Taking everything together that we have shown so far, multivariate Faà Di
Bruno now follows using a similar proof of induction as in the univariate
case.

######  Theorem 4 (Multivariate Faà Di Bruno).

For some $n,k\in\mathbb{N}$ with $n\geq k$ and function
$g:\mathbb{R}^{n_{x}}\rightarrow\mathbb{R}^{n_{y}}$ which is at least $n-k+1$
differentiable on some domain $D_{g}\subset\mathbb{R}^{n_{x}}$ and a function
$f:\mathbb{R}^{n_{y}}\rightarrow\mathbb{R}^{n_{f}}$ which is at least $n$
times differentiable on some domain $D_{f}\subset\mathbb{R}^{n_{y}}$, the
$n$-th order total derivative of $f\circ g$ is given by Faà Di Bruno’s
formula:

$\displaystyle
d^{n}f(g(x))=\sum_{k=1}^{n}f_{g^{k}}(g(x))\textbf{B}^{g}_{n,k}(dx)^{\otimes
n},$ (5.120)

where $f_{g^{k}}(g(x))$ denotes the $n_{f}\times n_{y}^{k}$ matrix of the
$k$th-order partial derivatives of $f(\cdot)$ and similarly $g_{x^{l}}(x)$
denotes the $n_{y}\times n_{x}^{l}$ matrix of the $l$th-order partial
derivatives of $g(\cdot)$.

###### Proof.

We prove the multivariate case again by induction. Using the apparatus set up
in this paper, the proof follows closely the univariate one. We leave the
induction beginning to the reader and directly start by assuming that (5.120)
holds for some arbitrary $n\in\mathbb{N}$. Then it follows using the
differential version of the product rule in Appendix A (7.167) and the
differential version of the chain rule in Appendix A (7.169):

$\displaystyle d^{n+1}f(g(x))$ $\displaystyle=d(d^{n}f(g(x)))$ (5.121)
$\displaystyle=d(\sum_{k=1}^{n}f_{g^{k}}(g(x))\textbf{B}^{g}_{n,k}(dx)^{\otimes
n})$ (5.122)
$\displaystyle=\sum_{k=1}^{n}d(f_{g^{k}}(g(x))\textbf{B}^{g}_{n,k}(dx)^{\otimes
n})$ (5.123)
$\displaystyle=\sum_{k=1}^{n}f_{g^{k+1}}\left(I_{n_{y}^{k}}\otimes
g_{x}\right)\left(\textbf{B}^{g}_{n,k}(dx)^{\otimes n}\otimes
dx\right)+f_{g^{k}}\frac{\partial\textbf{B}^{g}_{n,k}}{\partial
x^{\prime}}\left(dx\right)^{\otimes n+1}$ (5.124)
$\displaystyle=\sum_{k=1}^{n}f_{g^{k+1}}\left(\textbf{B}^{g}_{n,k}\otimes
g_{x}\right)\left(dx\right)^{\otimes
n+1}+f_{g^{k}}\frac{\partial\textbf{B}^{g}_{n,k}}{\partial
x^{\prime}}\left(dx\right)^{\otimes n+1}$ (5.125)

By factoring for $f_{g^{k}}\left(g(x)\right)$, re-parameterizing the sums and
noting that
$\textbf{B}^{g}_{n,0}\left(\cdot\right)=\textbf{B}^{g}_{n,n+1}\left(\cdot\right)$
= 0, we get

$\displaystyle=\sum_{k=1}^{n+1}f_{g^{k}}\left[\left(\textbf{B}^{g}_{n,k-1}\otimes
g_{x}\right)+\frac{\partial\textbf{B}^{g}_{n,k}}{\partial
x^{\prime}}\right]\left(dx\right)^{\otimes n+1}.$ (5.126)

Lemma 3 then provides, that for each $k$

$\sum_{k=1}^{n+1}f_{g^{k}}\textbf{B}^{g}_{n+1,k}\left(dx\right)^{\otimes
n+1}=\sum_{k=1}^{n+1}f_{g^{k}}\left[\left(\textbf{B}^{g}_{n,k-1}\otimes
g_{x}\right)+\frac{\partial\textbf{B}^{g}_{n,k}}{\partial
x^{\prime}}\right]\left(dx\right)^{\otimes n+1}$ (5.127)

holds, which completes the proof. ∎

##  6 Applications

###  6.1 Computing Moments of Multivariate Distributions

For a $k\times 1$ random variable $X$ we can compute its moments using the
moment generating function (MGF)

$M_{X}\left(t\right)=E_{X}\left[\exp{t^{\prime}X}\right],$ (6.128)

where $t$ is a $k\times 1$ vector. The resulting series expansion can then be
written as

$M_{X}\left(t\right)=\sum_{k=0}^{\infty}{1\over k!}m_{k}t^{\otimes k},$
(6.129)

with $m_{n}=E\left[\left(X^{\prime}\right)^{\otimes n}\right]$. The moments of
order $n$ are then given by

$m_{n}=\left.{\partial^{n}M_{X}\over\partial t^{\otimes
n}}\right|_{t=\textbf{0}}.$ (6.130)

#### Example: (Multivariate Normal Distribution)

For a $k\times 1$ random variable $X$ which is normally distributed with mean
vector $\mu$ and covariance matrix $\Sigma$, the multivariate MGF can be
computed to be

$M_{X}\left(t\right)=\exp{\left(t^{\prime}\mu+{1\over 2}t^{\prime}\Sigma
t\right)}.$ (6.131)

To compute its $n$-th derivative we again make use of our novel version of Faà
Di Bruno’s formula given by

$\displaystyle d^{n}(f\circ
g)=\sum_{k=1}^{n}f_{g^{k}}(g(x))\textbf{B}^{g}_{n,k}(dx)^{\otimes n}.$ (6.132)

However, in this case we will make use not of the differential form but of the
resulting matrix of partial derivatives given by

${\partial^{n}(f\circ g)\over\partial(x^{\prime})^{\otimes
n}}=\sum_{k=1}^{n}f_{g^{k}}(g(x))\textbf{B}^{g}_{n,k}.$ (6.133)

Note that this representation might not be a unique representation for the
matrix of derivatives.888For more details on this issue, see Appendix B. We
will deal with this issue towards the end of this application.

Now, setting $f(y)=\exp{y}$ with $f:\mathbb{R}\rightarrow\mathbb{R}$ and
$g(t)=t^{\prime}\mu+{1\over 2}t^{\prime}\Sigma t$ with
$g:\mathbb{R}^{k}\rightarrow\mathbb{R}$, we get $f_{g^{j}}=\exp{(g)}$ for all
$j$ and

$\displaystyle g_{t}$ $\displaystyle=\mu^{\prime}+{1\over
2}{d\left(t^{\prime}\Sigma t\right)\over dt}$ (6.134) $\displaystyle
g_{t^{j}}$ $\displaystyle={1\over 2}{d^{j}\left(t^{\prime}\Sigma t\right)\over
dt^{j}}\quad\text{for all}\quad j\geq 2.$ (6.135)

Evaluated at $t=0$, we get $f_{g^{j}}=1$ and therefore we can write moments
using

$m_{n}=\sum_{j=1}^{n}\left.\textbf{B}_{n,j}^{g\left(t\right)}\right|_{t=0}.$
(6.136)

Since only derivatives of $g$ up to the second order remain non-zero when
evaluated at $t=0$ , i.e.

$\displaystyle\left.g_{t}\right|_{t=0}$ $\displaystyle=\mu^{\prime}$ (6.137)
$\displaystyle\left.g_{t^{2}}\right|_{t=0}$
$\displaystyle=vec\left(\Sigma\right)^{\prime}:=\left(\sigma_{11},\sigma_{12},\sigma_{12},\sigma_{22}\right)$
(6.138) $\displaystyle\left.g_{t^{j}}\right|_{t=0}$
$\displaystyle=0\quad\text{for all}\quad j\geq 3,$ (6.139)

the sum reduces to

$m_{n}=\sum_{j=0}^{\left\lfloor{{n\over
2}}\right\rfloor}\textbf{B}_{n,n-j}\left(\mu^{\prime},vec\left(\Sigma\right)^{\prime},0,\cdots,0\right).$
(6.140)

Therefore we arrive at

$m_{n}=\sum_{j=0}^{\left\lfloor{{n\over
2}}\right\rfloor}{n!\over(n-2j)!(j)!2^{j}}\left((\mu^{\prime})^{\otimes(n-2j)}\otimes\left(vec\left(\Sigma\right)^{\prime}\right)^{\otimes
j}\right).$ (6.141)

Notice that for centered moments with $\mu=0$, we have to evaluate the
expression only for $j={n\over 2}$. Hence get for even $n$

$m_{n}={n!\over{n\over 2}!2^{n\over 2}}vec(\Sigma)^{\prime\otimes{n\over 2}}.$
(6.142)

These expressions are almost identical to Holmquist (1988) except for his use
of the symmetrizer matrix $S_{k,n}$ and the use of columns vectors. When
defining the representation based on row vectors, the moment vector in
Holmquist (1988) is given by

$m_{n}=\sum_{j=0}^{\left\lfloor{{n\over
2}}\right\rfloor}{n!\over(n-2j)!(j)!2^{j}}\left((\mu^{\prime})^{\otimes(n-2j)}\otimes\left(vec\left(\Sigma\right)^{\prime}\right)^{\otimes
j}\right)S_{k,n}.$ (6.143)

In essence the symmetrizer matrix ensures that each component of the moment
vector represents the exact value of one of the many possible scalar moments.
To illustrate this point and therefore the relevance of the symmetrizer
matrix, we compute $m_{4}$ with $\mu=0$ and $k=2$ using equation (6.142),
which yields

$m_{4}=3vec(\Sigma)^{\prime\otimes{2}}.$ (6.144)

In the series expansion of the moment generating function we get

$M_{X}\left(t\right)=\sum_{k=1}^{\infty}{1\over k!}m_{k}t^{\otimes k}.$
(6.145)

Therefore, we show the moment vector $m_{4}$ together with the corresponding
vector $t^{\otimes 4}$:

$m_{4}=\left(\begin{array}[]{c}3\,{\sigma_{11}}^{2}\\\
3\,\sigma_{11}\,\sigma_{12}\\\ 3\,\sigma_{11}\,\sigma_{12}\\\
3\,\sigma_{11}\,\sigma_{22}\\\ 3\,\sigma_{11}\,\sigma_{12}\\\
3\,{\sigma_{12}}^{2}\\\ 3\,{\sigma_{12}}^{2}\\\ 3\,\sigma_{12}\,\sigma_{22}\\\
3\,\sigma_{11}\,\sigma_{12}\\\ 3\,{\sigma_{12}}^{2}\\\ 3\,{\sigma_{12}}^{2}\\\
3\,\sigma_{12}\,\sigma_{22}\\\ 3\,\sigma_{11}\,\sigma_{22}\\\
3\,\sigma_{12}\,\sigma_{22}\\\ 3\,\sigma_{12}\,\sigma_{22}\\\
3\,{\sigma_{22}}^{2}\end{array}\right)^{\prime},\quad t^{\otimes
4}=\left(\begin{array}[]{c}{t_{1}}^{4}\\\ {t_{1}}^{3}\,t_{2}\\\
{t_{1}}^{3}\,t_{2}\\\ {t_{1}}^{2}\,{t_{2}}^{2}\\\ {t_{1}}^{3}\,t_{2}\\\
{t_{1}}^{2}\,{t_{2}}^{2}\\\ {t_{1}}^{2}\,{t_{2}}^{2}\\\ t_{1}\,{t_{2}}^{3}\\\
{t_{1}}^{3}\,t_{2}\\\ {t_{1}}^{2}\,{t_{2}}^{2}\\\ {t_{1}}^{2}\,{t_{2}}^{2}\\\
t_{1}\,{t_{2}}^{3}\\\ {t_{1}}^{2}\,{t_{2}}^{2}\\\ t_{1}\,{t_{2}}^{3}\\\
t_{1}\,{t_{2}}^{3}\\\ {t_{2}}^{4}\end{array}\right).$ (6.146)

The elements of $m_{4}$ corresponding to a specific element of $t^{\otimes 4}$
should deliver the corresponding individual scalar moment. So for example, the
entry of $m_{4}$ corresponding to $t_{1}^{2}t_{2}^{2}$ should be equal to the
scalar moment $E\left[x_{1}^{2}x_{2}^{2}\right]$. However, as the entries of
$t_{1}^{2}t_{2}^{2}$ appear multiple times, one can notice that the
corresponding entries in $m_{4}$ aren’t identical. To arrive at the correct
moment we average over all values appearing with $t_{1}^{2}t_{2}^{2}$. By
averaging over these entries, we get

$m_{4}=\left(\begin{array}[]{c}3\,{\sigma_{11}}^{2}\\\
3\,\sigma_{11}\,\sigma_{12}\\\ 3\,\sigma_{11}\,\sigma_{12}\\\
2\,{\sigma_{12}}^{2}+\sigma_{11}\,\sigma_{22}\\\
3\,\sigma_{11}\,\sigma_{12}\\\
2\,{\sigma_{12}}^{2}+\sigma_{11}\,\sigma_{22}\\\
2\,{\sigma_{12}}^{2}+\sigma_{11}\,\sigma_{22}\\\
3\,\sigma_{12}\,\sigma_{22}\\\ 3\,\sigma_{11}\,\sigma_{12}\\\
2\,{\sigma_{12}}^{2}+\sigma_{11}\,\sigma_{22}\\\
2\,{\sigma_{12}}^{2}+\sigma_{11}\,\sigma_{22}\\\
3\,\sigma_{12}\,\sigma_{22}\\\
2\,{\sigma_{12}}^{2}+\sigma_{11}\,\sigma_{22}\\\
3\,\sigma_{12}\,\sigma_{22}\\\ 3\,\sigma_{12}\,\sigma_{22}\\\
3\,{\sigma_{22}}^{2}\end{array}\right)^{\prime},$ (6.147)

which implies
$E\left[x_{1}^{2}x_{2}^{2}\right]=2\sigma_{12}^{2}+\sigma_{11}\sigma_{22}$.
This moment vector $m_{4}$ is now identical to the one computed using the
formula in Holmquist (1988). The symmetrizer matrix in Holmquist (1988)
ensures that all entries that correspond to the same individual scalar moment
deliver the exact same value. The fundamental issue here, is that the matrix
of derivative is not uniquely defined by

${\partial^{n}(f\circ g)\over\partial(x^{\prime})^{\otimes
n}}=\sum_{k=1}^{n}f_{g^{k}}(g(x))\textbf{B}^{g}_{n,k}.$ (6.148)

Instead by using

${\partial^{n}(f\circ g)\over\partial(x^{\prime})^{\otimes
n}}=\sum_{k=1}^{n}f_{g^{k}}(g(x))\textbf{B}^{g}_{n,k}S_{n_{x},n},$ (6.149)

where $S_{n,k}$ is the symmetrizer matrix for n vectors of size $n_{x}\times
1$, we would arrive at a representation identical to Holmquist (1988).999See
Appendix B for more details. In practice, the computation of the full
symmetrizer matrix might not be necessary. In our application, it would
suffice to identify the elements of $m_{n}$ which are not identical with
respect to the corresponding elements of $t^{\otimes n}$ and then averaging
over all realizations.

##  7 Concluding Remarks

In this paper, we introduce a novel generalization of the partial Bell
polynomial. This generalization is similar in structure to the univariate
version, except for the replacement of the scalar multiplication operation
with the Kronecker product. This extension present us with several challenges,
because the Kronecker product doesn’t generally commute. Nevertheless, we are
able to show that under a simple transformation well-known recurrences
relations for the partial Bell polynomial continue to hold. It is exactly
because of that insight that we can prove a generalization of Faà Di Bruno’s
formula which is again similar in structure to the univariate version. In our
application we show that using such a representation of higher-order composite
derivatives, the computation of higher-order moments of the multivariate
normal distribution becomes a simple task. This representation has further
potential applications in multivariate statistics and in the theory and
computation of dynamic stochastic general equilibrium models in economics,
which we attempt to explore in the future.

## References

  * (1)
  * Chacón and Duong (2021) Chacón, J. E., and T. Duong (2021): “Higher order differential analysis with vectorized derivatives,” arXiv:2011.01833 [math].
  * Charalambides (2018) Charalambides, C. A. (2018): _Enumerative combinatorics_. Chapman and Hall/CRC.
  * Comtet (1974) Comtet, L. (1974): _Advanced Combinatorics: The art of finite and infinite expansions_. Springer Science & Business Media.
  * Constantine and Savits (1996) Constantine, G., and T. Savits (1996): “A Multivariate Faa di Bruno Formula with Applications,” _Transactions of the American Mathematical Society_ , 348(2), 503–520.
  * D’Angeli and Donno (2017) D’Angeli, D., and A. Donno (2017): “Shuffling matrices, kronecker product and discrete fourier transform,” _Discrete Applied Mathematics_ , 233, 1–18.
  * Di Nardo, Guarino, and Senato (2011) Di Nardo, E., G. Guarino, and D. Senato (2011): “A new algorithm for computing the multivariate Faa di Bruno’s formula,” _Applied mathematics and computation_ , 217(13), 6286–6295.
  * Gzyl (1986) Gzyl, H. (1986): “Multidimensional extension of Faa di Bruno’s formula,” _Journal of Mathematical Analysis and Applications_ , 116(2), 450–455.
  * Hardy (2006) Hardy, M. (2006): “Combinatorics of Partial Derivatives,” arXiv:math/0601149.
  * Holmquist (1985) Holmquist, B. (1985): “The direct product permuting matrices,” _Linear and Multilinear Algebra_ , 17(2), 117–141.
  * Holmquist (1988) Holmquist, B. (1988): “Moments and cumulants of the multivariate normal distribution,” _Stochastic Analysis and Applications_ , 6(3), 273–278.
  * Holmquist (1996) Holmquist, B. (1996): “The d-variate vector hermite polynomial of order k,” _Linear Algebra and its Applications_ , 237-238, 155–190.
  * Johnson (2002) Johnson, W. P. (2002): “The curious history of Faà di Bruno’s formula,” _The American mathematical monthly_ , 109(3), 217–234.
  * Leipnik and Pearce (2007) Leipnik, R. B., and C. E. M. Pearce (2007): “The multivariate Faà di Bruno formula and multivariate Taylor expansions with explicit integral remainder term,” _The ANZIAM Journal_ , 48(3), 327–341.
  * Magnus (2010) Magnus, J. R. (2010): “On the concept of matrix derivative,” _Journal of Multivariate Analysis_ , 101(9), 2200–2206.
  * Magnus and Neudecker (1979) Magnus, J. R., and H. Neudecker (1979): “The commutation matrix: some properties and applications,” _The annals of statistics_ , 7(2), 381–394.
  * Magnus and Neudecker (1999) Magnus, J. R., and H. Neudecker (1999): _Matrix differential calculus with applications in statistics and econometrics_. John Wiley, New York, rev. ed edn.
  * Mishkov (2000) Mishkov, R. L. (2000): “Generalization of the formula of Faa di Bruno for a composite function with a vector argument,” _International Journal of Mathematics and Mathematical Sciences_ , 24(7), 481–491.
  * Noschese and Ricci (2003) Noschese, S., and P. E. Ricci (2003): “Differentiation of multivariable composite functions and Bell polynomials,” _Journal of Computational Analysis and Applications_ , 5, 333–340.
  * Rogers (1980) Rogers, G. (1980): _Matrix Derivatives_. Marcel Dekker, New York.
  * Schumann (2019) Schumann, A. (2019): “Multivariate Bell Polynomials and Derivatives of Composed Functions,” arXiv:1903.03899 [math].
  * Sultan and Tracy (1996) Sultan, S. A., and D. S. Tracy (1996): “Moments of the complex multivariate normal distribution,” _Linear algebra and its applications_ , 237, 191–204.
  * Tracy and Sultan (1993) Tracy, D. S., and S. A. Sultan (1993): “Higher order moments of multivariate normal distribution using matrix derivatives,” _Stochastic Analysis and Applications_ , 11(3), 337–348.

## Appendix

### A. Derivatives and Differentials

In this section, we present and derive representations and derivative rules
vector- and matrix-valued derivatives and differentials. In that we mostly
follow Magnus and Neudecker (1999) and Rogers (1980).

We start by defining the vector- and matrix-valued representations of
derivative and differentials. For this purpose, let matrices $F,G,H$ be
functions of a vector $x$ of size $n_{x}\times 1$. The matrix $F$ is of size
$s\times t$, $G$ of size $y\times 1$, and $H$ of size $p\times q$. We define a
differential operator ${\partial\over\partial x^{\prime}}$ w.r.t the transpose
of $x$ by

${\partial\over\partial x^{\prime}}=\left({\partial\over\partial
x_{1}},\ldots,{\partial\over\partial x_{n_{x}}}\right).$ (7.150)

First-order derivatives of matrices and vectors are then defined by

$F_{x}={\partial F\over\partial x^{\prime}}:=F\otimes{\partial\over\partial
x^{\prime}}.$ (7.151)

Analogously, higher-order derivatives are defined by

$F_{x^{k}}={\partial F\over\partial(x^{\prime})^{\otimes
k}}:=F\otimes\left({\partial\over\partial x^{\prime}}\right)^{\otimes k}.$
(7.152)

This notation implies some structure on the ordering of partial derivative
within ${\partial F\over\partial x^{\prime}}$.101010See Tracy and Sultan
(1993), Sultan and Tracy (1996), and Rogers (1980) who also make use of this
notation. More explicitly this structure corresponds to the replacement of
each component $F_{i,j}$ of $F$ by its respective row vector of derivatives,
i.e.

${\partial F_{i,j}\over\partial x^{\prime}}=\left({\partial
F_{i,j}\over\partial x_{1}},\ldots,{\partial F_{i,j}\over\partial
x_{n_{x}}}\right).$ (7.153)

When F is a vector, denoted by $f$ in the following, we can write the
corresponding first-order differential:

#### The First-order Differential

Let $f(x)$ be a differentiable $s\times 1$ vector function of an $n_{x}\times
1$ vector of real variables $x$. Then we can use the notation of the
derivative matrix to obtain

$\displaystyle df=\frac{\partial f}{\partial x}dx=f_{x}dx$

as the first differential of $f(\cdot)$ at $x$ with increment $dx$. In more
detail,

$\displaystyle df(x)=\begin{bmatrix}df_{1}(x)\\\ df_{2}(x)\\\ \vdots\\\
df_{s}(x)\end{bmatrix}$ $\displaystyle=\begin{bmatrix}\frac{\partial
f_{1}(x)}{\partial x_{1}}dx_{1}+\frac{\partial f_{1}(x)}{\partial
x_{2}}dx_{2}+\dots+\frac{\partial f_{1}(x)}{\partial x_{n_{x}}}dx_{n_{x}}\\\
\frac{\partial f_{2}(x)}{\partial x_{1}}dx_{1}+\frac{\partial
f_{2}(x)}{\partial x_{2}}dx_{2}+\dots+\frac{\partial f_{2}(x)}{\partial
x_{n_{x}}}dx_{n_{x}}\\\ \vdots\\\ \frac{\partial f_{s}(x)}{\partial
x_{1}}dx_{1}+\frac{\partial f_{s}(x)}{\partial
x_{2}}dx_{2}+\dots+\frac{\partial f_{s}(x)}{\partial x_{n_{x}}}dx_{n_{x}}\\\
\end{bmatrix}$ $\displaystyle=\begin{bmatrix}\frac{\partial f_{1}(x)}{\partial
x_{1}}&\frac{\partial f_{1}(x)}{\partial x_{2}}&\dots&\frac{\partial
f_{1}(x)}{\partial x_{n_{x}}}\\\ \frac{\partial f_{2}(x)}{\partial
x_{1}}&\ddots&&\vdots\\\ \vdots&&\ddots&\vdots\\\ \frac{\partial
f_{s}(x)}{\partial x_{1}}&\dots&\dots&\frac{\partial f_{s}(x)}{\partial
x_{n_{x}}}\end{bmatrix}\begin{bmatrix}dx_{1}\\\ dx_{2}\\\ \vdots\\\
dx_{n_{x}}\end{bmatrix}$

Using a more concise notation we get,

$\displaystyle df$ $\displaystyle=\begin{bmatrix}f_{1,x}\\\ f_{2,x}\\\
\vdots\\\ f_{s,x}\end{bmatrix}dx=f_{x}d_{x}.$

#### The Second-order Differential

Let $f(x)$ be a twice differentiable $s\times 1$ vector function of an
$n_{x}\times 1$ vector of real variables $x$. Then we can use the notation of
the derivative matrix to obtain

$\displaystyle d^{2}f(x)=\frac{\partial\left(\frac{\partial f(x)}{\partial
x}dx\right)}{\partial x}dx=\frac{\partial^{2}f(x)}{\partial x^{2}}dx^{\otimes
2}=f_{x^{2}}dx^{\otimes 2}$

as the second differential of $f(\cdot)$ at $x$ with increment $dx$. In more
detail,

$\displaystyle d^{2}f(x)$
$\displaystyle=d(df(x))=\begin{bmatrix}d(df_{1}(x))\\\ d(df_{2}(x))\\\
\vdots\\\ d(df_{s}(x))\end{bmatrix}=\begin{bmatrix}d(f_{1,x}dx)\\\
d(f_{2,x}dx)\\\ \vdots\\\
d(f_{s,x}dx)\end{bmatrix}=\begin{bmatrix}\begin{bmatrix}df_{1,x_{1}}&df_{1,x_{2}}&\dots&df_{1,x_{n_{x}}}\end{bmatrix}dx\\\
\begin{bmatrix}df_{2,x_{1}}&df_{2,x_{2}}&\dots&df_{2,x_{n_{x}}}\end{bmatrix}dx\\\
\vdots\\\
\begin{bmatrix}df_{s,x_{1}}&df_{s,x_{2}}&\dots&df_{s,x_{n_{x}}}\end{bmatrix}dx\end{bmatrix}$
$\displaystyle=\begin{bmatrix}f_{1,x_{1},x}dx&f_{1,x_{2},x}dx&\dots&f_{1,x_{n_{x}},x}dx\\\
f_{2,x_{1},x}dx&f_{2,x_{2},x}dx&\dots&f_{2,x_{n_{x}},x}dx\\\
\vdots&\vdots&\ddots&\vdots\\\
f_{s,x_{1},x}dx&f_{s,x_{2},x}dx&\dots&f_{s,x_{n_{x}},x}dx\\\ \end{bmatrix}dx$
$\displaystyle=\begin{bmatrix}f_{x_{1},x}dx&f_{x_{2},x}dx&\dots&f_{x_{n_{x}},x}dx\end{bmatrix}dx$
$\displaystyle=\begin{bmatrix}f_{x_{1},x}&f_{x_{2},x}&\dots&f_{x_{n_{x}},x}\end{bmatrix}\begin{bmatrix}dx&\operatorname*{0}\limits_{n_{x}\times
1}&\dots&\operatorname*{0}\limits_{n_{x}\times 1}\\\
\operatorname*{0}\limits_{n_{x}\times
1}&dx&\dots&\operatorname*{0}\limits_{n_{x}\times 1}\\\
\vdots&\vdots&\ddots&\cdots\\\ \operatorname*{0}\limits_{n_{x}\times
1}&\operatorname*{0}\limits_{n_{x}\times 1}&\dots&dx\end{bmatrix}dx$
$\displaystyle=f_{x,x}\left(I_{n_{x}}\otimes
dx\right)dx=f_{x^{2}}\left(dx\otimes dx\right)=f_{x^{2}}dx^{\otimes 2}.$

#### The $k$th-order Differential

Let $f(x)$ be a $k$ times differentiable $s\times 1$ vector function of an
$n_{x}\times 1$ vector of real variables $x$. Then

$\displaystyle
d^{k}f(x)=\frac{\partial\left(\frac{\partial^{k-1}f(x)}{\partial
x^{k-1}}dx^{\otimes k-1}\right)}{\partial
x}dx=\frac{\partial^{k}f(x)}{\partial x^{k}}dx^{\otimes
k}=f_{x^{k}}dx^{\otimes k}$

as the $k$-order differential of $f(\cdot)$ at $x$ with increment $dx$. In
more detail,

$\displaystyle d^{k}f(x)$
$\displaystyle=d(d^{k-1}f(x))=\begin{bmatrix}d(d^{k-1}f_{1}(x))\\\
d(d^{k-1}f_{2}(x))\\\ \vdots\\\
d(d^{k-1}f_{s}(x))\end{bmatrix}=\begin{bmatrix}d(f_{1,x^{k-1}}dx^{\otimes
k-1})\\\ d(f_{2,x^{k-1}}dx^{\otimes k-1})\\\ \vdots\\\
d(f_{s,x^{k-1}}dx^{\otimes k-1})\end{bmatrix}$
$\displaystyle=\begin{bmatrix}\begin{bmatrix}df_{1,x_{1},x^{k-1}}&df_{1,x_{2},x^{k-1}}&\dots&df_{1,x_{n_{x},x^{k-1}}}\end{bmatrix}dx^{\otimes
k-1}\\\
\begin{bmatrix}df_{2,x_{1},x^{k-1}}&df_{2,x_{2},x^{k-1}}&\dots&df_{2,x_{n_{x},x^{k-1}}}\end{bmatrix}dx^{\otimes
k-1}\\\ \vdots\\\
\begin{bmatrix}df_{s,x_{1},x^{k-1}}&df_{s,x_{2},x^{k-1}}&\dots&df_{s,x_{n_{x},x^{k-1}}}\end{bmatrix}dx^{\otimes
k-1}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}\begin{bmatrix}f_{1,x_{1},x^{k-1}}dx&f_{1,x_{2},x^{k-1}}dx&\dots&f_{1,x_{n_{x}},x^{k-1}}dx\end{bmatrix}\\\
\begin{bmatrix}f_{2,x_{1},x^{k-1}}dx&f_{2,x_{2},x^{k-1}}dx&\dots&f_{2,x_{n_{x}},x^{k-1}}dx\end{bmatrix}\\\
\vdots\\\
\begin{bmatrix}f_{s,x_{1},x^{k-1}}dx&f_{s,x_{2},x^{k-1}}dx&\dots&f_{s,x_{n_{x}},x^{k-1}}dx\end{bmatrix}\\\
\end{bmatrix}dx^{\otimes k-1}$
$\displaystyle=\begin{bmatrix}f_{x_{1},x^{k-1}}dx&f_{x_{2},x^{k-1}}dx&\dots&f_{x_{n_{x}},x^{k-1}}dx\end{bmatrix}dx^{\otimes
k-1}$
$\displaystyle=\begin{bmatrix}f_{x_{1},x^{k-1}}&f_{x_{2},x^{k-1}}&\dots&f_{x_{n_{x}},x^{k-1}}\end{bmatrix}\begin{bmatrix}dx&\operatorname*{0}\limits_{n_{x}\times
1}&\dots&\operatorname*{0}\limits_{n_{x}\times 1}\\\
\operatorname*{0}\limits_{n_{x}\times
1}&dx&\dots&\operatorname*{0}\limits_{n_{x}\times 1}\\\
\vdots&\vdots&\ddots&\cdots\\\ \operatorname*{0}\limits_{n_{x}\times
1}&\operatorname*{0}\limits_{n_{x}\times 1}&\dots&dx\end{bmatrix}dx^{\otimes
k-1}$ $\displaystyle=f_{x,x^{k-1}}\left(I_{n_{x}^{k-1}}\otimes
dx\right)dx^{\otimes k-1}=f_{x^{k}}\left(dx^{\otimes k-1}\otimes
dx\right)=f_{x^{k}}dx^{\otimes k}.$

Similarly, we can define a matrix differential by

#### Matrix Differential

$dF={\partial F\over\partial x^{\prime}}\left(I_{t}\otimes dx\right).$ (7.154)

We can now present rules of differentiation. We omit the addition of rule as
it is identical to the univariate case.

#### Product Rule

If the matrix product $FH$ exists, i.e. $t=p$, then

${\partial FH\over\partial x^{\prime}}={\partial F\over\partial
x^{\prime}}\left(H\otimes I_{n_{x}}\right)+F{\partial H\over\partial
x^{\prime}}.$ (7.155)

#### Kronecker Product Rule

$\displaystyle{\partial\left(F\otimes H\right)\over\partial x^{\prime}}=$
$\displaystyle F\otimes{\partial H\over\partial
x^{\prime}}+K_{s,p}\left(H\otimes{\partial F\over\partial
x^{\prime}}\right)\left(K_{q,t}\otimes I_{n_{x}}\right)$ (7.156)
$\displaystyle=$ $\displaystyle F\otimes{\partial H\over\partial
x^{\prime}}+\left({\partial F\over\partial x^{\prime}}\otimes
H\right)\left(I_{t}\otimes K_{n_{x},q}\right)$ (7.157)

The second equality follows either from the application of the commutation of
the matrices in the parentheses and some computation using commutation
matrices, or by applying the extended product rule for Kronecker products
below to the case of $n=2$.

#### Extended Kronecker Product Rule

$\displaystyle{\partial\left(A_{1}\otimes\ldots\otimes
A_{n}\right)\over\partial x^{\prime}}=$ $\displaystyle\left({\partial
A_{1}\over\partial
x^{\prime}}\otimes\bar{A}_{2:n}\right)\left(I_{s_{1}}\otimes
K_{n_{x},\prod_{j=2}^{n}s_{j}}\right)$ (7.158)
$\displaystyle\ldots+\sum_{k=2}^{n-1}\left(\bar{A}_{1:k-1}\otimes{\partial
A_{k}\over\partial
x^{\prime}}\otimes\bar{A}_{k+1:n}\right)\left(I_{\prod_{j=1}^{k}s_{j}}\otimes
K_{n_{x},\prod_{j=k+1}^{n}s_{j}}\right)$ (7.159)
$\displaystyle\ldots+A_{1:n-1}\otimes{\partial A_{n}\over\partial x^{\prime}}$
(7.160)

with $A_{k:l}:=A_{k}\otimes\cdots\otimes A_{l}$ and
$A_{i}\in\mathcal{M}_{r_{i}\times s_{i}}$ for all $i$.

###### Proof.

In the following $A_{i:j}$ will refer to the Kronecker product of all matrices
$A_{i}$ though $A_{j}$ with $i<j$. A typical element of $A_{1:n}$ is of the
form

$a^{(1)}_{i_{1},j_{1}}\cdot a^{(2)}_{i_{2},j_{2}}\cdot\ldots\cdot
a^{(n)}_{i_{n},j_{n}},$ (7.161)

where $\cdot$ is the multiplication operator acting on scalars and
$A_{k}=\left(a^{(k)}_{i_{k},j_{k}}\right)$. The definition of matrix
derivatives implies via

$A_{1:n}\otimes{\partial\over\partial x^{\prime}},$ (7.162)

that we can apply the product rule for scalars on each component individually.
Therefore we get

$\displaystyle{\partial(a^{(1)}_{i_{1},j_{1}}\cdot
a^{(2)}_{i_{2},j_{2}}\cdot\ldots\cdot a^{(n)}_{i_{n},j_{n}})\over\partial
x^{\prime}}$
$\displaystyle=\sum_{k=1}^{n}a^{(1)}_{i_{1},j_{1}}\cdot\ldots\cdot{\partial
a^{(k)}_{i_{k},j_{k}}\over\partial x^{\prime}}\cdot\ldots\cdot
a^{(n)}_{i_{n},j_{n}}.$ (7.163)

This representation implies that we can compute the derivate of $A_{1:n}$ by
computing $n$ individual term, where we keep all components constant, instead
for the matrix we want to the take the derivative of. Let $\bar{A_{\cdot}}$
denote a matrix which is to be treated as a constant when taking the
derivative.

${\partial\left(A_{1}\otimes\ldots\otimes A_{n}\right)\over\partial
x^{\prime}}={\partial(A_{1}\otimes\bar{A}_{2:n})\over\partial
x^{\prime}}+\sum_{k=2}^{n-1}{\partial\left(\bar{A}_{1:k-1}\otimes
A_{k}\otimes\bar{A}_{k+1:n-1}\right)\over\partial
x^{\prime}}+{\partial(\bar{A}_{1:n-1}\otimes A_{n})\over\partial x^{\prime}}.$
(7.164)

Using that, and by applying commutation results from Section 3, we can write

$\displaystyle{\partial\left(A_{1}\otimes\ldots\otimes
A_{n}\right)\over\partial x^{\prime}}=$
$\displaystyle{\partial(A_{1}\otimes\bar{A}_{2:n})\over\partial
x^{\prime}}+\sum_{k=2}^{n-1}{\partial\left(\bar{A}_{1:k-1}\otimes
A_{k}\otimes\bar{A}_{k+1:n}\right)\over\partial
x^{\prime}}+{\partial(\bar{A}_{1:n-1}\otimes A_{n})\over\partial x^{\prime}}.$
$\displaystyle=$
$\displaystyle\left(A_{1}\otimes\bar{A}_{2:n}\right)\otimes{\partial\over\partial
x^{\prime}}+\sum_{k=2}^{n-1}\left(\bar{A}_{1:k-1}\otimes
A_{k}\otimes\bar{A}_{k+1:n}\right)\otimes{\partial\over\partial
x^{\prime}}+\bar{A}_{1:n-1}\otimes{\partial A_{n}\over\partial x^{\prime}}$
$\displaystyle=$ $\displaystyle A_{1}\otimes\left({\partial\over\partial
x^{\prime}}\otimes\bar{A}_{2:n}\right)K_{d,\prod_{j=2}^{n}s_{j}}+\sum_{k=2}^{n-1}\bar{A}_{1:k-1}\otimes
A_{k}\otimes\left({\partial\over\partial
x^{\prime}}\otimes\bar{A}_{k+1:n}\right)K_{n_{x},\prod_{j=k+1}^{n}s_{j}}\ldots$
$\displaystyle\ldots+A_{1:n-1}\otimes{\partial A_{n}\over\partial x^{\prime}}$
$\displaystyle=$ $\displaystyle
A_{1}I_{s_{1}}\otimes\left({\partial\over\partial
x^{\prime}}\otimes\bar{A}_{2:n}\right)K_{d,\prod_{j=2}^{n}s_{j}}\ldots$
$\displaystyle\ldots+\sum_{k=2}^{n-1}\left(\bar{A}_{1:k-1}\otimes
A_{k}\right)I_{\prod_{j=1}^{k}s_{j}}\otimes\left({\partial\over\partial
x^{\prime}}\otimes\bar{A}_{k+1:n}\right)K_{n_{x},\prod_{j=k+1}^{n}s_{j}}\ldots$
$\displaystyle\ldots+A_{1:n-1}\otimes{\partial A_{n}\over\partial x^{\prime}}$
$\displaystyle=$ $\displaystyle\left({\partial A_{1}\over\partial
x^{\prime}}\otimes\bar{A}_{2:n}\right)\left(I_{s_{1}}\otimes
K_{n_{x},\prod_{j=2}^{n}s_{j}}\right)\ldots$
$\displaystyle\ldots+\sum_{k=2}^{n-1}\left(\bar{A}_{1:k-1}\otimes{\partial
A_{k}\over\partial
x^{\prime}}\otimes\bar{A}_{k+1:n}\right)\left(I_{\prod_{j=1}^{k}s_{j}}\otimes
K_{d,\prod_{j=k+1}^{n}s_{j}}\right)\ldots$ $\displaystyle+\ldots
A_{1:n-1}\otimes{\partial A_{n}\over\partial x^{\prime}}.$

∎

#### Chain Rule

${\partial\left(F\circ G\right)\over\partial x^{\prime}}={\partial
F\over\partial G^{\prime}}\left(I_{t}\otimes{\partial G\over\partial
x^{\prime}}\right).$ (7.165)

If $F$ is a vector of size $s\times 1$ this reduces to

${\partial\left(F\circ G\right)\over\partial x^{\prime}}={\partial
F\over\partial G^{\prime}}{\partial G\over\partial x^{\prime}}.$ (7.166)

The differentiation rules above can also be expressed on their differential
version:

#### Product Rule - Differential

$dFH={\partial FH\over\partial x^{\prime}}\left(I_{q}\otimes
dx\right)={\partial F\over\partial x^{\prime}}\left(H\otimes
dx\right)+F{\partial H\over\partial x^{\prime}}\left(I_{q}\otimes dx\right).$
(7.167)

#### Product Rule - Kronecker Product - Differential

$d\left(F\otimes H\right)={\partial\left(F\otimes H\right)\over\partial
x^{\prime}}\left(I_{tq}\otimes dx\right)=\left(F\otimes{\partial
H\over\partial x^{\prime}}+\left({\partial F\over\partial x^{\prime}}\otimes
H\right)\left(I_{t}\otimes K_{n_{x},q}\right)\right)\left(I_{tq}\otimes
dx\right)$ (7.168)

#### Chain Rule - Differential

$\displaystyle d\left(F\circ G\right)$ $\displaystyle={\partial\left(F\circ
G\right)\over\partial x^{\prime}}\left(I_{t}\otimes dx\right)$ (7.169)
$\displaystyle={\partial F\over\partial G^{\prime}}\left(I_{t}\otimes{\partial
G\over\partial x^{\prime}}\right)\left(I_{t}\otimes dx\right)$ (7.170)
$\displaystyle={\partial F\over\partial G^{\prime}}\left(I_{t}\otimes{\partial
G\over\partial x^{\prime}}dx\right).$ (7.171)

### B. Uniqueness and the Symmetrizer Matrix

For non-univariate function differentials the issue of non-uniqueness of the
corresponding matrix or vector of partial derivatives arises.

Given a matrix of partial derivatives of a vector function $f$ by $D_{f}$ we
can write the k-th order differential using

$d^{k}f=D_{f}\left(dx\right)^{\otimes k}.$ (7.172)

While the differential representation will be unique, the matrix of
derivatives $D_{f}$ might not be.111111See also Chacón and Duong (2021) for a
general discussion of this issue.

To see this, note the following example.

######  Example 2.

Consider the function

$\displaystyle f(x)$ $\displaystyle={1\over 2}a^{\prime}x^{\otimes 2}$ (7.173)
$\displaystyle={1\over
2}a_{1}x_{1}^{2}+\left(a_{2}+a_{3}\right)x_{1}x_{2}+{1\over 2}a_{4}x_{2}^{2}$
(7.174)

with $x=(x_{1},x_{2})^{\prime}$ and
$a=\left(a_{1},a_{2},a_{3},a_{4}\right)^{\prime}$. Writing down the second
order differential explicitly provides

$\displaystyle d^{2}f$
$\displaystyle=a_{1}dx_{1}^{2}+\left(a_{2}+a_{3}\right)dx_{1}dx_{2}+a_{2}dx_{2}^{2}$
(7.175) $\displaystyle=\left(a_{1},{a_{2}+a_{3}\over 2},{a_{2}+a_{3}\over
2},a_{4}\right)(dx)^{\otimes 2}.$ (7.176)

Using the second-order differential to identify the matrix of partial
derivatives then leads us to write

${\partial^{2}f\over\partial x^{\prime\otimes
2}}=\left(a_{1},{a_{2}+a_{3}\over 2},{a_{2}+a_{3}\over 2},a_{4}\right).$
(7.177)

This is also exactly the same vector we would have gotten if would have taken
partial derivatives of $f$ directly. However, we can write the second order
differential similarly by

$\displaystyle d^{2}f$
$\displaystyle=\left(a_{1},a_{2},a_{3},a_{4}\right)(dx)^{\otimes 2}$ (7.178)
$\displaystyle=a_{1}dx_{1}^{2}+\left(a_{2}+a_{3}\right)dx_{1}dx_{2}+a_{2}dx_{2}^{2},$
(7.179)

which would imply

${\partial^{2}f\over\partial x^{\prime\otimes
2}}=\left(a_{1},a_{2},a_{3},a_{4}\right).$ (7.180)

We can make the decision to choose a unique matrix of derivatives by requiring
that it fulfill some property of symmetry. Generally speaking, values of
$D_{f}$ corresponding to the same component of the differential should be
equal. For example, for a second order differential the values of $D_{f}$ that
are multiplied with $dx_{1}dx_{2}$ and $dx_{2}dx_{1}$ should be identical. The
symmetrizer matrix $S_{n,k}$ for vectors of size $n$ and Kronecker products
with $k$ factors introduced by Holmquist (1985) achieves exactly this.

######  Definition 5.

For the Kronecker product column vectors $x_{1}\otimes
x_{2}\otimes\cdots\otimes x_{m}$ with
$x_{1},x_{2},\ldots,x_{m}\in\mathbb{R}^{n}$, the symmetrizer matrix is given
by the following property

$S_{n,m}\left(x_{1}\otimes x_{2}\otimes\ldots\otimes x_{m}\right)={1\over
m!}\sum_{\sigma\in Sym(m)}x_{\sigma(1)}\otimes
x_{\sigma(2)}\otimes\ldots\otimes x_{\sigma(m)}.$ (7.181)

The symmetrizer matrix itself is symmetric. Therefore, for row vectors
$x_{1},x_{2},\ldots,x_{m}\in\mathbb{R}^{n}$ we get

$\left(x_{1}\otimes x_{2}\otimes\ldots\otimes x_{m}\right)S_{n,m}={1\over
m!}\sum_{\sigma\in Sym(m)}x_{\sigma(1)}\otimes
x_{\sigma(2)}\otimes\ldots\otimes x_{\sigma(m)}.$ (7.182)

Using this definition, it then follows that

$d^{k}f=D_{f}\left(dx\right)^{\otimes k}=D_{f}S_{n,k}\left(dx\right)^{\otimes
k},$ (7.183)

because $S_{n,k}\left(dx\right)^{\otimes k}=\left(dx\right)^{\otimes k}$.
While the differential does not change using the transformation, the matrix of
derivative is now given by $D_{f}S_{n,k}$.

######  Example 3.

In our example above, we can derive the symmetrizer matrix in the following
way: For vectors $x,y$ of size $d\times 1$ the symmetrizer fulfills the
following property

$S_{d,2}\left(x\otimes y\right)={x\otimes y+y\otimes x\over 2}$ (7.184)

for all vectors $x$ and $y$. We know that there exists a commutation matrix
$K_{d,d}$ such that $K_{d,d}\left(x\otimes y\right)=y\otimes x$. Therefore we
can write

$\displaystyle S_{d,2}\left(x\otimes y\right)$
$\displaystyle=I_{d^{2}}{x\otimes y\over 2}+K_{d,d}{x\otimes y\over 2}$
(7.185) $\displaystyle={1\over 2}\left(I_{d^{2}}+K_{d,d}\right)\left(x\otimes
y\right).$ (7.186)

Since this holds for arbitrary choices of $x$ and $y$ it follows that

$S_{d,2}={1\over 2}\left(I_{d^{2}}+K_{d,d}\right).$ (7.187)

For $d=2$ this becomes

$S_{2,2}=\begin{pmatrix}1\;0\;0\;0\\\ 0\;{1\over 2}\;{1\over 2}\;0\\\
0\;{1\over 2}\;{1\over 2}\;0\\\ 0\;0\;0\;1\end{pmatrix}={1\over
2}\left(I_{2}+K_{2,2}\right).$ (7.188)

Thusly, we can write

$\displaystyle d^{2}f$
$\displaystyle=\left(a_{1},a_{2},s_{3},a_{4}\right)S_{2,2}(dx)^{\otimes 2}$
(7.189) $\displaystyle=\left(a_{1},{a_{2}+a_{3}\over 2},{a_{2}+a_{3}\over
2},a_{4}\right)S_{2,2}(dx)^{\otimes 2}$ (7.190)
$\displaystyle={\partial^{2}f\over\partial x^{\prime\otimes 2}}(dx)^{\otimes
2}.$ (7.191)

#### Symmetrizing Faà Di Bruno

Similarly, for our generalization of Faà Di Bruno’s formula, this implies that
the differential representation

$\displaystyle
d^{n}f(g(x))=\sum_{k=1}^{n}f_{g^{k}}(g(x))\textbf{B}^{g}_{n,k}(dx)^{\otimes
n},$ (7.192)

implies a representation of the matrix of derivatives given by

${\partial^{n}(f\circ g)\over\partial(x^{\prime})^{\otimes
n}}=\sum_{k=1}^{n}f_{g^{k}}(g(x))\textbf{B}^{g}_{n,k}S_{n_{x},n}.$ (7.193)

instead of

${\partial^{n}(f\circ g)\over\partial(x^{\prime})^{\otimes
n}}=\sum_{k=1}^{n}f_{g^{k}}(g(x))\textbf{B}^{g}_{n,k}.$ (7.194)