# Random Tensors and their Normal Distributions

Changqing Xu Corresponding author. Email<EMAIL_ADDRESS>School of
Mathematics and Physics
Suzhou University of Science and Technology, Suzhou, China  Ziming Zhang
Email<EMAIL_ADDRESS>School of Mathematics and Physics
Suzhou University of Science and Technology, Suzhou, China

###### Abstract

The main purpose of this paper is to introduce the random tensor with normal
distribution, which promotes the matrix normal distribution to a higher order
case. Some basic knowledge on tensors are introduced before we focus on the
random tensors whose entries follow normal distribution. The random tensor
with standard normal distribution(SND) is introduced as an extension of random
normal matrices. As a random multi-array deduced from an affine transformation
on a SND tensor, the general normal random tensor is initialised in the paper.
We then investigate some equivalent definitions of a normal tensor and present
the description of the density function, characteristic function, moments, and
some other functions related to a random matrix. A general form of an even-
order multi-variance tensor is also introduced to tackle a random tensor. Also
presented are some equivalent definitions for the tensor normal distribution.
We initialize the definition of high order standard Gaussian tensors, general
Gaussian tensors, deduce some properties and their characteristic functions.

keywords: Random matrix; random tensor; Gaussian matrix; Gaussian tensor;
Characteristic function.
AMS Subject Classification: 53A45, 15A69.

##  1 Introduction

The systematic treatment of multivariate statistics through matrix theory,
treating the multivariate statistics with a compact way, has been developed
since 1970s [5, 29]. In [16] Kollo and Rosen introduce the matrix normal
distribution as the extension of the ordinary multivariate normal one. In
2013, Ohlson, Ahmad and Rosen [23] introduce the multilinear normal
distribution for the analysis of multiway data as an extension of the matrix-
variate normal distribution, but there the main tool is still limited to
matrix theory. Also in 2013, Manceur and Dutilleul[19] investigate the problem
of the maximum likelihood estimation (MLE) for the tensor normal distribution.
Actually early in 2003, P. J. Basser and S. Pajevic [6] ever proposed 4-order
tensor-form normal distribution to preserve and describe the algebraic form
and geometric structure of the tensor random variable and perform various
algebraic and geometric operations on it.
In multivariate statistics, the $k$-moment of a random vector
$\textbf{x}\in\mathbb{R}^{n}$ is conventionally depicted by a matrix $M_{k}$
of size $n\times n^{k-1}$ for any positive integer $k$. The _multivariate
normal distribution_ , also called the _joint normal distribution_ , usually
deals with random vectors with normal distributions. The study of random
matrices, motivated by quantum mechanics and wireless communications [10, 11,
21, 22] etc. in the past 70 years, mainly focuses on the spectrum
properties[22] and can be used in many areas such as the classical analysis
and number theory using enumerative combinatorics[22], Fredholm
determinants[20], diffusion processes[1], integrable systems[2], the Riemann
Hilbert problem[13] and the free probability theory in the 1990s. The theory
of random matrices in multivariate statistics basically focuses on the
distribution of the eigenvalues of the matrices[28] .
While there have already been some well-established work on tensor-form
multivariate random distributions as mentioned above, we noticed that all
those work still stays in matrix scope. For example, in [23] the authors
prefer the Kronecker product instead of tensor product and the commutation
matrix instead of the commutation tensor which may make things much simple, as
we will describe in the paper.
The need to use the high order tensors has manifested in many areas more than
half a century ago, and the recent growing development of multivariate
distribution theory poses new challenge for finding some novel tools to
describe classical statistical concepts e.g. moment, characteristic function
and covariance etc. This in turn has facilitated the development of higher
order tensor theory for multilinear regression model[8] and the higher order
derivatives of distribution functions. Meanwhile, the description of an
implicit multi-relationship among a family of random variables pose a
challenge to modern statisticians. The applications of the high order tensors
in statistics was initialized by Cook etc. [8, 9] when the envelope models
were established.
In this paper, we first use tensors to express the high order derivatives,
which in turn leads to the simplification of the high-order moments and the
covariances of a random matrix. We also introduce the normal distributions of
a random matrix as well as that of a random tensor. The Gaussian tensors are
investigated to extend the random matrix theory.
By a random vector $\textbf{x}:=(x_{1},\ldots,x_{n})^{\top}\in\mathbb{R}^{n}$
we mean that each component $x_{i}$ is a random variable (r.v.). Here we
usually do not distinguish a row and a column vector unless specifically
mentioned. There are several equivalent definitions for a random vector to be
Gaussian. Given a constant vector $\mu\in\mathbb{R}^{n}$ and a positive
semidefinite matrix $\Sigma\in\mathbb{R}^{n\times n}$. A random vector
$\textbf{x}\in\mathbb{R}^{n}$ is called a _Gaussian_ or _normal_ vector with
parameter $(\mu,\Sigma)$ if it is normally distributed with
$E[\textbf{x}]=\mu$ and $\texttt{Var}(\textbf{x})=\Sigma$. This is equivalent
to a single variable normal distribution of $\alpha^{\top}\textbf{x}$ for all
$\alpha\in\mathbb{R}^{n}$. It is obvious from this fact that each component of
a normal vector is normal. The converse is nevertheless not true. A random
tensor $\mathcal{A}=(A_{i_{1}i_{2}\ldots i_{m}})$ is an $m$-order tensor whose
entries are random variables. A random matrix is a matrix whose entries are
random variables.
The covariance matrix of a random vector x restores the variances and
covariances of its coordinates and plays a very important role in the
statistical analysis. However, it cannot demonstrate the multi-variances of a
group of variables. On the other hand, the $k$-moment of a random vector
$\textbf{x}\in\mathbb{R}^{n}$ can be defined through the high-order derivative
of its characteristic function, which can be depicted much more clearly by a
high-order tensor. Note that an $m$-order derivative of a function
$f(\textbf{x})=f(x_{1},x_{2},\ldots,x_{n})$, usually defined as an $n\times
n^{m-1}$ matrix, can be expressed as a symmetric $m$-order $n$-dimensional
tensor. A direct bonus of a tensor form of derivative is that we can locate
and label easily any entry of $\frac{\partial^{m}f}{\partial\textbf{x}^{m}}$.
We denote $[n]$ for the set $\left\\{1,2,\ldots,n\right\\}$ for any positive
integer $n$ and $\mathbb{R}$ the field of real numbers. A tensor (also called
a _multiway matrix_) $\mathcal{A}$ of size $\rm{I}:=d_{1}\times
d_{2}\times\ldots\times d_{m}$ is an _$m$ -array_. $\mathcal{A}$ is called an
$m$th order $n$-dimensional real tensor or an $m\times n$ tensor if
$d_{1}=\ldots=d_{m}=[n]$. Let $\mathcal{T}({\rm I})$ be the set of all the
$m$-order tensors indexed by $\rm{I}$, $\mathcal{T}_{m;n}$ be the set of all
$m$th order $n$-dimensional real tensors , and $\mathcal{T}_{m}$ be the set of
all $m$-order tensors. An $m$-order tensor $\mathcal{A}\in\mathcal{T}_{m;n}$
is called a _symmetric tensor_ if the value of each entry $A_{i_{1}i_{2}\ldots
i_{m}}$ is invariant under any permutation on its indices. A tensor
$\mathcal{A}=(A_{i_{1}i_{2}\ldots i_{m}})\in\mathcal{T}({\rm I})$ can be
regarded as a multiway array. Thus a scalar, a vector and a matrix is
respectively a tensor of order zero, one and two. A tensor
$\mathcal{A}\in\mathcal{T}_{m;n}$ is called a _symmetric tensor_ if each entry
of $\mathcal{A}$ is invariant for any permutation of its indices. An $m$th
order $n$-dimensional real tensor $\mathcal{A}\in\mathcal{T}_{m;n}$ is
associated with an $m$-order $n$-variate homogeneous polynomial in the
following

$f_{\mathcal{A}}(\textbf{x}):=\mathcal{A}\textbf{x}^{m}=\sum\limits_{i_{1},i_{2},\ldots,i_{m}}A_{i_{1}i_{2}\ldots
i_{m}}x_{i_{1}}x_{i_{2}}\ldots x_{i_{m}}$ (1.1)

We denote by $\mathcal{S\\!T}_{m;n}$ the set of all $m$th order
$n$-dimensional symmetric tensors . A symmetric tensor
$\mathcal{A}\in\mathcal{S\\!T}_{m;n}$ is called _positive semidefinite_
(acronym by PSD) if $f_{\mathcal{A}}(\textbf{x}):=\mathcal{A}\textbf{x}^{m}>0$
($\geq 0$) for all $0\neq\textbf{x}\in\mathbb{R}^{n}$. $\mathcal{A}$ is called
a copositive tensor if $f_{\mathcal{A}}(\textbf{x})\geq 0$ for all nonnegative
vector x. For the study of the symmetric tensors, PSD tensors, the copositive
tensors, including their spectrum, decompositions and other properties, we
refer the reader to [4, 7, 24, 25] and [7].
Let $I:=I_{1}\times\ldots\times I_{m}$ where each $I_{k}$ represents an index
set (usually a set in form $\left\\{1,2,\ldots,n_{k}\right\\}$) and
$\mathcal{A}\in\mathcal{T}({\rm I})$. For any $k\in[m]$ and $j\in I_{k}$,
$\mathcal{A}$’s $j$-slice along the $k$-mode (denoted by $A^{(k)}[j]$) is
defined as an $(m-1)-$order tensor $(A_{i_{1}\ldots i_{k-1}ji_{k+1}\ldots
i_{m}})$ where the $k$th subscript of each entry of $\mathcal{A}$ is fixed to
be $j\in[I_{k}]$. Thus by slicing an $m$-order tensor can be made into an
$(m-1)$-order tensor. A high order tensor can be flattened or unfolded into a
matrix by slicing iteratively. For example, an $m\times n\times p$ tensor
$\mathcal{A}$ can be unfolded to be a matrix $A[1]\in\mathbb{R}^{m\times np}$
along the first mode and $A[2]\in\mathbb{R}^{n\times pm}$ if along the second
mode. There are ten options to flatten a 4-order tensor $m\times n\times
p\times q$ into a matrix: four to reserve one mode and stack the other three
and six to group two modes together to form a matrix. The product of tensors
can be defined in many different ways. Given any two tensors
$\mathcal{A},\mathcal{B}$ with appropriate size, we usually choose one (or
several directions), say $k$, and define their product w.r.t. the chosen
mode(s). This can be regarded as a generalisation of the matrix product. For
more detail, we refer the reader to [27] and [3].
Given a random vector $\textbf{x}\in\mathbb{R}^{n}$. The characteristic
function (CF) of x, is defined by

$\phi_{\textbf{x}}(\textbf{t})=E[exp(\imath\textbf{t}^{\prime}\textbf{x})],\forall\textbf{t}\in\mathbb{R}^{n}$

For any positive integer $k$ , the $k$-_moment_ of a random vector x is
defined by

$m_{k}(\textbf{x})=\frac{1}{\imath^{k}}\frac{d^{k}}{d\textbf{t}^{k}}\phi_{\textbf{x}}(\textbf{t})|_{\textbf{t}=0}$
(1.2)

Note that
$m_{1}(\textbf{x})=\frac{d}{d\textbf{t}}\phi_{\textbf{x}}(\textbf{t})\in\mathbb{R}^{n}$
is a vector and
$m_{k}(\textbf{x})=\frac{d^{k}}{d\textbf{t}^{k}}\phi_{\textbf{x}}(\textbf{t})\in\mathbb{R}^{n\times
n^{k-1}}$ is a matrix for each $k>1$ by the conventional definition.
To simplify the definition of the $k$-moments of a random vector, we present a
tensor form of the high order multivariate function. Let
$\textbf{y}=f(\textbf{x})$ be a mapping from $\mathcal{C}^{n}$ to
$\mathcal{C}^{m}$, i.e.,
$\textbf{y}=(y_{1},\ldots,y_{m})^{\top}\in\mathcal{C}^{m}$ with each
components $y_{i}=f_{i}(\textbf{x})$ sufficiently differentiable. Then we
denote $H(\textbf{y},\textbf{x})=(h_{ij})$ as the Jacobi matrix of y w.r.t. x
defined by $h_{ij}:=\frac{dy_{i}}{dx_{j}}$ for all $i\in[m],j\in[n]$. Thus
$H(\textbf{y},\textbf{x})\in\mathbb{R}^{m\times n}$. Now we define the
$k$-order differentiation tensor by

$\mathcal{H}^{k}(\textbf{y},\textbf{x})=(h_{ij_{1}j_{2}\ldots j_{k}})$ (1.3)

which is an $(k+1)$-order tensor of size
$m\times\overbrace{n\times\ldots\times n}^{k}$ where

$h_{ij_{1}j_{2}\ldots j_{k}}=\frac{\partial^{k}y_{i}}{\partial
x_{j_{1}}\partial x_{j_{2}}\ldots\partial x_{j_{k}}}$ (1.4)

for any $i\in[m],j_{1},\ldots,j_{k}\in[n]$. Recall that the conventional form
for the $k$-order differentiation of a mapping $\textbf{y}=f(\textbf{x})$
produces an $m\times n^{k}$ matrix $H^{k}(\textbf{y},\textbf{x})=(h_{ij})$
where $h_{ij}=\frac{\partial^{k}y_{i}}{\partial x_{j_{1}}\partial
x_{j_{2}}\ldots\partial x_{j_{k}}}$ with $j=\sum\limits_{s=1}^{k}j_{s}n^{k-s}$
($j=1,2,\ldots,n^{k}$), making the location of each entry ambiguous.
In this paper, we use tensor form to simplify high order differentiations. The
$k$-moment of a random vector $\textbf{x}\in\mathbb{R}^{n}$ is defined by

$m_{k}[\textbf{x}]=E[\textbf{x}^{k}]=E[\overbrace{\textbf{x}\times\textbf{x}\times\ldots\times\textbf{x}}^{k}]$
(1.5)

This is a natural extension of the $k$-moment in the univariate case since

$m_{1}[\textbf{x}]=E[\textbf{x}]\in\mathbb{R}^{n},m_{2}[\textbf{x}]=E[\textbf{x}\textbf{x}^{\top}]\in\mathbb{R}^{n\times
n},\cdots,m_{k}[\textbf{x}]=E[\textbf{x}^{k}]\in\mathcal{T}_{m;n}.$

The definition is identical to the one through characteristic function in the
tensor form, as in the following.

######  Lemma 1.1.

Let $k$ be any positive integer and $\textbf{x}\in\mathbb{R}^{n}$ be a random
vector with characteristic function $\phi_{\textbf{x}}(\textbf{t})$ (with
$\textbf{t}\in\mathbb{R}^{n}$). Then

$E[\textbf{x}^{k}]=\frac{1}{\imath^{k}}\frac{d^{k}}{d\textbf{t}^{k}}\phi_{\textbf{x}}(\textbf{t})|_{\textbf{t}=0}$
(1.6)

where $0$ is a zero vector in $\mathbb{R}^{n}$.

###### Proof.

It is easy to see from the definition of the characteristic function
$\phi_{\textbf{x}}(\textbf{t})$ and (1.3) and (1.4) that

$\frac{d^{k}}{d\textbf{t}^{k}}\phi_{\textbf{x}}(\textbf{t})=\imath^{k}E[\exp\\{\imath\textbf{t}^{\top}\textbf{x}\\}\textbf{x}^{k}]$
(1.7)

Thus (1.6) holds. ∎

Similarly we can also simplify the definition of the $k$-_central moment_
$\bar{m}_{k}[\textbf{x}]$ by
$\bar{m}_{k}[\textbf{x}]=E[(\textbf{x}-E[\textbf{x}])^{k}]$. Note that our
definition is consistent with the traditional one for $k\leq 2$. In the
following section, we will extend the $k$-moment, the characteristic function,
and the related terminology to the case for the random matrices.

##  2 The tensor forms of derivatives of matrices

Given two matrices $A\in\mathbb{R}^{m_{1}\times
n_{1}},B\in\mathbb{R}^{m_{2}\times n_{2}}$. We may form an 4-order tensor by
the tensor product (or outer product) of $A,B$. For any subset
$\theta:=\left\\{s,t\right\\}\subset\left\\{1,2,3,4\right\\}$ with $s<t$, we
denote $\left\\{p,q\right\\}=\left\\{s,t\right\\}^{c}$ be the complement of
$\theta$ with $p<q$. We define the tensor product of $A,B$ along $\theta$,
denoted $A\times_{(s,t)}B$, as the 4-order tensor where

$(A\times_{(s,t)}B)_{i_{1}i_{2}i_{3}i_{4}}=A_{i_{p}i_{q}}B_{i_{s}i_{t}}$

For example,
$(A\times_{(1,2)}B)_{i_{1}i_{2}i_{3}i_{4}}=A_{i_{3}i_{4}}B_{i_{1}i_{2}}$,
$(A\times_{(3,4)}B)_{i_{1}i_{2}i_{3}i_{4}}=A_{i_{1}i_{2}}B_{i_{3}i_{4}}$, and
$(A\times_{(2,3)}B)_{i_{1}i_{2}i_{3}i_{4}}=A_{i_{1}i_{4}}B_{i_{2}i_{3}}$.
Similarly we can also define
$A\times_{(1,3)}B,A\times_{(1,4)}B,A\times_{2,4)}B$. When
$\left\\{(s,t)\right\\}=\left\\{(3,4)\right\\}$, we just write $A\times B$
instead of $A\times_{(3,4)}B$. Note that

$A\times B\in\mathbb{R}^{m_{1}\times n_{1}\times m_{2}\times
n_{2}},A\times_{(2,3)}B\in\mathbb{R}^{m_{1}\times m_{2}\times n_{2}\times
n_{1}}.$

For a 4-order tensor $\mathcal{A}\in\mathbb{R}^{m\times n\times p\times q}$
and a matrix $B\in\mathbb{R}^{p\times q}$, the product $C=\mathcal{A}\times B$
is usually referred to $\mathcal{A}\times_{(3,4)}B$ which, sometimes written
simply as $\mathcal{A}B$, is a matrix $C\in\mathbb{R}^{m\times n}$ where
$C_{ij}=\sum\limits_{i^{\prime},j^{\prime}}A_{iji^{\prime}j^{\prime}}B_{i^{\prime}j^{\prime}}$.
Note that the product $\mathcal{A}\times_{4}B$ is still a 4-order tensor since

$(\mathcal{A}\times_{4}B)_{i_{1}i_{2}i_{3}i_{4}}=\sum\limits_{k}A_{i_{1}i_{2}i_{3}k}B_{ki_{4}}.$

For $A\in\mathbb{R}^{m\times n}$, denote $A^{[2]}:=A\times A$ and
$A^{(2)}=A\times_{(2,4)}A$, i.e.,

$A^{[2]}_{i_{1}i_{2}i_{3}i_{4}}=A_{i_{1}i_{2}}A_{i_{3}i_{4}},\quad
A^{(2)}_{i_{1}i_{2}i_{3}i_{4}}=A_{i_{1}i_{3}}A_{i_{2}i_{4}}.$

It is obvious that $A^{[2]}$ has size $m\times n\times m\times n$ and
$A^{(2)}$ has size $m\times m\times n\times n$.

######  Proposition 2.1.

Let $A,B,C,D$ be any matrices of appropriate sizes and
$\left\\{s,t\right\\}\subset[4]$ and
$\left\\{p,q\right\\}=\left\\{s,t\right\\}^{c}$. Then

(1)

$(A\times_{(s,t)}B)\times_{(s,t)}C=(B,C)A$ where $B,C$ are of same size and
$(X,Y)$ stands for the inner product of $X,Y$.

(2)

$(A\times_{(s,t)}B)\times_{(p,q)}C=A\times_{(s,t)}(B\times_{(p,q)}C)$.

(3)

$(A\times B)\times_{4}C=A\times(BC)$.

(4)

$(A\times_{(s,t)}B)(C\times_{(s,t)}D)=(AC)\times_{(s,t)}(B\times_{(s,t)}D)$.

###### Proof.

(1). We may assume that $A\in\mathbb{R}^{m\times n},B,C\in\mathbb{R}^{p\times
q}$ (the equation is valid only if $B,C$ are of the same size). For
simplicity, we let $\left\\{s,t\right\\}=\left\\{3,4\right\\}$ and denote
$D=(A\times_{(3,4)}B)\times_{(3,4)}C$. Then for any pair $(i,j)$, we have by
definition

$\displaystyle D_{ij}$ $\displaystyle=$
$\displaystyle\sum\limits_{i^{\prime},j^{\prime}}(A\times_{(3,4)}B)_{iji^{\prime}j^{\prime}}C_{i^{\prime}j^{\prime}}$
$\displaystyle=$
$\displaystyle\sum\limits_{i^{\prime},j^{\prime}}A_{ij}B_{i^{\prime}j^{\prime}}C_{i^{\prime}j^{\prime}}=(B,C)A_{ij}$

which immediately implies (1) for $(s,t)=(3,4)$. Similarly we can also show
its validity for other cases. The second and the third item can also be
checked using the same technique. To show the last item, we consider the case
$(s,t)=(2,4)$ and rewrite it in form

$(A_{1}\times A_{2})(B_{1}\times B_{2})=(A_{1}B_{1})\times(A_{2}B_{2})$ (2.1)

where $\times:=\times_{(2,4)}$ and $A_{i}\in\mathbb{R}^{m_{i}\times
n_{i}},B_{i}\in\mathbb{R}^{n_{i}\times p_{i}}$ for $i=1,2$. Denote the tensor
of the left hand side and the right hand side resp. by $L$ and $R$. Then we
have

$\displaystyle R_{i_{1}i_{2}i_{3}i_{4}}$ $\displaystyle=$
$\displaystyle(A_{1}B_{1})_{i_{1}i_{3}}(A_{2}B_{2})_{i_{2}i_{4}}$
$\displaystyle=$
$\displaystyle(\sum\limits_{j=1}^{n_{1}}a_{i_{1}j}^{(1)}b_{ji_{3}}^{(1)})(\sum\limits_{k=1}^{n_{2}}a_{i_{2}k}^{(2)}b_{ki_{4}}^{(2)})$
$\displaystyle=$
$\displaystyle\sum\limits_{j,k}a_{i_{1}j}^{(1)}a_{i_{2}k}^{(2)}b_{ji_{3}}^{(1)}b_{ki_{4}}^{(2)}$
$\displaystyle=$ $\displaystyle\sum\limits_{j,k}(A_{1}\times
A_{2})_{i_{1}i_{2}jk}(B_{1}\times
B_{2})_{jki_{3}i_{4}}=L_{i_{1}i_{2}i_{3}i_{4}}$

for all possible $i_{1},i_{2},i_{3},i_{4}$. Thus (2.1) holds. This argument
can also be extended to other cases. ∎

Proposition 2.1 can be extended to multiplications of tensors, which allows us
to define $A^{[k]}$ recursively as $A^{[k+1]}=A^{[k]}\times_{(2k+1,2k+2)}A$
and $A^{(k+1)}=A^{(k)}\times_{(2k,2k+2)}A$. Note that the size of $A^{[k]}$ is
$(m\times n)^{[k]}$ and the size of $A^{(k)}$ is $m^{[k]}\times n^{[k]}$.
Recall that the product of a 4-order tensor $\mathcal{A}\in\mathbb{R}^{m\times
n\times m\times n}$ and a matrix $P\in\mathbb{R}^{m\times n}$ is defined as

$(\mathcal{A}P)_{ij}=\sum\limits_{i^{\prime},j^{\prime}}A_{iji^{\prime}j^{\prime}}P_{i^{\prime}j^{\prime}},\quad(P\mathcal{A})_{ij}=\sum\limits_{i^{\prime},j^{\prime}}P_{i^{\prime}j^{\prime}}A_{i^{\prime}j^{\prime}ij}$

######  Corollary 2.2.

(1)

$(A\times_{(s,t)}I_{n})\times_{(s,t)}I_{n}=nA$.

(2)

$(I_{m}\times I_{n})A=\texttt{Tr}(A)I_{m}$ for any $A\in\mathbb{R}^{n\times
n}$.

(3)

$A\times(I_{m}\times_{c}I_{n})=(I_{m}\times_{c}I_{n})A=A$ for any
$A\in\mathbb{R}^{m\times n}$.

(4)

$A^{\top}(I_{n}\times_{(2,3)}I_{m})=A,(I_{m}\times_{(2,3)}I_{n})A^{\top}=A$
for any $A\in\mathbb{R}^{m\times n}$.

Recall that the commutation tensor $\mathcal{K}_{n,m}=(K_{ijkl})$ is an
$n\times m\times m\times n$ tensor which transforms a matrix
$A\in\mathbb{R}^{m\times n}$ into its transpose [30], i.e.,
$\mathcal{K}_{m,n}A^{\top}=A$. It is easy to see that
$\mathcal{K}_{m,n}=I_{m}\times_{(2,3)}I_{n}$. Note also that tensor
$I_{m}\times_{c}I_{n}$ can be regarded as the identity tensor in the linear
space $\mathbb{R}^{m\times n\times m\times n}$ by (3) of Corollary 2.2.
Let $X=(x_{ij})\in\mathbb{R}^{m\times n}$ be a random matrix whose entries are
independent variables. Let $Y=(y_{ij})\in\mathbb{R}^{p\times q}$ be a matrix
each of whose entries $y_{ij}$ can be regarded as a function of $X$. The
derivative $\frac{dY}{dX}$ is interpreted as the 4-order tensor
$A=(A_{i_{1}i_{2}i_{3}i_{4}})$ of size $m\times n\times p\times q$ whose
entries are defined by

$A_{i_{1}i_{2}i_{3}i_{4}}=\frac{dY_{i_{3}i_{4}}}{dX_{i_{1}i_{2}}}.$

In order to simplify the expressions of high order moments of random matrices,
we now use tensors to describe the derivatives of matrices. In the following,
we will present some known results in tensor forms other than in the
conventional matrix versions. The following lemma is the derivative chain rule
in the matrix version.

######  Lemma 2.3.

Let $X\in\mathbb{R}^{m_{1}\times n_{1}},Y\in\mathbb{R}^{m_{2}\times
n_{2}},Z\in\mathbb{R}^{m_{3}\times n_{3}}$, and $Z=Z(Y),Y=Y(X)$. Then we have

$\frac{dZ}{dX}=\frac{dY}{dX}\times\frac{dZ}{dY}$ (2.2)

###### Proof.

We denote $\mathcal{A}=\frac{dY}{dX},\mathcal{B}=\frac{dZ}{dY}$ and
$\mathcal{C}=\frac{dZ}{dX}$. By definition we have

$\mathcal{A}\in\mathbb{R}^{m_{1}\times n_{1}\times m_{2}\times
n_{2}},\mathcal{B}\in\mathbb{R}^{m_{2}\times n_{2}\times m_{3}\times
n_{3}},\mathcal{C}\in\mathbb{R}^{m_{1}\times n_{1}\times m_{3}\times n_{3}}.$

Then for any given
$(i_{1},i_{2},i_{3},i_{4})\in[m_{1}]\times[n_{1}]\times[m_{3}]\times[n_{3}]$,
we have

$\displaystyle A_{i_{1}i_{2}i_{3}i_{4}}$ $\displaystyle=$
$\displaystyle(\frac{dZ}{dX})_{i_{1}i_{2}i_{3}i_{4}}=\frac{dz_{i_{3}i_{4}}}{dx_{i_{1}i_{2}}}$
$\displaystyle=$
$\displaystyle\sum\limits_{j_{1},j_{2}}\frac{dy_{j_{1}j_{2}}}{dx_{i_{1}i_{2}}}\frac{dz_{i_{3}i_{4}}}{dy_{j_{1}j_{2}}}$
$\displaystyle=$
$\displaystyle(\frac{dY}{dX}\times\frac{dZ}{dY})_{i_{1}i_{2}i_{3}i_{4}}$

Thus (2.2) holds. ∎

Lemma 2.3 can be extended to a more general case:

######  Lemma 2.4.

(1)

Let $Z=Z(Y_{1},Y_{2},\ldots,Y_{n})$ be the matrix-valued function of
$Y_{1},Y_{2},\ldots,Y_{n}$ where $Y_{k}=Y_{k}(X)$ for all $k\in[n]$. Then

$\frac{dZ}{dX}=\sum\limits_{k=1}^{n}\frac{dY_{k}}{dX}\times\frac{dZ}{dY_{k}}$
(2.3)

(2)

Let $X,Y,Z,U$ be matrix forms of variables and $U=U(Z),Z=Z(Y),Y=Y(X)$. Then we
have the chain

$dU=dX\times(dY/dX)\times(dZ/dY)\times(dU/dZ)$ (2.4)

where $dX=(dx_{ij})$.

It is easy to verify the results in Lemma 2.4 so that we omit its proof here.
The following results on the matrix derivatives are useful and will be used in
the next section.

######  Theorem 2.5.

Let $X=(X_{ij})\in\mathbb{R}^{m\times n}$ be a matrix whose elements are
independent variables. Then

(1)

$\frac{dX}{dX}=I_{m}\times_{c}I_{n}$.

(2)

$\frac{dX^{\top}}{dX}=I_{m}\times_{(2,3)}I_{n}=\mathcal{K}_{m,n}$.

(3)

$\frac{d(YZ)}{dX}=\frac{dY}{dX}\times_{4}Z+\frac{dZ}{dX}\times_{3}Y^{\top}$.

(4)

$\frac{dX^{2}}{dX}=I_{n}\times_{c}X+X^{\top}\times_{c}I_{n}$ when $m=n$.

(5)

$\frac{dX^{k}}{dX}=\sum\limits_{p=0}^{k-1}\left[(X^{\top})^{p}\times_{c}X^{k-1-p}\right]$
where $X\in\mathbb{R}^{n\times n}$.

###### Proof.

Let $\mathcal{A}=\frac{dX}{dX}$. Then we have by definition and the
independency of the elements of $X$ that

$A_{i_{1}i_{2}i_{3}i_{4}}=\frac{dx_{i_{3}i_{4}}}{dx_{i_{1}i_{2}}}=\delta_{i_{1}i_{3}}\delta_{i_{2}i_{4}}=(I_{m}\times_{c}I_{n})_{i_{1}i_{2}i_{3}i_{4}}$

Thus (1) is proved. Similarly we can prove (2) by noticing that

$\left[\frac{dX^{\top}}{dX}\right]_{i_{1}i_{2}i_{3}i_{4}}=\delta_{i_{1}i_{4}}\delta_{i_{2}i_{3}}$

which implies

$\frac{dX^{\top}}{dX}=I_{m}\times_{(2,3)}I_{n}=\mathcal{K}_{m,n}.$

To prove (3), we let $Y\in\mathbb{R}^{p\times r},Z\in\mathbb{R}^{r\times q}$.
Then $\frac{d(YZ)}{dX}\in\mathbb{R}^{m\times n\times p\times q}$ whose
elements are

$\displaystyle(\frac{d(YZ)}{dX})_{i_{1}i_{2}i_{3}i_{4}}$ $\displaystyle=$
$\displaystyle\frac{d[(YZ)_{i_{3}i_{4}}]}{dX_{i_{1}i_{2}}}=\frac{d[\sum\limits_{k}y_{i_{3}k}z_{ki_{4}}]}{dx_{i_{1}i_{2}}}=\sum\limits_{k}\frac{d(y_{i_{3}k}z_{ki_{4}})}{dx_{i_{1}i_{2}}}$
$\displaystyle=$
$\displaystyle\sum\limits_{k}\left[\frac{d(y_{i_{3}k})}{dx_{i_{1}i_{2}}}z_{ki_{4}}+y_{i_{3}k}\frac{d(z_{ki_{4}})}{dx_{i_{1}i_{2}}}\right]$
$\displaystyle=$
$\displaystyle(\frac{dY}{dX}\times_{4}Z)_{i_{1}i_{2}i_{3}i_{4}}+(\frac{dZ}{dX}\times_{3}Y^{\top})_{i_{1}i_{2}i_{3}i_{4}}$

To prove (4), we let $X\in\mathbb{R}^{n\times n}$ and take $Y=Z=X$. By (3) and
(1), we have

$\displaystyle\frac{dX^{2}}{dX}$ $\displaystyle=$
$\displaystyle\frac{dX}{dX}\times_{4}X+X\times_{3}\frac{dX}{dX}$
$\displaystyle=$
$\displaystyle(I_{n}\times_{c}I_{n})\times_{4}X+X\times_{3}(I_{n}\times_{c}I_{n})$
$\displaystyle=$ $\displaystyle I_{n}\times_{c}X+X^{\top}\times_{c}I_{n}$

To prove (5), we use the induction method to $k$. For $k=1$, the result is
immediate since both sides of (5) are identical to $I_{n}\times_{c}I_{n}$ by
(1). The result is also valid for $k=2$ by (4). Now suppose it is valid for a
positive integer $k>2$. We come to show its validity for $k+1$. By (3) we have

$\displaystyle\frac{dX^{k+1}}{dX}$ $\displaystyle=$
$\displaystyle\frac{dX^{k}}{dX}\times_{4}X+(X^{\top})^{k}\times_{3}\frac{dX}{dX}$
$\displaystyle=$
$\displaystyle\sum\limits_{p=0}^{k-1}\left[(X^{\top})^{p}\times_{c}X^{k-p}\right]+(I_{n}\times_{c}I_{n})\times_{3}(X^{k})^{\top}$
$\displaystyle=$
$\displaystyle\sum\limits_{p=0}^{k-1}\left[(X^{\top})^{p}\times_{c}X^{k-p}\right]+(X^{\top})^{k}\times
I_{n}$ $\displaystyle=$
$\displaystyle\sum\limits_{p=0}^{k}\left[(X^{\top})^{p}\times X^{k-p}\right]$

Thus we complete the proof of (5). ∎

######  Corollary 2.6.

Let $X=(X_{ij})\in\mathbb{R}^{m\times n}$ be a matrix whose elements are
independent, and $A\in\mathbb{R}^{p\times m},B\in\mathbb{R}^{n\times q}$ be
the constant matrices. Then

(1)

$\frac{d(AXB)}{dX}=A^{\top}\times_{c}B$.

(2)

$\frac{d(\det(X))}{dX}=\det(X)X^{-\top}$.

(3)

$\frac{d(\texttt{Tr}(X))}{dX}=I_{n}$ for $X\in\mathbb{R}^{n\times n}$.

(4)

$\frac{dX^{-1}}{dX}=-X^{-\top}\times X^{-1}$ when $X\in\mathbb{R}^{n\times n}$
is invertible.

###### Proof.

To prove (1), we take $Y=A,Z=XB$. By (4) of Theorem 2.5, we get

$\frac{d(AXB)}{dX}=\frac{d(XB)}{dX}\times_{3}A^{\top}=A^{\top}\times_{3}(\frac{d(X)}{dX}\times_{4}B)=A^{\top}\times_{3}[(I_{m}\times
I_{n})\times_{4}B]=A^{\top}\times B$ (2.5)

To prove (2), we denote $A=\frac{d(detX)}{dX}=(A_{ij})$. Then for any given
pair $(i,j)\in[m]\times[n]$, we have by the expansion of the determinant

$\displaystyle A_{ij}$ $\displaystyle=$
$\displaystyle\frac{d(\det(X))}{X_{ij}}=\frac{d}{dX_{ij}}(\sum\limits_{k=1}^{n}(-1)^{i+k}X_{ik}\det(X(i|k)))$
$\displaystyle=$ $\displaystyle(-1)^{i+j}\det(X(i|j))=[\det(X)X^{-1}]_{ji}$

where $X(i|j)$ represents the submatrix of $X$ obtained by the removal of the
$i$th row and the $j$th column of $X$. Note that Thus we have
$\frac{d(\det(X))}{dX}=\det(X)X^{-\top}$. Now (3) can be verified by noticing
the fact that for all $(i,j)$

$[\frac{d(\texttt{Tr}(X))}{dX}]_{ij}=\sum\limits_{k=1}^{n}\frac{d(X_{kk})}{dX_{ij}}=\sum\limits_{k=1}^{n}\delta_{ik}\delta_{jk}=(e_{i},e_{j})=\delta_{ij}$

where $e_{i}\in\mathbb{R}^{n}$ is the $i$th row of the identity matrix
$I_{n}$. Now we prove (4). Using (4) of Theorem 2.5 on the equation
$XX^{-1}=I_{n}$ (here $X=(X_{ij})\in\mathbb{R}^{n\times n}$), we have

$(\frac{d(X)}{dX})\times_{4}X^{-1}+X\times_{2}\frac{dX^{-1}}{dX}$

It follows that

$\displaystyle\frac{dX^{-1}}{dX}$ $\displaystyle=$
$\displaystyle-X^{-1}\times_{2}[\frac{d(X)}{dX}\times_{4}X^{-1}]$
$\displaystyle=$ $\displaystyle-X^{-1}\times_{2}[(I_{n}\times
I_{n})\times_{4}X^{-1}]$ $\displaystyle=$
$\displaystyle-X^{-1}\times_{2}[(I_{n}\times X^{-1}]$ $\displaystyle=$
$\displaystyle-X^{-\top}\times X^{-1}$

∎

Given any two 4-order tensors, say,

$\mathcal{A}=(A_{i_{1}i_{2}j_{1}j_{2}})\in\mathbb{R}^{m_{1}\times m_{2}\times
n_{1}\times
n_{2}},\mathcal{B}=(B_{i_{1}i_{2}j_{1}j_{2}})\in\mathbb{R}^{n_{1}\times
n_{2}\times q_{1}\times q_{2}}.$

The product of $\mathcal{A},\mathcal{B}$, denoted by $\mathcal{A}\mathcal{B}$,
is referred to as the 4-order tensor of size $m_{1}\times m_{2}\times
q_{1}\times q_{2}$ whose entries are defined by

$(\mathcal{A}\mathcal{B})_{i_{1}i_{2}j_{1}j_{2}}=\sum\limits_{k_{1},k_{2}}A_{i_{1}i_{2}k_{1}k_{2}}B_{k_{1}j_{2}j_{1}j_{2}}$

This definition can also be carried over along other pair of directions. We
will not go into detail at this point in this paper, and want to point out
that all the results concerning the tensor forms of the derivatives can be
transformed into the conventional matrix forms, which can be achieved by
Kronecker product.

##  3 On Gaussian matrices

In this section, we introduce and study the random matrices with Gaussian
distributions and investigate the tensor products of such matrices. We denote
$\parallel\\!\cdot\\!\parallel$ for the Euclidean norm and
$S^{k-1}:=\left\\{s\in\mathbb{R}^{k}:\parallel\\!s\\!\parallel=1\right\\}$ for
the _unit sphere_ in $\mathbb{R}^{k}$ for any positive integer $k>1$. Let
$m,n>1$ be two positive integers. Then $\alpha\in S^{m-1},\beta\in S^{n-1}$
implies $\alpha\otimes\beta\in S^{mn-1}$ since
$\parallel\\!\alpha\otimes\beta\\!\parallel=\parallel\\!\alpha\\!\parallel\cdot{}\parallel\\!\beta\\!\parallel=1$.
Let $\textbf{X}=(x_{ij})\in\mathbb{R}^{m\times n}$ be a random matrix. The
characteristic function (CF) of X is defined by

$\phi_{\textbf{X}}(T)=E[exp(\imath\texttt{Tr}(T^{\prime}\textbf{X}))],\quad\forall
T\in\mathbb{R}^{m\times n}$

Note that $\phi_{\textbf{X}}(T)=\phi_{\textbf{x}}(\textbf{t})$ where
$\textbf{x}=\texttt{vec}(\textbf{X})$ and $\textbf{t}=\texttt{vec}(T)$ are
respectively the vectorization of X and $T$. While vectorization allows us to
treat all derivative (and thus the high order moments) of random matrices, it
also pose a big challenge for identifying the $k$-moment corresponding to each
coordinate of X. We introduce the tensor expression for all these basic
terminology thereafter. For our convenience, we denote by $X_{i\cdot}$
($X_{\cdot j}$) the $i$th row (resp. $j$th column) of a random matrix X. X is
called a _standard normally distributed_ (_SND_) or a _SND_ matrix if

(1)

$X_{i\cdot}$’s i.i.d. with $X_{i\cdot}\sim\mathcal{N}_{n}(0,I_{n})$;

(2)

$X_{\cdot j}$’s i.i.d. with $X_{\cdot j}\sim\mathcal{N}_{m}(0,I_{m})$.

that is, all the rows (and columns) of $X$ are i.i.d. with standard normal
distribution. This is denoted by $X\sim\mathcal{N}_{m,n}(0,I_{m},I_{n})$.
The following lemma concerning a necessary and sufficient condition for a SND
random vector will be frequently used in the paper.

######  Lemma 3.1.

Let $\textbf{x}\in\mathbb{R}^{n}$ be a random vector. Then
$\textbf{x}\sim\mathcal{N}_{n}(0,I_{n})$ if and only if
$\alpha^{\top}\textbf{x}\sim\mathcal{N}(0,1)$ for all unit vectors $\alpha\in
S^{n-1}$.

Lemma 3 is immediate from the fact (see e.g. [16]) that
$\textbf{x}\sim\mathcal{N}_{n}(\mu,\Sigma)$ if and only if

$\alpha^{\top}\textbf{x}\sim\mathcal{N}(\alpha^{\top}\mu,\alpha^{\top}\Sigma\alpha),\quad\forall\alpha\in
S^{n-1}.$

The following lemma presents equivalent conditions for a random matrix to be
SND:

######  Lemma 3.2.

Let $X=(x_{ij})\in\mathbb{R}^{m\times n}$ be a random matrix. Then the
following conditions are equivalent:

(1)

$X\sim\mathcal{N}_{m,n}(0,I_{m},I_{n})$.

(2)

$\texttt{vec}(X)\sim\mathcal{N}_{mn}(0,I_{mn})$.

(3)

All $x_{ij}$ are i.i.d. with $x_{ij}\sim\mathcal{N}(0,1)$.

(4)

$\alpha^{\top}X\beta\sim\mathcal{N}(0,1),\quad\forall\alpha\in
S^{m-1},\beta\in\mathbb{R}^{n-1}$.

###### Proof.

The equivalence of (2) and (3) is obvious. We now show
$(1)\Leftrightarrow(2)\Leftrightarrow(4)$. To show $(2)\implies(1)$, we denote
$\textbf{x}:=\texttt{vec}(X)\in\mathbb{R}^{mn}$ and suppose that
$\textbf{x}\sim\mathcal{N}_{mn}(0,I_{mn})$. Then $\texttt{cov}(X_{\cdot
i},X_{\cdot,j})=\Sigma_{ij}=0$ for all distinct $i,j\in[n]$. So the columns of
$X$ are independent. Furthermore, we have by Lemma 3 that

$\alpha^{\top}\textbf{x}\sim\mathcal{N}(0,1),\quad\forall\alpha\in S^{mn-1}$
(3.1)

Now set $\alpha=e_{j}\otimes\beta\in\mathbb{R}^{mn}$ ($\forall j\in[n]$) where
$e_{j}\in\mathbb{R}^{n}$ is the $j$th coordinate vector of $\mathbb{R}^{n}$
and $\beta\in S^{m-1}$. Then $\alpha\in S^{mn-1}$ and by (3.1) we have

$\beta^{\top}X_{\cdot
j}=\beta^{\top}Xe_{j}=(e_{j}^{\top}\otimes\beta^{\top})\textbf{x}=(e_{j}\otimes\beta)^{\top}\textbf{x}=\alpha^{\top}\textbf{x}\sim\mathcal{N}(0,1)$

It follows by Lemma 3 that $X_{\cdot j}\sim\mathcal{N}_{m}(0,I_{m})$ for all
$j\in[n]$. Consequently (2) implies (1) by definition.
$(2)\implies(4)$: Denote $\gamma:=\beta\otimes\alpha$ for any given $\alpha\in
S^{m-1},\beta\in S^{n-1}$. Then $\gamma\in S^{mn-1}$. Since
$\texttt{vec}(X)\sim\mathcal{N}_{mn}(0,I_{mn})$, we have, by Lemma , that
$\alpha^{\top}X\beta=(\beta^{\top}\otimes\alpha^{\top})\texttt{vec}(X)=\gamma^{\top}\texttt{vec}(X)\sim\mathcal{N}(0,1)$,
which proves (4).
To show $(4)\implies(2)$, we let $\gamma\in S^{mn-1}$. Then there is a unique
matrix $A=(a_{ij})\in\mathbb{R}^{m\times n}$ such that
$\gamma=\texttt{vec}(A)$ and
$\parallel\\!A\\!\parallel_{F}^{2}=\sum_{i,j}a_{ij}^{2}=\parallel\\!\alpha\\!\parallel^{2}=1$
($\parallel\\!A\\!\parallel_{F}$ denotes the Frobenius norm of matrix $A$).
Since
$\parallel\\!\texttt{vec}(A)\\!\parallel_{2}=\parallel\\!A\\!\parallel_{F}=1$,
$\texttt{vec}(A)\in S^{mn-1}$. By (4), we have
$\gamma^{\top}\texttt{vec}(X)=\texttt{vec}(A)^{\top}\texttt{vec}(X)\sim\mathcal{N}(0,1)$.
Consequently we get $X\sim\mathcal{N}_{mn}(0,I_{mn})$ by Lemma 3. ∎

The density function and the characteristic function of a SND random matrix
[2, 16, 14, 17] can also be obtained by Lemma 3.2.

######  Proposition 3.3.

Let $X\sim\mathcal{N}_{m,n}(0,I_{m},I_{n})$. Then

(1)

$f_{X}(T)=(2\pi)^{-mn/2}\exp\left\\{-\frac{1}{2}\texttt{Tr}(T^{\top}T)\right\\}$
where $T\in\mathbb{R}^{m\times n}$.

(2)

$\phi_{X}(T)=\exp\left\\{-\frac{1}{2}\texttt{Tr}(T^{\top}T)\right\\}$ where
$T\in\mathbb{R}^{m\times n}$.

Let $M=(m_{ij})\in\mathbb{R}^{n_{1}\times n_{2}}$, and
$\Sigma_{k}=(\sigma_{ij}^{(k)})\in\mathbb{R}^{n_{k}\times n_{k}}$ be positive
definite for $k=1,2$. A random matrix $X=(X_{ij})\in\mathbb{R}^{n_{1}\times
n_{2}}$ is called a _Gaussian matrix_ 111We do not use the term normal matrix
since it is referred to a matrix satisfying $XX^{\top}=X^{\top}X$ with
parameters $(M,\Sigma_{1},\Sigma_{2})$, written as
$X\sim\mathcal{N}_{m,n}(M,\Sigma_{1},\Sigma_{2})$, if

(a)

Each row $X_{i\cdot}$ follows a Gaussian distribution with

$X_{i\cdot}\sim\mathcal{N}_{n}(M_{i\cdot},\sigma_{ii}^{(1)}\Sigma_{2}),\quad\forall
i\in[m],$ (3.2)

(b)

Each column vector $X_{\cdot j}$ follows a Gaussian distribution with

$X_{\cdot j}\sim\mathcal{N}_{m}(M_{\cdot
j},\sigma_{jj}^{(2)}\Sigma_{1}),\quad\forall j\in[n]$ (3.3)

Such a random matrix $X$ is called a _Gaussian matrix_. It follows that the
vectorization of a Gaussian matrix $X$ is a Gaussian vector, i.e.,

$\texttt{vec}(X)\sim\mathcal{N}_{mn}(\texttt{vec}(M),\Sigma_{2}\otimes\Sigma_{1})$
(3.4)

A Gaussian vector cannot be shaped into a Gaussian matrix if its covariance
matrix possesses no Kronecker decomposition of two PSD matrices. For any two
random matrices (vectors) $X,Y$ of the same size, we denote $X=Y$ if their
distributions are identical. The following statement, which can be found in
[16], shows that an affine transformation preserves the Gaussian distribution.

######  Lemma 3.4.

Let $X\sim\mathcal{N}_{n_{1},n_{2}}(\mu,\Sigma_{1},\Sigma_{2})$ and
$Y=B_{1}XB_{2}^{\top}+C$ with $B_{i}\in\mathbb{R}^{m_{i}\times n_{i}}$
($i=1,2$) being constant matrices. Then

$Y\sim\mathcal{N}_{m_{1},m_{2}}(C+B_{1}\mu
B_{2}^{\top},B_{1}\Sigma_{1}B_{1}^{\top},B_{2}\Sigma_{2}B_{2}^{\top})$

The following statement can be regarded as an alternative definition of a
Gaussian matrix.

######  Lemma 3.5.

Let $X=(x_{ij})\in\mathbb{R}^{n_{1}\times n_{2}}$ be a random matrix and
$\Sigma_{i}=A_{i}A_{i}^{\top}$ with each $A_{i}\in\mathbb{R}^{n_{i}\times
n_{i}}$ nonsingluar ($i=1,2$). Then
$X\sim\mathcal{N}_{n_{1},n_{2}}(M,\Sigma_{1},\Sigma_{2})$ if and only if there
exist a SND random matrix $Z\in\mathbb{R}^{p\times q}$ such that

$X=A_{1}ZA_{2}^{\top}+M$ (3.5)

where $M\in\mathbb{R}^{m\times n}$ is a constant matrix.

The proof of Lemma 3.5 can be found in [2]. Now we let
$X\sim\mathcal{N}_{n_{1},n_{2}}(M,\Sigma_{1},\Sigma_{2})$ be a Gaussian matrix
where $M=(m_{ij})\in\mathbb{R}^{n_{1}\times n_{2}}$ and
$\Sigma_{k}=(\sigma_{ij}^{(k)})\in\mathbb{R}^{n_{k}\times n_{k}}$ ($k=1,2$)
being positive definite. Write

$\omega_{T}=\texttt{Tr}[(T-M)^{\top}\Sigma_{1}^{-1}(T-M)\Sigma_{2}^{-1}]$
(3.6)

where $T\in\mathbb{R}^{n_{1}\times n_{2}}$ is arbitrary. The characteristic
function (CF) of $X$ is defined as

$\phi_{X}(T):=E[\exp(\imath\left<X,T\right>)],\quad
T\in\mathbb{R}^{n_{1}\times n_{2}}$ (3.7)

We have

######  Corollary 3.6.

Let $X\in\mathbb{R}^{n_{1}\times n_{2}}$ be a random matrix. Then the density
and the characteristic function of $X$ are respectively

$f_{X}(T)=(2\pi)^{-n_{1}n_{2}/2}(\det(\Sigma_{1}))^{-n_{2}/2}(\det(\Sigma_{2}))^{-n_{1}/2}\exp\left\\{-\frac{1}{2}\omega_{T}\right\\}$
(3.8)

and

$\phi_{X}(T)=\exp\left\\{\imath\texttt{Tr}(T^{\top}M)-\frac{1}{2}\texttt{Tr}(T^{\top}\Sigma_{1}T\Sigma_{2})\right\\}$
(3.9)

where $T$ takes values in $\mathbb{R}^{n_{1}\times n_{2}}$.

Lemma 3.4 can be used to justify the definition of Gaussian matrices if we
take Lemma 3.5 as the original one. Let
$X\sim\mathcal{N}_{n_{1},n_{2}}(\mu,\Sigma_{1},\Sigma_{2})$, $A=I_{m}$ and
$B=e_{j}$ ($\forall j\in[n]$) is the $j$th coordinate vector of
$\mathbb{R}^{n}$. Then we have

$AXB=X_{\cdot j},AMB=\mu_{\cdot
j},A\Sigma_{1}A^{\top}=\Sigma_{1},B\Sigma_{2}B^{\top}=\sigma_{jj}^{2},$

Thus $\textbf{x}_{j}\sim\mathcal{N}_{m,1}(\mu_{\cdot
j},\Sigma_{1},\sigma_{jj}^{2})$ by Lemma 3.4, which is equivalent to (3.3).
This argument also applies to prove (3.2). Furthermore,
$x_{ij}\sim\mathcal{N}(\mu_{ij},(\sigma^{(1)}_{ii}\sigma^{(2)}_{jj})^{2})$ for
all $i\in[n_{1}],j\in[n_{2}]$.
For any matrix $A\in\mathbb{R}^{m\times n}$, we use $A[S_{1}|S_{2}]$ to denote
the submatrix of $A$ whose entries $a_{ij}$’s are confined in $i\in S_{1},j\in
S_{2}$ where $\emptyset\neq S_{1}\subset[m],\emptyset\neq S_{2}\subset[n]$.
This is denoted by $A[S]$ when $S_{1}=S_{2}=S$. It follows from Lemma 3.4 that
any submatrix of a Gaussian matrix is also Gaussian.

######  Corollary 3.7.

Let $X\sim\mathcal{N}_{n_{1},n_{2}}(\mu,\Sigma_{1},\Sigma_{2})$. and $\neq
S_{i}\subset[n_{i}]$ with cardinality $\left|S_{i}\right|=r_{i}$ for $i=1,2$.
Then

$X[S_{1}|S_{2}]\sim\mathcal{N}_{r_{1},r_{2}}(\mu[S_{1}|S_{2}],\Sigma_{1}[S_{1}],\Sigma_{2}[S_{2}])$
(3.10)

###### Proof.

We may assume that
$S_{1}=\left\\{i_{1}<i_{2}<\ldots<i_{r_{1}}\right\\},S_{2}=\left\\{j_{1}<j_{2}<\ldots<j_{r_{2}}\right\\}$,
and for $i=1,2$, choose matrix

$P_{1}^{\top}=[e_{i_{1}},e_{i_{2}},\ldots,e_{i_{r_{1}}}],\quad
P_{2}^{\top}=[f_{j_{1}},f_{j_{2}},\ldots,f_{j_{r_{2}}}]$

where $e_{k}\in\mathbb{R}^{n_{1}}$ is the $k$th coordinate (column) vector in
$\mathbb{R}^{n_{1}}$, and $f_{k}\in\mathbb{R}^{n_{2}}$ is the $k$th coordinate
(column) vector in $\mathbb{R}^{n_{2}}$, thus we have
$P_{i}\in\mathbb{R}^{r_{i}\times n_{i}}$ for $i=1,2$. Since
$X[S_{1}|S_{2}]=P_{1}XP_{2}^{\top}$, we have

$X[S_{1}|S_{2}]\sim\mathcal{N}_{r_{1},r_{2}}(P_{1}\mu
P_{2}^{\top},P_{1}\Sigma_{1}P_{1}^{\top},P_{2}\Sigma_{2}P_{2}^{\top})$

Then (3.10) follows by noticing that $\mu[S_{1}|S_{2}]=P_{1}\mu P_{2}^{\top}$
and $\Sigma_{i}[S_{i}]=P_{i}\Sigma_{i}P_{i}^{\top}$ ($i=1,2$). ∎

For a random vector $\textbf{x}\sim\mathcal{N}_{n}(\mu,\Sigma)$, we have
$\textbf{x}=\mu+A\textbf{y}$ with
$\mu\in\mathbb{R}^{n},A\in\mathbb{R}^{n\times n}$ satisfying
$AA^{\top}=\Sigma$ and $\textbf{y}\sim\mathcal{N}_{n}(0,I_{n})$. It follows
that $m_{1}[\textbf{x}]=\mu$ and

$\displaystyle m_{2}[\textbf{x}]$ $\displaystyle=$ $\displaystyle
E[\textbf{x}\textbf{x}^{\top}]=E[(\mu+A\textbf{y})(\mu+A\textbf{y})^{\top}]$
$\displaystyle=$
$\displaystyle\mu\mu^{\top}+AE[\textbf{y}\textbf{y}^{\prime}]A^{\top}$
$\displaystyle=$ $\displaystyle\mu\mu^{\top}+AA^{\top}=\mu\mu^{\top}+\Sigma$

since $E[\textbf{y}\textbf{y}^{\top}]=m_{2}[\textbf{y}]=I_{n}$. The $k$-moment
of a random matrix $X\in\mathbb{R}^{m\times n}$ is defined as the $2k$-order
tensor $E[X^{(k)}]$ which is of size $m^{[k]}\times n^{[k]}$. Write
$m_{k}[X]=(\mu_{i_{1}i_{2}\ldots i_{k}j_{1}j_{2}\ldots j_{k}})$. By definition

$\mu_{i_{1}i_{2}\ldots i_{k}j_{1}j_{2}\ldots
j_{k}}=E[x_{i_{1}j_{1}}x_{i_{2}j_{2}}\ldots x_{i_{k}j_{k}}]$

######  Lemma 3.8.

Let $Y\in\mathbb{R}^{m\times n}$ be a standard normal matrix (SNM)
$Y\sim\mathcal{N}_{m,n}(0,I_{m},I_{n})$. Then $m_{1}[Y]=0$ and

$m_{2}[Y]=E[Y\times Y]=I_{m}\times_{c}I_{n}$ (3.11)

###### Proof.

The hypothesis in the Lemma is equivalent to
$\textbf{y}:=\texttt{vec}(Y)\sim\mathcal{N}_{mn}(0,I_{mn})$ which implies
$m_{1}[\textbf{y}]=0,m_{2}[\textbf{y}]=I_{mn}$ by Lemma 3.2. Now we let
$Z:=(z_{i_{1}i_{2}i_{3}i_{4}})=Y\times Y$ be defined by
$z_{i_{1}i_{2}i_{3}i_{4}}=y_{i_{1}i_{2}}y_{i_{3}i_{4}}$. Thus we have

$M^{(2)}_{i_{1}i_{2}i_{3}i_{4}}=E[y_{i_{1}i_{2}}y_{i_{3}i_{4}}]=\texttt{cov}(y_{i_{1}i_{2}},y_{i_{3}i_{4}})=\delta_{i_{1}i_{3}}\delta_{i_{2}i_{4}}$
(3.12)

which is immediately followed by (3.11). The holding of the last equality is
due to Lemma 3.2. ∎

Now we let $X$ be a random matrix following matrix normal distribution
$X\sim\mathcal{N}_{n_{1},n_{2}}(\mu,\Sigma_{1},\Sigma_{2})$ for some constant
matrix $\mu\in\mathbb{R}^{n_{1}\times n_{2}}$ and positive semidefinite
matrices $\Sigma_{i}\in\mathbb{R}^{n_{i}\times n_{i}}$ satisfying
$\Sigma_{i}=A_{i}A_{i}^{\top},A_{i}\in\mathbb{R}^{n_{i}\times n_{i}}$
($i=1,2$) where each $A_{i}$ is nonsingular. Write
$Z=A_{1}^{-1}(X-\mu)A_{2}^{-\top}$. Then $Z\in\mathbb{R}^{n_{1}\times n_{2}}$
is a SND matrix and we have

$X=\mu+A_{1}ZA_{2}^{\top}$ (3.13)

It follows that $m_{1}[X]=E[X]=\mu\in\mathbb{R}^{n_{1}\times n_{2}}$ and
$m_{2}[X]=E[X\times X]$ whose entries are defined by

$m^{(2)}_{i_{1}i_{2}i_{3}i_{4}}=E[x_{i_{1}i_{2}}x_{i_{3}i_{4}}]=\texttt{cov}(X_{i_{1}i_{2}},X_{i_{3}i_{4}})$
(3.14)

Thus $m_{2}[X]\in\mathbb{R}^{n_{1}\times n_{2}\times n_{1}\times n_{2}}$. An
_$k$ -moment_ of $X$ is defined as an $2k$-order tensor $M_{k}[X]:=E[X^{[k]}]$
(with size $(n_{1}\times n_{2})^{[k]}$). Then each entry of $M_{k}[X]$ can be
described as

$(M_{k}[X])_{i_{1}i_{2}\ldots i_{k}j_{1}j_{2}\ldots
j_{k}}=E[X_{i_{1}j_{1}}X_{i_{2}j_{2}}\ldots X_{i_{k}j_{k}}]$

For any matrices $A,B,C,D$ and non-overlapped subset
$\left\\{s_{i},t_{i}\right\\}\subset[8]$ with $s_{i}<t_{i}$. The tensor

$\mathcal{T}=(T_{i_{1}i_{2}\ldots i_{8}})\equiv
A\times_{(s_{2},t_{2})}B\times_{(s_{3},t_{3})}C\times_{(s_{4},t_{4})}D$

yields an 8-order tensor whose entries are defined by

$T_{i_{1}i_{2}\ldots
i_{8}}=A_{i_{s_{1}}i_{t_{1}}}B_{i_{s_{2}}i_{t_{2}}}C_{i_{s_{3}}i_{t_{3}}}D_{i_{s_{4}}i_{t_{4}}}$

where
$\left\\{s_{1},t_{1}\right\\}=(\cup_{k=2}^{4}\left\\{s_{k},t_{k}\right\\})^{c}$.

Given any tensors $\mathcal{A}=(A_{i_{1}i_{2}\ldots
i_{p}}),\mathcal{B}=(B_{j_{1}j_{2}\ldots j_{q}})$, resp. of order $p$ and
order $q$, we define the _tensor product_
$\mathcal{C}:=\mathcal{A}\times\mathcal{B}$ of $\mathcal{A}$ and $\mathcal{B}$
as the $(p+q)$-order tensor whose elements are defined by

$C_{i_{1}i_{2}\ldots i_{p}j_{1}j_{2}\ldots j_{q}}=A_{i_{1}i_{2}\ldots
i_{p}}B_{j_{1}j_{2}\ldots j_{q}}$

For $p=q$, we can also define the _cross tensor product_ of
$\mathcal{A},\mathcal{B}$, as the $2p$-order tensor
$\mathcal{D}=\mathcal{A}\times_{c}\mathcal{B}$ defined by

$D_{i_{1}j_{1}i_{2}j_{2}\ldots i_{p}j_{p}}=A_{i_{1}i_{2}\ldots
i_{p}}B_{j_{1}j_{2}\ldots j_{p}}$

Note that when $A\in\mathbb{R}^{m\times n},B\in\mathbb{R}^{m\times n}$, we
have

$A\times B=A\times_{(3,4)}B\in\mathbb{R}^{m\times n\times m\times n},\quad
A\times_{c}B=A\times_{(2,4)}B\in\mathbb{R}^{m\times m\times n\times n}.$

In the following when we write $A\times B$ we usually mean $A\times_{(3,4)}B$,
the tensor product of $A$ and $B$, if not mentioned the other way.

######  Lemma 3.9.

Let $A(T):=\imath\mu-\Sigma_{1}T\Sigma_{2}$. Then we have

(1)

$\frac{dA}{dT}=-\Sigma_{1}\times\Sigma_{2}$.

(2)

$\frac{d(A\times_{c}A)}{dT}=\frac{dA}{dT}\times_{(3,6)}A+\frac{dA}{dT}\times_{c}A$.

###### Proof.

We write $\mathcal{B}=(B_{i_{1}i_{2}j_{1}j_{2}})=\frac{dA(T)}{dT}$ i.e., the
derivative of $A(T)$ w.r.t. $T$, which, according to the definition, is of
size $n_{1}\times n_{1}\times n_{2}\times n_{2}$. Then

$B_{i_{1}i_{2}j_{1}j_{2}}=\frac{dA_{i_{2}j_{2}}}{dT_{i_{1}j_{1}}}=-\frac{d}{dT_{i_{1}j_{1}}}\sum\limits_{k,l}\sigma_{i_{2}k}^{(1)}\sigma_{lj_{2}}^{(2)}T_{kl}=-\sum\limits_{k,l}\sigma_{i_{2}k}^{(1)}\sigma_{lj_{2}}^{(2)}\delta_{i_{1}k}\delta_{j_{1}l}=-\sigma_{i_{2}i_{1}}^{(1)}\sigma_{j_{1}j_{2}}^{(2)}$
(3.15)

It follows that $\mathcal{B}=\frac{dA}{dT}=-\Sigma_{1}\times\Sigma_{2}$.
To prove (2), we first note that $A\times_{c}A\in\mathbb{R}^{n_{1}\times
n_{1}\times n_{2}\times n_{2}}$. Thus $\frac{d(A\times_{c}A)}{dT}$ is an
6-order tensor. Denote
$\mathcal{C}=(C_{i_{1}i_{2}i_{3}j_{1}j_{2}j_{3}})=\frac{d(A\times_{c}A)}{dT}$.
Then $\mathcal{C}$ is of size $n_{1}^{[3]}\times n_{2}^{[3]}\equiv n_{1}\times
n_{1}\times n_{1}\times n_{2}\times n_{2}\times n_{2}$, and

$C_{i_{1}i_{2}i_{3}j_{1}j_{2}j_{3}}=\frac{d[A_{i_{2}j_{2}}A_{i_{3}j_{3}}]}{dT_{i_{1}j_{1}}}=(\frac{dA}{dT})_{i_{1}i_{2}j_{1}j_{2}}A_{i_{3}j_{3}}+(\frac{dA}{dT})_{i_{1}i_{3}j_{1}j_{3}}A_{i_{2}j_{2}}$
(3.16)

for all possible $i_{k},j_{k}$, which completes the proof of (2). ∎

Let $\textbf{X}\sim\mathcal{N}_{n_{1},n_{2}}(\mu,\Sigma_{1},\Sigma_{2})$ and
let $\phi:=\phi_{X}(T)$ be its characteristic function. For our convenience,
we denote $\phi^{\prime}$ ($A^{\prime}$) for the first order derivative of
$\phi$ ($A(T)$) w.r.t. $T$, $\phi^{\prime\prime}$ ($A^{\prime\prime}$) for the
second order derivative of $\phi$ ($A(T)$) w.r.t. $T$, and $\phi^{(k)}$
($A^{(k)}$) for the $k$-order derivative of $\phi$ ($A(T)$)w.r.t. $T$. By
Lemma 3.9, we can characterize the derivatives of the characteristic function
of a Gaussian matrix $X$, as in the following:

######  Theorem 3.10.

Let $\textbf{X}\sim\mathcal{N}_{n_{1},n_{2}}(\mu,\Sigma_{1},\Sigma_{2})$,
$\phi:=\phi_{X}(T)$ be its characteristic function and $A=A(T)$ be defined as
above. Then

(1)

$\phi^{\prime}=\phi(\imath\mu-\Sigma_{1}T\Sigma_{2})=\phi A$.

(2)

$\phi^{\prime\prime}=-\phi\left[\mu^{[2]}+\Sigma_{1}\times\Sigma_{2}-\Sigma_{1}^{[2]}T^{[2]}\Sigma_{2}^{[2]}+\imath(I_{n_{1}}\times\Sigma_{1})(\mu\times
T+T\times\mu)(I_{n_{2}}\times\Sigma_{2})\right]$

(3)

$\phi^{(3)}=A\times_{(3,6)}\phi^{\prime\prime}+\phi(A\times A)^{\prime}$.

(4)

$\phi^{(k+1)}=\sum\limits_{i+j=k}A^{(i)}\times_{(3,6)}\phi^{(j)}+\sum\limits_{i+j=k-1}\phi^{(i)}\times_{(2,5)}A^{(j)}$.

###### Proof.

In order to prove (1), we denote
$f=\texttt{Tr}(\Sigma_{1}T\Sigma_{2}T^{\prime})$. Then
$\frac{df}{dT}\in\mathbb{R}^{n_{1}\times n_{2}}$. Since

$\texttt{Tr}(\Sigma_{1}T\Sigma_{2}T^{\prime})=(\Sigma_{1}T,(T\Sigma_{2})^{\prime})=\sum\limits_{i,j,k,l}\sigma^{(1)}_{ij}\sigma^{(2)}_{lk}t_{jk}t_{il},$

It follows that for any pair $(u,v)$ where $u\in[n_{1}],v\in[n_{2}]$, we have

$\displaystyle(df/dT)_{uv}=df/dt_{uv}$ $\displaystyle=$
$\displaystyle\sum\limits_{i,j,k,l}\sigma^{(1)}_{ij}\sigma^{(2)}_{lk}(\delta_{ju}\delta_{kv}t_{il}+\delta_{iu}\delta_{lv}t_{jk})$
$\displaystyle=$
$\displaystyle\sum\limits_{i,l}\sigma^{(1)}_{iu}\sigma^{(2)}_{lv}t_{il}+\sum\limits_{j,k}\sigma^{(1)}_{uj}\sigma^{(2)}_{vk}t_{jk}$
$\displaystyle=$ $\displaystyle 2(\Sigma_{1}T\Sigma_{2})_{uv}$

Thus we have

$\frac{d(\texttt{Tr}(\Sigma_{1}T\Sigma_{2}T^{\prime}))}{dT}=2\Sigma_{1}T\Sigma_{2}$
(3.17)

It follows that
$\phi^{\prime}(T):=\frac{d\phi}{dT}=\phi(\imath\mu-\Sigma_{1}T\Sigma_{2})$ due
to (LABEL:eq:cf4X). Thus (1) is proved.
To prove (2), we denote $A:=A(T)=\imath\mu-\Sigma_{1}T\Sigma_{2}$ as in Lemma
3.9. Then again by Lemma 3.9 we get

$\displaystyle\frac{d^{2}\phi}{dT^{2}}$ $\displaystyle=$
$\displaystyle\frac{d}{dT}(\phi^{\prime})=\frac{d}{dT}(\phi
A)=\frac{d\phi}{dT}\times_{(2,4)}A+\phi(\frac{dA}{dT})$ $\displaystyle=$
$\displaystyle-\phi[\mu\times_{(2,4)}\mu+\Sigma_{1}\times_{(2,4)}\Sigma_{2}-U\times_{(2,4)}U+\imath(\mu\times_{(2,4)}U+U\times_{(2,4)}\mu)]$

where $U=\Sigma_{1}T\Sigma_{2}$. By (4) of Lemma 2.1, $U\times
U=(\Sigma_{1}\times\Sigma_{1})(T\times T)(\Sigma_{1}\times\Sigma_{1})$ and
thus (2) holds.

Now (3) can be verified by using Lemma 3.9, and (4) is also immediate if we
use the induction approach to take care of it. ∎

######  Corollary 3.11.

Let $\textbf{X}\sim\mathcal{N}_{n_{1},n_{2}}(0,\Sigma_{1},\Sigma_{2})$. Then

(1)

$m_{2}[\textbf{X}]=\Sigma_{1}\times\Sigma_{2}$ .

(2)

$m_{k}[\textbf{X}]=0$ for all odd $k$.

(3)

$m_{4}[\textbf{X}]=\Sigma_{1}\times_{(2,4)}\Sigma_{2}\times_{(5,7)}\Sigma_{1}\times_{(6,8)}\Sigma_{2}+\Sigma_{1}\times_{(2,4)}\Sigma_{2}\times_{(3,7)}\Sigma_{1}\times_{(6,8)}\Sigma_{2}+\Sigma_{1}\times_{(2,6)}\Sigma_{2}\times_{(5,7)}\Sigma_{1}\times_{(4,8)}\Sigma_{2}$

###### Proof.

(1). By definition we have

$m_{2}[\textbf{X}]=\frac{1}{\imath^{2}}\phi^{\prime\prime}(T)_{T=0}=-\phi^{\prime\prime}(T)_{T=0}.$

The result is followed by (2) of Theorem 3.10.

(2). It is obvious that $m_{1}[\textbf{X}]=0$. By the hypothesis, we have
$A(0)=0,A^{\prime}(0)=-\Sigma_{1}\times\Sigma_{2}$. Thus we have

$(A\times
A)^{\prime}\mid_{T=0}=A^{\prime}(0)\times_{(3,6)}A(0)+A^{\prime}(0)\times_{c}A(0)=0,$

We now use the induction to $k$ to prove (2). By Theorem 3.10, we have

$\displaystyle m_{3}[\textbf{X}]$ $\displaystyle=$
$\displaystyle\frac{1}{\imath^{3}}\phi^{(3)}(T)_{T=0}$ $\displaystyle=$
$\displaystyle-\imath\left[A(0)\times_{(3,6)}\phi^{\prime\prime}(0)+\phi(A\times
A)^{\prime}\right]|_{T=0}$ $\displaystyle=$ $\displaystyle-\imath
A(0)\times_{(3,6)}\phi^{\prime\prime}(0)=0$

Now we assume the result holds for an odd number $k$. Then by (4) of Theorem
3.10, we have

$\displaystyle m_{k+2}$ $\displaystyle=$
$\displaystyle\frac{1}{\imath^{k+2}}\phi^{(k+2)}\mid_{T=0}$ $\displaystyle=$
$\displaystyle\frac{1}{\imath^{k+2}}\left(A\times_{(3,6)}\phi^{(k+1)}+A^{\prime}\times_{(3,6)}\phi^{(k)}+A\times_{(3,6)}\phi^{(k+1)}\right)\mid_{T=0}$
$\displaystyle=$
$\displaystyle\frac{1}{\imath^{k+2}}\left(A(0)\times_{(3,6)}\phi^{(k+1)}(0)+A^{\prime}(0)\times_{(3,6)}\phi^{(k)}(0)+A(0)\times_{(3,6)}\phi^{(k+1)}(0)\right)$
$\displaystyle=$ $\displaystyle 0$

since $A(0)=0$ and $\phi^{(k)}(0)$ by the hypothesis. Thus the result is
proved.

(3). This can be shown by using (4) of Theorem 3.10. But we can also prove it
by comparing the item (iv) in Theorem 2.2.7 (Page 203) in [16]. ∎

##  4 Random tensors with Gaussian distributions

A tensor $\mathcal{A}=(A_{i_{1}i_{2}\ldots i_{m}})$ is called a _random
tensor_ provided that each of its entries is a random variable. In this
section, we first consider 3-order random tensors
$\mathcal{A}\in\mathbb{R}^{n_{1}\times n_{2}\times n_{3}}$. Unfolding
$\mathcal{A}$ along any mode $k$ generates a random matrix
$A[k]\in\mathbb{R}^{n_{k}\times n_{i}n_{j}}$ where
$\left\\{i,j,k\right\\}=[3]$. We wll investigate the Gaussian tensor which is
defined as a random tensor with a Gaussian-type distribution in each
direction.
We need some new notations. For any $j\in[n_{2}],k\in[n_{3}]$, we use
$A(:,j,k)$(the notation is borrowed from MATLAB) to denote the $(j,k)$-fibre
obtained from $\mathcal{A}$, and define $A(i,:,k)$ and $A(i,j,:)$
analogically. We state an equivalent definition for standard normally
distributed (SND) tensor.

######  Definition 4.1.

A random tensor $\mathcal{A}\in\mathbb{R}^{n_{1}\times n_{2}\times n_{3}}$ is
said to follow a _standard normal distribution_ (_SND_), denoted
$\mathcal{A}\sim\mathcal{N}(0,I_{n_{1}},I_{n_{2}},I_{n_{3}})$, if the
following three conditions hold

(1)

$A(:,j,k)\in\mathbb{R}^{n_{1}}$’s i.i.d. for all $j\in[n_{2}],k\in[n_{3}]$
with $\mathcal{N}_{n_{1}}(0,I_{n_{1}})$.

(2)

$A(i,:,k)\in\mathbb{R}^{n_{2}}$’s i.i.d. for all $i\in[n_{1}],k\in[n_{3}]$
with $\mathcal{N}_{n_{2}}(0,I_{n_{2}})$.

(3)

$A(i,j,:)\in\mathbb{R}^{n_{3}}$’s i.i.d. for all $i\in[n_{1}],j\in[n_{2}]$
with $\mathcal{N}_{n_{3}}(0,I_{n_{3}})$.

As in the random matrix case, we have

######  Theorem 4.2.

Let $\mathbb{X}=(X_{ijk})\in\mathbb{R}^{n_{1}\times n_{2}\times n_{3}}$ be a
random tensor. Denote $n=n_{1}n_{2}n_{3}$ and $m_{l}=n/n_{l}$ for $l=1,2,3$.
The following items are equivalent:

(1)

$\mathbb{X}\sim\mathcal{N}(0,I_{n_{1}},I_{n_{2}},I_{n_{3}})$.

(2)

$X[l]\sim\mathcal{N}(0,I_{n_{l}},I_{m_{l}})$ for all $l=1,2,3$.

(3)

$\texttt{vec}(\mathbb{X})\sim\mathcal{N}_{n}(0,I_{n})$.

(4)

All $X_{ijk}$’s are i.i.d. with $X_{ijk}\sim\mathcal{N}(0,1)$.

(5)

$\mathbb{X}\times_{1}\alpha^{(1)}\times_{2}\alpha^{(2)}\times_{3}\alpha^{(3)}\sim\mathcal{N}(0,1),\forall\alpha^{(i)}\in
S^{n_{i}-1},i=1,2,3$.

###### Proof.

We first show that $(1)\implies(2)$. It suffices to show that
$X[1]\sim\mathcal{N}_{n_{1},m_{1}}(0,I_{n_{1}},I_{m_{1}})$ since
$X[1]\in\mathbb{R}^{n_{1}\times m_{1}}$. By definition of a SND tensor, we
know that $X(:,j,k)\sim\mathcal{N}_{n_{1}}(0,I_{n_{1}})$ for all $j,k$. On the
other hand, for any $i\in[n_{1}]$, we have
$\texttt{vec}(X(i,:,:))\sim\mathcal{N}_{m_{1}}(0,I_{m_{1}})$ by
$X(i,:,:)\sim\mathcal{N}_{n_{2},n_{3}}(0,I_{n_{2}},I_{n_{3}})$ and Lemma 3.2.
Thus
$X[1](i,:)=(\texttt{vec}(X(i,:,:)))^{\top}\sim\mathcal{N}_{m_{1}}(0,I_{m_{1}})$.
The implication $(2)\implies(1)$ is obvious.
The equivalence $(2)\Leftrightarrow(3)$ is directly from Lemma 3.2 and
$(3)\Leftrightarrow(4)$ is obvious. We now prove $(3)\implies(5)$. Let
$\alpha^{(l)}\in S^{n_{l}-1}$ for $l=1,2,3$ and let
$\beta=\alpha^{(2)}\otimes\alpha^{(3)}$. Then we have $\beta\in S^{m_{1}-1}$.
Furthermore, we have

$\mathbb{X}\times_{1}\alpha^{(1)}\times_{2}\alpha^{(2)}\times_{3}\alpha^{(3)}=(\alpha^{(1)})^{\top}X[1]\beta$

which is SND by Lemma 3.2. Thus we have
$\mathbb{X}\times_{1}\alpha^{(1)}\times_{2}\alpha^{(2)}\times_{3}\alpha^{(3)}\sim\mathcal{N}(0,1)$
for all $\alpha^{(l)}\in S^{n_{l}-1}$ ($l=1,2,3$). Thus (5) holds.
Conversely, suppose (5) hold and we want to show (2). For any given $\alpha\in
S^{n_{1}-1}$, we denote $A(\alpha)=\mathbb{X}\times_{1}\alpha$. Then
$A(\alpha)\in\mathbb{R}^{n_{2}\times n_{3}}$. Furthermore, for any
$\alpha^{(2)}\in S^{n_{2}-1},\alpha^{(3)}\in S^{n_{3}-1}$, we have

$\mathbb{X}\times_{1}\alpha^{(1)}\times_{2}\alpha^{(2)}\times_{3}\alpha^{(3)}=(\alpha^{(2)})^{\top}A(\alpha)\alpha^{(3)}\in\mathcal{N}(0,1)$

It follows by Lemma 3.2 that
$A(\alpha)\in\mathcal{N}_{n_{2},n_{3}}(0,I_{n_{2}},I_{n_{3}})$ for every
$\alpha\in S^{n_{1}-1}$. Specifically, if we take
$\alpha=e_{i}\in\mathbb{R}^{n_{1}}$ to be the $i$th coordinate vector in
$\mathbb{R}^{n_{1}}$, then we have $A(\alpha)=X(i,:,:)$ ($i\in[n_{1}]$). Hence
we have

$X(i,:,:)\in\mathcal{N}_{n_{2},n_{3}}(0,I_{n_{2}},I_{n_{3}}),\forall
i\in[n_{1}].$

This shows that all the slices of $\mathbb{X}$ along the mode-1 is a SND
matrix. We can also show that all the slices (along the other two directions)
are SND matrices by employing the same technique. This complete the proof that
all the five items are equivalent. ∎

From Definition 4.1, we can see that a hypercubic random tensor
$\mathcal{A}\in\mathcal{T}_{3;n}$ ($\mathcal{T}_{3;n}:=\mathbb{R}^{n\times
n\times n}$) is SND if $\mathcal{A}\textbf{x}^{3}\sim\mathcal{N}(0,1)$ for any
unit vector $\textbf{x}\in\mathbb{R}^{n}$. It is easy to see from Definition
4.1 that

######  Lemma 4.3.

Let $\mathcal{A}=(A_{ijk})\in\mathbb{R}^{n_{1}\times n_{2}\times n_{3}}$ be a
random tensor. Then $\mathcal{A}$ is a SND tensor if and only if all
$A_{ijk}$’s are $i.i.d.$ with $\mathcal{N}(0,1)$.

###### Proof.

For each $(i,j,k)\in[n_{1}]\times[n_{2}]\times[n_{3}]$, we take

$\alpha^{1}=e^{i,1}\in\mathbb{R}^{n_{1}},\alpha^{2}=e^{j,2}\in\mathbb{R}^{n_{2}},\alpha^{3}=e^{k,3}\in\mathbb{R}^{n_{3}}$

where $e^{s,t}$ denotes the $s$th coordinate vector in $\mathbb{R}^{n_{t}}$
for any $s\in[n_{t}],t\in[3]$. Then
$\textbf{x}=\mathcal{A}\times_{1}\alpha^{1}\times_{2}\alpha^{2}\times_{3}\alpha^{3}=A_{ijk}$.
Thus $A_{ijk}\sim\mathcal{N}(0,1)$ since $\mathcal{A}$ is SND by Theorem 4.2.
The independency of $A_{ijk}$’s can be deduced from Lemma 3.2. ∎

Denote $n:=n_{1}n_{2}\ldots n_{m}$ and $m_{k}:=n/n_{k}$ where each $n_{k}>1$
is a positive integer. Let $\mu\in\mathbb{R}^{n_{1}\times\ldots\times n_{m}}$
be a constant tensor and $\Sigma_{k}\in\mathbb{R}^{n_{k}\times n_{k}}$ be
positive semidefinite for each $k\in[m]$. Suppose that
$\mathcal{A}\in\mathbb{R}^{n_{1}\times\ldots\times n_{m}}$. We denote
$A^{(k)}(:,j)$ for the $j$th column (fiber) vector of $A^{(k)}$ where $A[k]$
is the flattened matrix of $\mathcal{A}$ along mode-$k$. $\mathcal{A}$ is said
to be a _Gaussian tensor_ with parameters
$(\mathcal{M},\Sigma_{1},\ldots,\Sigma_{m})$ if
$A^{(k)}(:,1),A^{(k)}(:,2),\ldots,A^{(k)}(:,m_{k})$ are independent with

$A^{(k)}(:,j)\sim\mathcal{N}_{n_{k},m_{k}}(M^{(k)}(:,j),\lambda_{k}\Sigma_{k}),\quad\forall
k\in[m].$ (4.1)

Definition 4.1 reduces to a general multivariate Gaussian distribution when
$m=1$ and to a Gaussian matrix when $m=2$. We also note that a tensor
$\mathcal{A}\in\mathcal{T}_{m,n}$ is a Gaussian tensor if
$\mathcal{A}\textbf{x}^{m}$ follows a Gaussian distribution for every nonzero
vector $\textbf{x}\in\mathbb{R}^{n}$.

The following theorem tells that each flattened matrix of an 3-order SND
tensor along any direction is a SND matrix.

######  Theorem 4.4.

Let $\mathcal{A}\in\mathbb{R}^{n_{1}\times\ldots\times n_{m}}$ be a Gaussian
tensor with $\textrm{I}:=n_{1}\times\ldots\times n_{m}$. Then
$A[k]\in\mathbb{R}^{n_{k}\times N_{k}}$ is a Gaussian matrix for each
$k\in[m]$, where $A[k]$ is the flattened matrix of $\mathcal{A}$ along the
$k$-mode.

Theorem 4.4 is directly from Definition 4.1.
Now we consider an $m$-order tensor $\mathcal{A}\in\mathcal{T}({\rm I})$ of
size $\rm{I}:=d_{1}\times d_{2}\times\ldots\times d_{m}$ and denote $a[k,j]$
the $j$th fibre of $\mathcal{A}$ along the $k$-mode where $k\in[m]$ and $j$
ranges from 1 to $N_{k}:=d_{1}d_{2}\ldots d_{m}/d_{k}$. We call $\mathcal{A}$
a _standard Gaussian tensor_ if $a[k,j]\sim\mathcal{N}_{d_{k}}(0,I_{d_{k}})$
for each $k,j$, and denote
$\mathcal{A}\sim\mathcal{N}_{\rm{I}}(0,I_{d_{1}},\ldots,I_{d_{m}})$. A random
tensor $\mathcal{A}\in\mathcal{T}({\rm I})$ is said to follow a _Gaussian_ (or
_normal_) distribution if
$A_{[k,j]}\sim\mathcal{N}_{I_{k}}(M_{[k,j]},\Sigma_{k})$ for each $k,j$. The
following result also applies to a general case.

######  Theorem 4.5.

Let $\mathbb{Y}=(Y_{ijk})\in\mathbb{R}^{n_{1}\times n_{2}\times n_{3}}$ be a
3-order random tensor, $\mathcal{M}=(M_{ijk})\in\mathbb{R}^{n_{1}\times
n_{2}\times n_{3}}$ be a constant tensor and $U_{k}\in\mathbb{R}^{n_{k}\times
n_{k}}$ be invertible matrices. If

$\mathbb{Y}=\mathcal{M}+\mathbb{X}\times_{1}U_{1}\times_{2}U_{2}\times_{3}U_{3}$
(4.2)

where $\mathbb{X}\sim\mathcal{N}(0,I_{n_{1}},I_{n_{2}},I_{n_{3}})$ is a
standard normal tensor. Then $\mathbb{Y}$ follows a Gaussian distribution with
parameters $(\mathcal{M},\Sigma_{1},\Sigma_{2},\Sigma_{3})$ where
$\Sigma_{k}=U_{k}^{\top}U_{k}$ for $k=1,2,3$.

###### Proof.

Since

$E[\mathbb{Y}]=E[\mathcal{M}+\mathbb{X}\times_{1}U_{1}\times_{2}U_{2}\times_{3}U_{3}]=\mathcal{M}+E[\mathbb{X}]\times_{1}U_{1}\times_{2}U_{2}\times_{3}U_{3},$

which follows that $E[\mathbb{Y}]=\mathcal{M}$ due to the hypothesis that
$\mathbb{X}\sim\mathcal{N}(0,I_{n_{1}},I_{n_{2}},I_{n_{3}})$. Now we want to
show that

$Y[k]\sim\mathcal{N}_{n_{k},m_{k}}(M[k],\Sigma_{k},\Sigma_{l}\otimes\Sigma_{p})$

for each $k\in[3]$ where
$\left\\{l,p\right\\}=\left\\{1,2,3\right\\}\backslash\left\\{k\right\\}$. In
the following we need only to show this holds for $k=1$, and the other two
situations ($k=2,3$) can be handled the same way. By (4.2), we have

$Y[1]=M[1]+U_{1}^{\top}X[1](U_{2}\otimes U_{3})$ (4.3)

It follows by Lemma 3.5 that
$Y[1]\sim\mathcal{N}_{n_{1},n_{2}n_{3}}(M[1],\Sigma_{1},\Sigma_{2}\otimes\Sigma_{3})$
since $\Sigma_{1}=U_{1}^{\top}U_{1}$ and all the rows of $Y[1]$, denoted
$Y[1](i,:)$ follows $\mathcal{N}_{n_{2}n_{3}}(M[1](i,:),\Sigma_{(2,3)})$ with

$\Sigma_{(2,3)}=(U_{2}\otimes U_{3})^{\top}(U_{2}\otimes
U_{3})=(U_{2}^{\top}U_{2})\otimes(U_{3}^{\top}U_{3})=\Sigma_{2}\otimes\Sigma_{3}.$

Similarly, we can show that

$Y[2]\sim\mathcal{N}_{n_{2},n_{1}n_{3}}(M[2],\Sigma_{2},\Sigma_{1}\otimes\Sigma_{3}),\quad
Y[3]\sim\mathcal{N}_{n_{3},n_{1}n_{2}}(M[3],\Sigma_{3},\Sigma_{1}\otimes\Sigma_{2})$

By the definition, we get the result. ∎

The above result can also be easily generalized to an $m$-order Gaussian
tensor as in the following theorem.

######  Theorem 4.6.

A random tensor $\mathbb{Y}\in\mathcal{T}_{m;n}$ follows a normal distribution
$\mathbb{Y}\sim\mathcal{N}_{m;n}(\mathcal{M},\Sigma_{1},\ldots,\Sigma_{m})$
iff there exist some matrices $U_{k}$ ($k\in[m]$) such that
$\mathbb{Y}=\mathbb{X}\times_{1}U_{1}\times_{2}U_{2}\ldots\times_{m}U_{m}$ and
$\mathbb{X}$ obeys a standard Gaussian distribution.

######  Theorem 4.7.

A random tensor $\mathbb{Y}\in\mathcal{T}_{\rm{I}}$ follows a normal
distribution
$\mathbb{Y}\sim\mathcal{N}_{\rm{I}}(\mathcal{M},\Sigma_{1},\ldots,\Sigma_{m})$
iff

$\mathbb{Y}[k]\sim\mathcal{N}_{n_{k},m_{k}}(\mathcal{M}_{k},\Sigma_{k},\Omega_{k})$

where $\mathbb{Y}[k]$ is the unfolding of $\mathbb{Y}$ along mode-$k$ and
$\Omega_{k}:=\Sigma_{m}\otimes\ldots\otimes\Sigma_{k+1}\otimes\Sigma_{k-1}\ldots\otimes\Sigma_{1}$
for each $k\in[m]$.

###### Proof.

This is true for $m=1,2$ by the result on random vector and random matrix
cases. Using the unfolding of tensor $\mathbb{Y}$ and induction on $m$, we can
easily get the result. ∎

Let $\mathbb{X}\sim\mathcal{N}_{\rm{I}}(0,I_{d_{1}},\ldots,I_{d_{m}})$ be a
random following a standard normal distribution(SND). The density function of
$\mathbb{X}$ is defined by

$f_{\mathbb{X}}(\mathcal{T})=(2\pi)^{-\frac{1}{2}d}\exp(-\frac{1}{2}\left<\mathcal{T},\mathcal{T}\right>)$

where $d=d_{1}d_{2}\ldots d_{m}$.

######  Lemma 4.8.

Let $\mathbb{X}\sim\mathcal{N}_{\rm{I}}(0,I_{d_{1}},\ldots,I_{d_{m}})$ be a
random following a SND. Then the CF of $\mathbb{X}$ is

$\phi_{\mathbb{X}}(\mathcal{T})=\exp\left\\{-\frac{1}{2}\left<\mathcal{T},\mathcal{T}\right>\right\\}$
(4.4)

where $\mathcal{T}=(T_{i_{1}\ldots i_{m}})\in\mathbb{R}^{d_{1}\times
d_{2}\times\ldots\times d_{m}}$.

###### Proof.

By definition of CF, we have

$\displaystyle\phi_{\mathbb{X}}(\mathcal{T})$ $\displaystyle=$ $\displaystyle
E[\exp(\imath\left<\mathcal{T},\mathbb{X}\right>)]$ $\displaystyle=$
$\displaystyle E[\exp(\imath\sum\limits_{i_{1},\ldots,i_{m}}T_{i_{1}\ldots
i_{m}}X_{i_{1}\ldots i_{m}})]$ $\displaystyle=$ $\displaystyle
E[\prod\limits_{i_{1},\ldots,i_{m}}\exp(\imath(T_{i_{1}\ldots
i_{m}}X_{i_{1}\ldots i_{m}}))]$ $\displaystyle=$
$\displaystyle\prod\limits_{i_{1},\ldots,i_{m}}E[\exp(\imath(T_{i_{1}\ldots
i_{m}}X_{i_{1}\ldots i_{m}}))]$ $\displaystyle=$
$\displaystyle\prod\limits_{i_{1},\ldots,i_{m}}\exp(-\frac{1}{2}T_{i_{1}\ldots
i_{m}}^{2})$ $\displaystyle=$
$\displaystyle\exp[-\frac{1}{2}\sum\limits_{i_{1},\ldots,i_{m}}T_{i_{1}\ldots
i_{m}}^{2}]$ $\displaystyle=$
$\displaystyle\exp[-\frac{1}{2}\left<\mathcal{T},\mathcal{T}\right>]$

∎

An $m$-order random tensor $\mathbb{X}=(X_{i_{1}i_{2}\ldots i_{m}})$ of size
$\mathcal{I}:=n_{1}\times n_{2}\times\ldots\times n_{m}$ is called a _Gaussian
tensor_ , with parameters $(\mu,\Sigma_{1},\ldots,\Sigma_{m})$ where
$\mu\in\mathcal{T}_{\mathcal{I}}$ is of same size as that of $\mathbb{X}$, and
each $\Sigma_{k}\in\mathbb{R}^{n_{k}\times n_{k}}$ is a positive semidefinite
matrix, denoted
$\mathbb{X}\sim\mathcal{N}_{\rm{I}}(\mu,\Sigma_{1},\ldots,\Sigma_{m})$, if
there exists a SND tensor $\mathbb{Z}$ such that
$\mathbb{Z}\in\mathbb{R}^{n_{1}\times\ldots\times n_{m}}$ and

$\mathbb{X}=\mu+\mathbb{Z}\times_{1}U_{1}\times_{2}U_{2}\times_{3}\ldots\times_{m}U_{m}$
(4.5)

where $U_{k}\in\mathbb{R}^{n_{k}\times n_{k}}$ satisfies
$U_{k}U_{k}^{\top}=\Sigma_{k}$ for each $k\in[m]$. This is obviously a
generalization of a third order Gaussian tensor, and thus some results on the
third order Gaussian tensor can also apply.

Given positive semidefinite matrices $\Sigma_{1},\ldots,\Sigma_{m}$ with
$\Sigma_{k}\in\mathbb{R}^{n_{k}\times n_{k}}$. We denote $\Omega_{k}$ as the
Kronecker product
$\Sigma_{1}\otimes\Sigma_{k-1}\otimes\Sigma_{k+1}\otimes\ldots\otimes\Sigma_{m}$
where the item $\Omega_{k}$ is removed. Note that matrix
$\Omega_{k}\in\mathbb{R}^{m_{k}\times m_{k}}$ ($m_{k}=n/n_{k}$ with
$n:=n_{1}n_{2}\ldots n_{m}$ as we defined before) is also psd. We have

######  Theorem 4.9.

An $m$-order random tensor $\mathbb{X}=(X_{i_{1}i_{2}\ldots
i_{m}})\in\mathcal{T}_{\mathcal{I}}$ with $\mathcal{I}:=n_{1}\times
n_{2}\times\ldots\times n_{m}$ is a Gaussian tensor with
$\mathbb{X}\sim\mathcal{N}_{\mathcal{I}}(\mu,\Sigma_{1},\ldots,\Sigma_{m})$ if
each of its unfolded matrix
$X[k]\sim\mathcal{N}_{n_{k},m_{k}}(\mu[k],\Sigma_{k},\Omega_{k})$ for
$k\in[m]$.

Theorem 4.9 can be proved by the same technique as in the case of order three.
We need some preparations before presenting the characteristic function of a
Gaussian tensor. Recall that for any matrices $A_{k}\in\mathbb{R}^{p_{k}\times
n_{k}}$ ($k\in[m]$), the tensor product $\mathcal{A}:=A_{1}\times
A_{2}\times\ldots\times A_{m}$ is an $2m$-order tensor of size $p_{1}\times
n_{1}\times p_{2}\times n_{2}\times\ldots\times p_{m}\times n_{m}$. for any
given $m$-order tensor $\mathcal{A}=(A_{i_{1}i_{2}\ldots i_{m}})$ of size
$\mathcal{I}:=n_{1}\times\ldots\times n_{m}$, the tensor
$\mathcal{A}^{[2]}=\mathcal{A}\times\mathcal{A}$, which is called the square
of $\mathcal{A}$ (in the sense of tensor product), is an $2m$-order tensor
with size $\mathcal{I}^{[2]}$ whose entries indexed by $B_{i_{1}i_{2}\ldots
i_{m};j_{1}j_{2}\ldots j_{m}}$, are defined by

$B_{i_{1}i_{2}\ldots i_{m};j_{1}j_{2}\ldots j_{m}}=A_{i_{1}i_{2}\ldots
i_{m}}A_{j_{1}j_{2}\ldots j_{m}}$

For an $m$-order tensor $\mathcal{A}=(A_{i_{1}i_{2}\ldots i_{m}})$ of size
$\mathcal{I}:=n_{1}\times\ldots\times n_{m}$ and any matrices
$U_{k}\in\mathbb{R}^{n_{k}\times p_{k}}$ ($k\in[m]$), the product
$\hat{\mathcal{A}}:=\mathcal{A}\times_{1}U_{1}\times_{2}U_{2}\times_{3}\ldots\times_{m}U_{m}$
is defined as the $m$-order tensor whose entries are defined by

$\hat{A}_{i_{1}i_{2}\ldots
i_{m}}=\sum\limits_{j_{1},j_{2},\ldots,j_{m}}A_{j_{1}j_{2}\ldots
j_{m}}u_{j_{1}i_{1}}^{(1)}u_{j_{2}i_{2}}^{(2)}\ldots u_{j_{m}i_{m}}^{(m)}$

This is also denoted by $\mathcal{A}\cdot[U_{1},U_{2},\ldots,U_{m}]$ (see e.g.
[15]). Similarly we can define $[U_{1},U_{2},\ldots,U_{m}]\cdot\mathcal{A}$.
The following two lemmas will be used to prove our main result:

######  Lemma 4.10.

Let $\mathcal{A}$ be an $m$-order tensor of size
$\mathcal{I}:=n_{1}\times\ldots\times n_{m}$ and let
$U_{k}\in\mathbb{R}^{n_{k}\times p_{k}}$ ($k\in[m]$). Then we have

$\mathcal{A}\cdot[U_{1},U_{2},\ldots,U_{m}]=\mathcal{A}\mathcal{U}$ (4.6)

where $\mathcal{U}:=U_{1}\times U_{2}\times\ldots\times U_{m}$. Similarly,

$\mathcal{U}\mathcal{A}=[U_{1},U_{2},\ldots,U_{m}]\cdot\mathcal{A}$ (4.7)

###### Proof.

Denote $\hat{\mathcal{A}}=\mathcal{A}\cdot[U_{1},U_{2},\ldots,U_{m}]$. Then by
definition

$\hat{A}_{i_{1}i_{2}\ldots
i_{m}}=\sum\limits_{j_{1},j_{2},\ldots,j_{m}}A_{j_{1}j_{2}\ldots
j_{m}}u_{j_{1}i_{1}}^{(1)}u_{j_{2}i_{2}}^{(2)}\ldots u_{j_{m}i_{m}}^{(m)}$
(4.8)

On the other hand, we have

$\displaystyle(\mathcal{A}\mathcal{U})_{i_{1}i_{2}\ldots i_{m}}$
$\displaystyle=$
$\displaystyle\sum\limits_{j_{1},j_{2},\ldots,j_{m}}A_{j_{1}j_{2}\ldots
j_{m}}U_{j_{1}j_{2}\ldots j_{m};i_{1}i_{2}\ldots i_{m}}$ $\displaystyle=$
$\displaystyle\sum\limits_{j_{1},j_{2},\ldots,j_{m}}A_{j_{1}j_{2}\ldots
j_{m}}u_{j_{1}i_{1}}^{(1)}u_{j_{2}i_{2}}^{(2)}\ldots u_{j_{m}i_{m}}^{(m)}$

Thus (4.6) holds. The equation (4.7) can be proved similarly. ∎

Recall that the inner product of two tensors $\mathcal{A},\mathcal{B}$ of same
size, say, $\mathcal{I}:=n_{1}\times\ldots\times n_{m}$, is defined as

$\left<\mathcal{A},\mathcal{B}\right>=\sum\limits_{i_{1},i_{2},\ldots,i_{m}}A_{i_{1}i_{2}\ldots
i_{m}}B_{i_{1}i_{2}\ldots i_{m}}$

We have

######  Lemma 4.11.

Let $\mathcal{A},\mathcal{B}$ be $m$-order tensors of size
$\mathcal{I}:=n_{1}\times\ldots\times n_{m}$ and
$\mathcal{U}=(U_{i_{1}i_{2}\ldots i_{m};j_{1}j_{2}\ldots j_{m}})$ be the
tensor of size $n_{1}\times n_{2}\times\ldots\times n_{m}\times n_{1}\times
n_{2}\times\ldots\times n_{m}$. Then we have

$\left<\mathcal{A},\mathcal{B}\mathcal{U}\right>=\left<\mathcal{U}\mathcal{A},\mathcal{B}\right>$
(4.9)

Similarly,

$\left<\mathcal{A}\mathcal{U},\mathcal{B}\right>=\left<\mathcal{A},\mathcal{U}\mathcal{B}\right>$
(4.10)

###### Proof.

We first consider the left hand side of (4.9) which, by definition of the
inner product, equals

$\displaystyle\left<\mathcal{A},\mathcal{B}\mathcal{U}\right>$
$\displaystyle=$
$\displaystyle\sum\limits_{i_{1},i_{2},\ldots,i_{m}}A_{i_{1}i_{2}\ldots
i_{m}}(\mathcal{B}\mathcal{U})_{i_{1}i_{2}\ldots i_{m}}$ $\displaystyle=$
$\displaystyle\sum\limits_{i_{1},\ldots,i_{m}}A_{i_{1}\ldots
i_{m}}\left(\sum\limits_{j_{1},\ldots,j_{m}}B_{j_{1}\ldots
j_{m}}U_{j_{1}\ldots j_{m};i_{1}\ldots i_{m}}\right)$ $\displaystyle=$
$\displaystyle\sum\limits_{j_{1},\ldots,j_{m}}B_{j_{1}\ldots
j_{m}}\left(\sum\limits_{i_{1},\ldots,i_{m}}U_{j_{1}\ldots j_{m};i_{1}\ldots
i_{m}}A_{i_{1}\ldots i_{m}}\right)$ $\displaystyle=$
$\displaystyle\left<\mathcal{U}\mathcal{A},\mathcal{B}\right>$

Thus (4.9) holds. With the same technique, (4.10) can also be proved. ∎

Now we present the characteristic function of a general Gaussian tensor.

######  Theorem 4.12.

Let $\mathbb{X}=(X_{i_{1}i_{2}\ldots i_{m}})\in\mathcal{T}_{\mathcal{I}}$ be a
Gaussian tensor with
$\mathbb{X}\sim\mathcal{N}_{\mathcal{I}}(\mu,\Sigma_{1},\ldots,\Sigma_{m})$.
Then the CF of $\mathbb{X}$ is

$\phi_{\mathbb{X}}(\mathcal{T})=\exp\left\\{\imath\left<\mathcal{T},\mu\right>-\frac{1}{2}\left<\mathcal{T}^{[2]},\Sigma\right>\right\\}$
(4.11)

where $\mathcal{T}=(T_{i_{1}\ldots i_{m}})\in\mathcal{T}_{\mathcal{I}}$ and

$\Sigma=\Sigma_{1}\times\Sigma_{2}\times\ldots\times\Sigma_{m}$ (4.12)

###### Proof.

Denote $\mathcal{U}:=U_{1}\times U_{2}\times\ldots\times U_{m}$ where
$U_{k}=(u_{ij}^{(k)})\in\mathbb{R}^{n_{k}\times n_{k}}$ is a matrix satisfying
$U_{k}^{\top}U_{k}=\Sigma_{k}$ for each $k\in[m]$. Then
$\mathcal{U}=(U_{i_{1}i_{2}\ldots i_{m};j_{1}j_{2}\ldots j_{m}})$ is an
$2m$-order tensor whose entries are defined by

$U_{i_{1}i_{2}\ldots i_{m};j_{1}j_{2}\ldots
j_{m}}=u_{i_{1}j_{1}}^{(1)}u_{i_{2}j_{2}}^{(2)}\ldots u_{i_{m}j_{m}}^{(m)}$

Now that
$\mathbb{X}\sim\mathcal{N}_{\mathcal{I}}(\mu,\Sigma_{1},\ldots,\Sigma_{m})$,
it follows by definition that there exists a SND tensor $\mathbb{Z}$ with size
$n_{1}\times\ldots\times n_{m}$ and

$\mathbb{X}=\mu+\mathbb{Z}\times_{1}U_{1}\times_{2}U_{2}\times_{3}\ldots\times_{m}U_{m}$

where $U_{k}\in\mathbb{R}^{n_{k}\times n_{k}}$ satisfies
$U_{k}U_{k}^{\top}=\Sigma_{k}$. It follows that
$\mathbb{X}=\mu+\mathbb{Z}\mathcal{U}$. By definition of CF, we have

$\displaystyle\phi_{\mathbb{X}}(\mathcal{T})$ $\displaystyle=$ $\displaystyle
E[\exp(\imath\left<\mathcal{T},\mathbb{X}\right>)]$ $\displaystyle=$
$\displaystyle
E[\exp(\imath\left<\mathcal{T},\mu+\mathbb{Z}\mathcal{U}\right>)]$
$\displaystyle=$
$\displaystyle\exp\left\\{\imath\left<\mathcal{T},\mu\right>\right\\}E[\exp(\imath\left<\mathcal{T},\mathbb{Z}\mathcal{U}\right>)]$
$\displaystyle=$
$\displaystyle\exp\left\\{\imath\left<\mathcal{T},\mu\right>\right\\}E[\exp(\imath\left<\mathcal{U}\mathcal{T},\mathbb{Z}\right>)]$
$\displaystyle=$
$\displaystyle\exp\left\\{\imath\left<\mathcal{T},\mu\right>\right\\}\phi_{\mathbb{Z}}(\mathcal{U}\mathcal{T})$
$\displaystyle=$
$\displaystyle\exp\left\\{\imath\left<\mathcal{T},\mu\right>\right\\}\exp\left\\{-1/2\left<\mathcal{U}\mathcal{T},\mathcal{U}\mathcal{T}\right>\right\\}$

The last equation is due to Lemma 4.8 since $\mathbb{Z}$ is a SND tensor. The
result (4.11) will be proved if we can prove
$\left<\mathcal{U}\mathcal{T},\mathcal{U}\mathcal{T}\right>=\left<\mathcal{T}^{[2]},\Sigma\right>$.
In fact, we have

$\displaystyle\left<\mathcal{U}\mathcal{T},\mathcal{U}\mathcal{T}\right>$
$\displaystyle=$
$\displaystyle\sum\limits_{i_{1},\ldots,i_{m}}(\mathcal{U}\mathcal{T})^{2}_{i_{1}\ldots
i_{m}}$ $\displaystyle=$
$\displaystyle\sum\limits_{i_{1},\ldots,i_{m}}\left(\sum\limits_{j_{1},\ldots,j_{m}}U_{i_{1}\ldots
i_{m};j_{1}\ldots j_{m}}T_{j_{1}\ldots
j_{m}}\right)\left(\sum\limits_{k_{1},\ldots,k_{m}}U_{i_{1}\ldots
i_{m};k_{1}\ldots k_{m}}T_{k_{1}\ldots k_{m}}\right)$ $\displaystyle=$
$\displaystyle\left<\mathcal{T}^{[2]},\Sigma\right>$

The last equality is due to the definition (4.12) and
$\Sigma_{s}=U_{s}^{\top}U_{s}$, which together implies

$\displaystyle\Sigma_{j_{1}\ldots j_{m};k_{1}\ldots k_{m}}$ $\displaystyle=$
$\displaystyle\prod\limits_{s=1}^{m}\sigma^{(s)}_{j_{s}k_{s}}$
$\displaystyle=$
$\displaystyle\prod\limits_{s=1}^{m}\left(\sum\limits_{i_{s}}u_{i_{s}j_{s}}^{(s)}u_{i_{s}k_{s}}^{(s)}\right)$
$\displaystyle=$
$\displaystyle\sum\limits_{i_{1},\ldots,i_{m}}\left(\prod\limits_{s=1}^{m}u_{i_{s}j_{s}}^{(s)}u_{i_{s}k_{s}}^{(s)}\right)$
$\displaystyle=$ $\displaystyle\sum\limits_{i_{1},\ldots,i_{m}}U_{i_{1}\ldots
i_{m};j_{1}\ldots j_{m}}U_{i_{1}\ldots i_{m};k_{1}\ldots k_{m}}$

∎

We conclude the paper by pointing out that there are several equivalent
definitions for a Gaussian tensor as in the matrix case and the regular
multivariate statistics. These will be investigated in our future research and
not be covered in the current manuscript. We will also study the probability
density function of a Gaussian tensor , present some tensor Gaussian models
and give the estimation to the parameters involved in the models.

## References

  * [1] Austin Benson, David F. Gleich, and Lek-Heng Lim, The spacey random walk: a stochastic process for higher-order data, SIAM Review, 59(2017), pp. 321-345.
  * [2] M. Bilodeau and D. Brenner, Theory of Multivariate Statistics, Springer-Verlag, New York, 1999.
  * [3] U. Chatterjee and N. Chatterjee, Vector & Tensor Analysis, Academic Publishers, NY, 2010.
  * [4] Dustin Cartwright, Bernd Sturmfels, The number of eigenvalues of a tensor, Linear Algebra Appl., 438(2013), pp. 942-952.
  * [5] J. D. Carroll and J. Chang, Analysis of individual differences in multidimensional scaling via an n-way generalisation of ’Eckart Young’ decomposition, Psychometrika, 35(1970), pp. 283-319.
  * [6] P. J. Basser and S. Pajevic, A normal distribution for tensor-valued random variables: applications to diffusion tensor MRI, IEEE Transactions on Medical Imaging, 22(2003) 785-794.
  * [7] P. Comon, G. Golub, L.-H. Lim, and B. Mourrain, Symmetric tensors and symmetric tensor rank, SIAM. J. Matrix Analysis and Applications, 30(2008) 1254-1279.
  * [8] R. D. Cook, B. Li and F. Charomonte, Envelope models for parsimonious and efficient multivariate linear regression, Statistica Sinica, 20(2010) 927-1010.
  * [9] R. D. Cook, I. S. Helland and Z. Su, Envelope and partial least squares regression, J. R. Stat. Soc. Ser. B. Stat. Methodol., 75(2013) 851-877.
  * [10] R. Couillet and M. Debbah, Random matrix methods for wireless communications, Cambridge Univ Press, 2011.
  * [11] P. J. Forrester, Log-gases and random matrices, Number 34, Princeton Univ Press, 2010.
  * [12] S. Hu, Z. Huang, C. Ling and L. Qi, On determinants and eigenvalue theory of tensors, Journal of Symbolic Computation, 50(2013) 508-531.
  * [13] X. Jiang, L.-H. Lim, Y. Yao and Y. Ye, Statistical ranking and combinatorial Hodge theory, Math. Program., Ser. B 127(2011) 203-244.
  * [14] N. L. Johnson, S. Kotz and N. Balakrishnan, Discrete Multivariate Distributions, Wiley, New York, 1997.
  * [15] T. G. Kolda, B. W. Bader, Tensor Decompositions and Applications, SIAM Review, 2009.
  * [16] T. Kollo and D. von Rosen, Advanced Multivariate Statistics with Matrices, Springer, New York, 2005.
  * [17] S. Kotz, N. Balakrishnan and N. L. Johnson, Continuous Multivariate Distributions, Wiley, New York, 2000.
  * [18] L. -H. Lim, Singular values and eigenvalues of tensors: A variational approach, Proc. of 1st IEEE Int’l Workshop on Computational Advances in Multi-Sensor Adaptive Processing (CAMSAP), 2005, pp. 129-132.
  * [19] A. M. Manceura, and P. Dutilleul, Maximum likelihood estimation for the tensor normal distribution: Algorithm, minimum sample size, and empirical bias and dispersion, Journal of Computational and Applied Math., 239(2013) 37-49.
  * [20] M. L. Mehta, A non-linear differential equation and a Fredholm determinant, J. de Phys. I France 2(1992): 1721-1729.
  * [21] M. L. Mehta, Random matrices, volume 142, Academic Press, 2004.
  * [22] F. Mezzadri and N. C. Snaith, Recent perspectives in random matrix theory and number theory, Cambridge Univ Press, 2005.
  * [23] M. Ohlson, M. R. Ahmad, D. von Rosen, The multilinear normal distribution: Introduction and some basic properties, Journal of Multivariate Analysis, 113(2013) 37-47.
  * [24] L. Qi, Eigenvalues of a real supersymmetric tensor, Journal of Symbolic Computation, 40(2005): 1302-1324.
  * [25] L. Qi, Symmetric nonnegative tensors and copositive tensors, Linear Algebra Appl., 439(2013): 228-238.
  * [26] L. Qi, Hankel tensors: Associated Hankel matrices and Vandermonde decomposition, Communications in Mathematical Sciences, 13 (2015): 113-125.
  * [27] L. Qi and Z. Luo, Tensor Analysis: Spectral theory and Special tensors, SIAM press, Philadelpha, 2017.
  * [28] T. Tao, Topics in Random Matrix Theory, Graduate studies in mathematics, American Mathematical Society, 2012.
  * [29] S. Weinberg, Gravitation and Cosmology: Principles and Applications of the General Theory of Relativity, Wiley, New York, 1972.
  * [30] C. Xu, L.He and Z. Lin, Commutation matrices and commutation tensors, Linear and Multilinear Algebra, DOI: 10.1080/03081087.2018.1556242.