# Random Double Tensors Integrals

Shih Yu Chang Shih Yu Chang is with the Department of Applied Data Science,
San Jose State University, San Jose, CA, U. S. A. (e-mail:
shihyu.chang@sjsu.edu).

###### Abstract

In this work, we try to build a theory for random double tensor integrals
(DTI). We begin with the definition of DTI and discuss how randomness
structure is built upon DTI. Then, the tail bound of the unitarily invariant
norm for the random DTI is established and this bound can help us to derive
tail bounds of the unitarily invariant norm for various types of two tensors
means, e.g., arithmetic mean, geometric mean, harmonic mean, and general mean.
By associating DTI with perturbation formula, i.e., a formula to relate the
tensor-valued function difference with respect the difference of the function
input tensors, the tail bounds of the unitarily invariant norm for the
Lipschitz estimate of tensor-valued function with random tensors as arguments
are derived for vanilla case and quasi-commutator case, respectively. We also
establish the continuity property for random DTI in the sense of convergence
in the random tensor mean, and we apply this continuity property to obtain the
tail bound of the unitarily invariant norm for the derivative of the tensor-
valued function.

Index terms— Einstein product, double tensor integrals (DTI), random DTI, tail
bound, Lipschitz estimate, convergence in the random tensor mean, derivative
of tensor-valued function

##  1 Introduction

In recent years, tensors have been applied to different applications in
science and engineering [1, 2]. However, most of these applications assume
that systems modelled by tensors are deterministic and such assumption is not
always true and practical in problems involving tensor formulations. In recent
years, more research results have pioneered some theories about random tensors
[3, 4, 5]. One important question in random tensors is about concentration
behavior of random tensors. In [6], we extend Lapalace transform method and
Lieb’s concavity theorem from matrices to tensors, and apply these tools to
generalize the classical bounds associated with the names Chernoff, Bennett,
and Bernstein from the scalar to the tensor setting. In [7], this work extends
previous work by considering the tail behavior of the top $k$-largest singular
values of a function of the tensors summation, instead of the largest/smallest
singular value of the tensors summation directly (identity function) explored
in [6]. Majorization and antisymmetric tensor product tools are main
techniques utilized to establish inequalities for unitarily norms of
multivariate tensors. Random tensors summation form discussed in [6, 7] is
linear form, i.e., each summand of random tensors with degree one. In works
[8, 9], we extend the Hanson-Wright inequality for the maximum eigenvalue of
the quadratic form of random Hermitian tensors under Einstein product. We
separate the quadratic form of random tensors into diagonal summation and
coupling (non-diagonal) summation parts. For the diagonal part, we can apply
Bernstein inequality to bound the tail probability of the maximum eigenvalue
of the summation of independent random Hermitian tensors directly. For
coupling summation part, we have to apply decoupling method first, i.e.,
decoupling inequality to bound expressions with dependent random Hermitian
tensors with independent random Hermitian tensors, before applying Bernstein
inequality again to bound the tail probability of the maximum eigenvalue of
the coupling summation of independent random Hermitian tensors. Previous works
are based on tensors with Einstein products. Since Kilmer et al. introduced
the new multiplication method between two third-order tensors around 2008 and
third-order tensors with such multiplication structure are also called as
T-product tensors [10], T-product tensors have been applied to many fields in
science and engineering, such as low-rank tensor approximation, signal
processing, image feature extraction, machine learning, computer vision, and
the multi-view clustering problem, etc. The discussion about concentration
behaviors based on T-product tensors can also be found in [11, 12].

Inspired by operator mean theory (also called Kubo–Ando theory), we try to
consider other operations besides $+$ (arithmetic mean) among tensors [13].
The matrix mean for double operators can be expressed by Eq. (5:1:2) in [13],
which has the same formation of double operator integral theory discussed in
[14]. In this work, we begin to define double tensor integrals (DTI) and
consider the tail bound for the unitarily invariant norm of random DTI, see
Theorem 1. This bound can help us to establish tail bounds for various types
of tensor means besides arithmetic mean. Since DTI can be used to express
perturbation formula, i.e., an formula to relate the tensor-valued function
difference with respect the difference of the function input tensors, we
establish Lipschitz estimate for random tensors by the tail bound format, see
Theorem 2. We also generalize Lipschitz estimate for random tensors with
another quasi-commutator tensor, $\mathcal{D}$ by providing the tail bound for
the unitarily invariant norm of
$\mathcal{D}\star_{N}f(\mathcal{A})-f(\mathcal{B})\star_{N}\mathcal{D}$ in
Theorem 3, where $\mathcal{A},\mathcal{B}$ are random Hermitian tensors. We
also establish a continuity for random DTI in the sense of convergence in the
random tensor mean. This continuity property helps us to obtain the tail
bounds for the unitarily invariant norm of the derivative of the tensor-valued
function under vanilla case and quasi-commutator case.

We define double tensor integrals (DTI) and randomness of DTI in Section 2.
The tail bound for the unitarily invariant norm of random DTI and its
applications to obtain various tail bounds for different types of double
tensors means like arithmetic mean, geometric mean, harmonic mean, and general
mean, are discussed in Section 3. In Section 4, we establish Lipschitz
estimates for random tensors by the tail bound format for vanilla case and
quasi-commutator case. We will establish continuity of random DTI based on the
convergence in tensor mean of random Hermitian tensors in Section 5. The
application of DTI theory to acquire the tail bound for the unitarily
invariant norm of the derivative of the tensor-valued function is presented by
Section 6. Finally, the conclusions are given in Section 7.

##  2 Random Double Tensor Integrals

The purpose of this section is to define random double tensor integrals (DTI).
We begin with the definition of DTI in Section 2.1. In Section 2.2, we will
present what are randomness objects at DTI discussed at this work.

Without loss of generality, one can partition the dimensions of a tensor into
two groups, say $M$ and $N$ dimensions, separately. Thus, for two
order-($M$+$N$) tensors: $\mathcal{X}\stackrel{{\scriptstyle\mbox{\tiny
def}}}{{=}}(x_{i_{1},\cdots,i_{M},j_{1},\cdots,j_{N}})\in\mathbb{C}^{I_{1}\times\cdots\times
I_{M}\times J_{1}\times\cdots\times J_{N}}$ and
$\mathcal{Y}\stackrel{{\scriptstyle\mbox{\tiny
def}}}{{=}}(y_{i_{1},\cdots,i_{M},j_{1},\cdots,j_{N}})\in\mathbb{C}^{I_{1}\times\cdots\times
I_{M}\times J_{1}\times\cdots\times J_{N}}$, according to [15], the _tensor
addition_ $\mathcal{X}+\mathcal{Y}\in\mathbb{C}^{I_{1}\times\cdots\times
I_{M}\times J_{1}\times\cdots\times J_{N}}$ is given by

$\displaystyle(\mathcal{X}+\mathcal{Y})_{i_{1},\cdots,i_{M},j_{1},\cdots,j_{N}}$
$\displaystyle\stackrel{{\scriptstyle\mbox{\tiny def}}}{{=}}$ $\displaystyle
x_{i_{1},\cdots,i_{M},j_{1},\cdots,j_{N}}$ (1)
$\displaystyle+y_{i_{1},\cdots,i_{M},j_{1},\cdots,j_{N}}.$

On the other hand, for tensors $\mathcal{X}\stackrel{{\scriptstyle\mbox{\tiny
def}}}{{=}}(x_{i_{1},\cdots,i_{M},j_{1},\cdots,j_{N}})\in\mathbb{C}^{I_{1}\times\cdots\times
I_{M}\times J_{1}\times\cdots\times J_{N}}$ and
$\mathcal{Y}\stackrel{{\scriptstyle\mbox{\tiny
def}}}{{=}}(y_{j_{1},\cdots,j_{N},k_{1},\cdots,k_{L}})\in\mathbb{C}^{J_{1}\times\cdots\times
J_{N}\times K_{1}\times\cdots\times K_{L}}$, according to [15], the _Einstein
product_ (or simply referred to as _tensor product_ in this work)
$\mathcal{X}\star_{N}\mathcal{Y}\in\mathbb{C}^{I_{1}\times\cdots\times
I_{M}\times K_{1}\times\cdots\times K_{L}}$ is given by

$\displaystyle(\mathcal{X}\star_{N}\mathcal{Y})_{i_{1},\cdots,i_{M},k_{1},\cdots,k_{L}}\stackrel{{\scriptstyle\mbox{\tiny
def}}}{{=}}$ (2)
$\displaystyle\sum\limits_{j_{1},\cdots,j_{N}}x_{i_{1},\cdots,i_{M},j_{1},\cdots,j_{N}}y_{j_{1},\cdots,j_{N},k_{1},\cdots,k_{L}}.$

One can find more preliminary facts about tensors based on Einstein product in
[6, 15]. In the remaining of this paper, we will represent the scalar value
$I_{1}\times\cdots\times I_{N}$ by $\mathbb{I}_{1}^{N}$.

###  2.1 Double Tensor Integrals

From Theorem 3.2 in [16], every Hermitian tensor
$\mathcal{H}\in\mathbb{C}^{I_{1}\times\cdots\times I_{N}\times
I_{1}\times\cdots\times I_{N}}$ has the following decomposition

$\displaystyle\mathcal{H}$ $\displaystyle=$
$\displaystyle\sum\limits_{i=1}^{\mathbb{I}_{1}^{N}}\lambda_{i}\mathcal{U}_{i}\star_{1}\mathcal{U}^{H}_{i}\mbox{
\leavevmode\nobreak\ with\leavevmode\nobreak\ \leavevmode\nobreak\
$\langle\mathcal{U}_{i},\mathcal{U}_{i}\rangle=1$ and
$\langle\mathcal{U}_{i},\mathcal{U}_{j}\rangle=0$ for $i\neq j$,}$ (3)
$\displaystyle\stackrel{{\scriptstyle\mbox{\tiny def}}}{{=}}$
$\displaystyle\sum\limits_{i=1}^{\mathbb{I}_{1}^{N}}\lambda_{i}\mathcal{P}_{\mathcal{U}_{i}}$

where $\mathcal{U}_{i}\in\mathbb{C}^{I_{1}\times\cdots\times I_{N}\times 1}$,
and the tensor $\mathcal{P}_{\mathcal{U}_{i}}$ is defined as
$\mathcal{U}_{i}\star_{1}\mathcal{U}^{H}_{i}$. The values $\lambda_{i}$ are
named as _eigevalues_. A Hermitian tensor with the decomposition shown by Eq.
(3) is named as _eigen-decomposition_. A Hermitian tensor $\mathcal{H}$ is a
positive definite (or positive semi-definite) tensor if all its eigenvalues
are positive (or nonnegative).

Let $\mathcal{A},\mathcal{B}\in\mathbb{C}^{I_{1}\times\cdots\times I_{N}\times
I_{1}\times\cdots\times I_{N}}$ be Hermitian tensors with the following eigen-
decompositions:

$\displaystyle\mathcal{A}$ $\displaystyle=$
$\displaystyle\sum\limits_{i=1}^{\mathbb{I}_{1}^{N}}\lambda_{i}\mathcal{U}_{i}\star_{1}\mathcal{U}^{H}_{i}\stackrel{{\scriptstyle\mbox{\tiny
def}}}{{=}}\sum\limits_{i=1}^{\mathbb{I}_{1}^{N}}\lambda_{i}\mathcal{P}_{\mathcal{U}_{i}},$
(4)

and

$\displaystyle\mathcal{B}$ $\displaystyle=$
$\displaystyle\sum\limits_{j=1}^{\mathbb{I}_{1}^{N}}\mu_{j}\mathcal{V}_{i}\star_{1}\mathcal{V}^{H}_{j}\stackrel{{\scriptstyle\mbox{\tiny
def}}}{{=}}\sum\limits_{j=1}^{\mathbb{I}_{1}^{N}}\mu_{j}\mathcal{P}_{\mathcal{V}_{j}}.$
(5)

We define _double tensor integrals_ (DTI) with respect to tensors
$\mathcal{A},\mathcal{B}$ and the function
$\psi:\mathbb{R}^{2}\rightarrow\mathbb{R}$, denoted as
$T_{\mathcal{A},\mathcal{B},\psi}(\mathcal{X})$, which can be expressed as

$\displaystyle T_{\mathcal{A},\mathcal{B},\psi}(\mathcal{X})$ $\displaystyle=$
$\displaystyle\sum\limits_{i=1}^{\mathbb{I}_{1}^{N}}\sum\limits_{j=1}^{\mathbb{I}_{1}^{N}}\psi(\lambda_{i},\mu_{j})\mathcal{P}_{\mathcal{U}_{i}}\star_{N}\mathcal{X}\star_{N}\mathcal{P}_{\mathcal{V}_{j}}.$
(6)

###### Lemma 1

Let $\psi,\phi:\mathbb{R}^{2}\rightarrow\mathbb{R}$ be two functions, we have
following relationships about $T_{\mathcal{A},\mathcal{B},\psi}(\mathcal{X})$:
Given $\psi$ is a constant function to one, we have

$\displaystyle T_{\mathcal{A},\mathcal{B},1}=\mathcal{I}.$ (7)

We also have:

$\displaystyle
T_{\mathcal{A},\mathcal{B},\psi\phi}(\mathcal{X})=T_{\mathcal{A},\mathcal{B},\psi}(\mathcal{X})\circ
T_{\mathcal{A},\mathcal{B},\phi}(\mathcal{X}),$ (8)

where $\circ$ is the entrywise product (Hadamard product). Finally, we have

$\displaystyle
T_{\mathcal{A},\mathcal{B},a\psi+b\phi}(\mathcal{X})=aT_{\mathcal{A},\mathcal{B},\psi}(\mathcal{X})+bT_{\mathcal{A},\mathcal{B},\phi}(\mathcal{X}),$
(9)

where $a,b$ are two complex numbers.

Proof:

Given a tensor $\mathcal{X}\in\mathbb{C}^{I_{1}\times\cdots\times I_{N}\times
I_{1}\times\cdots\times I_{N}}$ with orthgonal unitary tensors
$\mathcal{U}_{i}$ and $\mathcal{V}_{j}$ such that
$\mathcal{A}=\sum\limits_{i=1}^{\mathbb{I}_{1}^{N}}\lambda_{i}\mathcal{U}_{i}\star_{1}\mathcal{U}^{H}_{i}$
and
$\mathcal{B}=\sum\limits_{j=1}^{\mathbb{I}_{1}^{N}}\mu_{j}\mathcal{U}_{j}\star_{1}\mathcal{U}^{H}_{j}$,
we define the scaler $x_{i,j}$ assocoated to $\mathcal{X}$ as

$\displaystyle x_{i,j}$ $\displaystyle=$
$\displaystyle\langle\mathcal{X}\star_{N}\mathcal{V}_{j},\mathcal{U}_{i}\rangle.$
(10)

After selecting two specific indices $i^{\prime}$ and $j^{\prime}$, we have

$\displaystyle\langle
T_{\mathcal{A},\mathcal{B},\psi}(\mathcal{X})\star_{N}\mathcal{V}_{j^{\prime}},\mathcal{U}_{i^{\prime}}\rangle$
$\displaystyle=$
$\displaystyle\psi(\lambda_{i^{\prime}},\mu_{j^{\prime}})\langle\mathcal{X}\star_{N}\mathcal{V}_{j^{\prime}},\mathcal{U}_{i^{\prime}}\rangle=\psi(\lambda_{i^{\prime}},\mu_{j^{\prime}})x_{i^{\prime},j^{\prime}}.$
(11)

From Eq. (11) and the Hadamard product properties, we have properties provided
by Eqs. (7), (8) and (9). $\hfill\Box$

###  2.2 Random DTI

According to the DTI definition provided by Eq. (6), the random DTI considered
in this work is to assume that tensors $\mathcal{A},\mathcal{B}$ are random
Hermitian tensors and the remaining parameters $\psi$ and $\mathcal{X}$ are
deterministics. Therefore, we have the randomness at the following terms in
Eq. (6): $\psi(\lambda_{i},\mu_{j}),\mathcal{P}_{\mathcal{U}_{i}}$ and
$\mathcal{P}_{\mathcal{V}_{j}}$. If we are provided more detailed probability
density functions for entries of random Hermitian tensors
$\mathcal{A},\mathcal{B}$, all bounds derived in this work can be improved
with more dedicated expressions associated with parameters of probability
density functions.

##  3 Tail Bound for Random Tensor Integral Norms

###  3.1 Unitarily Invariant Tensor Norms

Let us represent the Hermitian eigenvalues of a Hermitian tensor
$\mathcal{H}\in\mathbb{C}^{I_{1}\times\cdots\times I_{N}\times
I_{1}\times\cdots\times I_{N}}$ in decreasing order by the vector
$\vec{\lambda}(\mathcal{H})=(\lambda_{1}(\mathcal{H}),\cdots,\lambda_{\mathbb{I}_{1}^{N}}(\mathcal{H}))$.
We use $\mathbb{R}_{\geq 0}(\mathbb{R}_{>0})$ to represent a set of
nonnegative (positive) real numbers. Let $\left\|\cdot\right\|_{\rho}$ be a
unitarily invariant tensor norm, i.e.,
$\left\|\mathcal{H}\star_{N}\mathcal{U}\right\|_{\rho}=\left\|\mathcal{U}\star_{N}\mathcal{H}\right\|_{\rho}=\left\|\mathcal{H}\right\|_{\rho}$,
where $\mathcal{U}$ is any unitary tensor. Let $\rho:\mathbb{R}_{\geq
0}^{\mathbb{I}_{1}^{N}}\rightarrow\mathbb{R}_{\geq 0}$ be the corresponding
gauge function that satisfies H$\ddot{o}$lder’s inequality so that

$\displaystyle\left\|\mathcal{H}\right\|_{\rho}=\left\||\mathcal{H}|\right\|_{\rho}=\rho(\vec{\lambda}(|\mathcal{H}|)),$
(12)

where $|\mathcal{H}|\stackrel{{\scriptstyle\mbox{\tiny
def}}}{{=}}\sqrt{\mathcal{H}^{H}\star_{N}\mathcal{H}}$.

We will provide several popular tensor norm examples which can be treated as
special cases of unitarily invariant tensor norm. The first one is Schatten
$p$-norm for tensors, denoted as $\left\|\mathcal{X}\right\|_{p}$, is defined
as:

$\displaystyle\left\|\mathcal{X}\right\|_{p}\stackrel{{\scriptstyle\mbox{\tiny
def}}}{{=}}(\mathrm{Tr}|\mathcal{X}|^{p})^{\frac{1}{p}},$ (13)

where $p\geq 1$. If $p=1$, it is the trace norm.

The second one is $k$-trace norm, denoted as $\mathrm{Tr}_{k}[\mathcal{X}]$,
defined by [17]. It is

$\displaystyle\mathrm{Tr}_{k}[\mathcal{X}]\stackrel{{\scriptstyle\mbox{\tiny
def}}}{{=}}\sum\limits_{1\leq i_{1}<i_{2}<\cdots i_{k}\leq
r}\lambda_{i_{1}}\lambda_{i_{1}}\cdots\lambda_{i_{k}}$ (14)

where $1\leq k\leq r$. If $k=1$, $\mathrm{Tr}_{k}[\mathcal{X}]$ is reduced as
trace norm.

The third one is Ky Fan like $k$-norm [18] for tensors. For
$k\in\\{1,2,\cdots,\mathbb{I}_{1}^{N}\\}$, the Ky Fan $k$-norm [18] for
tensors $\mathcal{X}\in\mathbb{C}^{I_{1}\times\cdots\times I_{N}\times
I_{1}\times\cdots\times I_{N}}$, denoted as
$\left\|\mathcal{X}\right\|_{(k)}$, is defined as:

$\displaystyle\left\|\mathcal{X}\right\|_{(k)}\stackrel{{\scriptstyle\mbox{\tiny
def}}}{{=}}\sum\limits_{i=1}^{k}\lambda_{i}(|\mathcal{X}|).$ (15)

If $k=1$, the Ky Fan $k$-norm for tensors is the tensor operator norm, denoted
as $\left\|\mathcal{X}\right\|$. In this work, we will apply the symbol
$\left\|\mathcal{X}\right\|_{\rho}$ to represent any unitarily invariant
tensor norm for the tensor $\mathcal{X}$.

In the following theorem, we will present the tail bound of unitarily
invariant tensor norm for a given tensor integral
$T_{\mathcal{A},\mathcal{B},\psi}(\mathcal{X})$.

###### Theorem 1

Let $\mathcal{A},\mathcal{B}\in\mathbb{C}^{I_{1}\times\cdots\times I_{N}\times
I_{1}\times\cdots\times I_{N}}$ be independent random Hermitian tensors with
$\mathcal{A}=\sum\limits_{i=1}^{\mathbb{I}_{1}^{N}}\lambda_{i}\mathcal{U}_{i}\star_{1}\mathcal{U}^{H}_{i}$
and
$\mathcal{B}=\sum\limits_{j=1}^{\mathbb{I}_{1}^{N}}\mu_{j}\mathcal{U}_{j}\star_{1}\mathcal{U}^{H}_{j}$,
then, for any $\theta>0$, we have

$\displaystyle\mathrm{Pr}\left(\left\|T_{\mathcal{A},\mathcal{B},\psi}(\mathcal{X})\right\|_{\rho}\geq\theta\right)\leq\frac{\left(\mathbb{I}_{1}^{N}\right)^{2}\left\|\mathcal{X}\right\|_{\rho}}{\theta}\sum\limits_{i=1}^{\mathbb{I}_{1}^{N}}\sum\limits_{j=1}^{\mathbb{I}_{1}^{N}}\mathbb{E}\left(\left|\psi(\lambda_{i},\mu_{j})\right|\right),$
(16)

where $\mathbb{E}$ is the expectation.

Proof:

Since we have the following norm estimation for
$\left\|T_{\mathcal{A},\mathcal{B},\psi}(\mathcal{X})\right\|_{\rho}$:

$\displaystyle\left\|T_{\mathcal{A},\mathcal{B},\psi}(\mathcal{X})\right\|_{\rho}$
$\displaystyle=$
$\displaystyle\left\|\sum\limits_{i=1}^{\mathbb{I}_{1}^{N}}\sum\limits_{j=1}^{\mathbb{I}_{1}^{N}}\psi(\lambda_{i},\mu_{j})\mathcal{P}_{\mathcal{U}_{i}}\star_{N}\mathcal{X}\star_{N}\mathcal{P}_{\mathcal{V}_{j}}\right\|_{\rho}$
(17) $\displaystyle\leq_{1}$
$\displaystyle\sum\limits_{i=1}^{\mathbb{I}_{1}^{N}}\sum\limits_{j=1}^{\mathbb{I}_{1}^{N}}\left|\psi(\lambda_{i},\mu_{j})\right|\left\|\mathcal{P}_{\mathcal{U}_{i}}\star_{N}\mathcal{X}\star_{N}\mathcal{P}_{\mathcal{V}_{j}}\right\|_{\rho}$
$\displaystyle=_{2}$
$\displaystyle\sum\limits_{i=1}^{\mathbb{I}_{1}^{N}}\sum\limits_{j=1}^{\mathbb{I}_{1}^{N}}\left|\psi(\lambda_{i},\mu_{j})\right|\left\|\mathcal{X}\right\|_{\rho},$

where $\leq_{1}$ comes from triangle inequality of the unitarily invariant
norm and $=_{2}$ comes from the definition of the unitarily invariant norm.

Then, we have the following bound for
$\mathrm{Pr}\left(\left\|T_{\mathcal{A},\mathcal{B},\psi}(\mathcal{X})\right\|_{\rho}\geq\theta\right)$

$\displaystyle\mathrm{Pr}\left(\left\|T_{\mathcal{A},\mathcal{B},\psi}(\mathcal{X})\right\|_{\rho}\geq\theta\right)$
$\displaystyle\leq_{1}$
$\displaystyle\mathrm{Pr}\left(\sum\limits_{i=1}^{\mathbb{I}_{1}^{N}}\sum\limits_{j=1}^{\mathbb{I}_{1}^{N}}\left|\psi(\lambda_{i},\mu_{j})\right|\left\|\mathcal{X}\right\|_{\rho}\geq\theta\right)$
(18) $\displaystyle\leq$
$\displaystyle\sum\limits_{i=1}^{\mathbb{I}_{1}^{N}}\sum\limits_{j=1}^{\mathbb{I}_{1}^{N}}\mathrm{Pr}\left(\left|\psi\left(\lambda_{i},\mu_{j}\right)\right|\geq\frac{\theta}{\left\|\mathcal{X}\right\|_{\rho}\left(\mathbb{I}_{1}^{N}\right)^{2}}\right)$
$\displaystyle\leq_{2}$
$\displaystyle\frac{\left\|\mathcal{X}\right\|_{\rho}\left(\mathbb{I}_{1}^{N}\right)^{2}}{\theta}\sum\limits_{i=1}^{\mathbb{I}_{1}^{N}}\sum\limits_{j=1}^{\mathbb{I}_{1}^{N}}\mathbb{E}\left(\left|\psi\left(\lambda_{i},\mu_{j}\right)\right|\right),$

where $\leq_{1}$ comes from the inequality obtained by Eq. (17), and the
$\leq_{2}$ is based on Markov inequality. $\hfill\Box$

We will consider several important examples of $\phi$, which will represent
different tensor means. Given two random Hermitian tensors
$\mathcal{A}=\sum\limits_{i=1}^{\mathbb{I}_{1}^{N}}\lambda_{i}\mathcal{U}_{i}\star_{1}\mathcal{U}^{H}_{i}$
and
$\mathcal{B}=\sum\limits_{j=1}^{\mathbb{I}_{1}^{N}}\mu_{j}\mathcal{U}_{j}\star_{1}\mathcal{U}^{H}_{j}$,
we use $p_{\lambda_{i}}(\leavevmode\nobreak\ )$ and
$p_{\mu_{j}}(\leavevmode\nobreak\ )$ to represent the probability density
functions for eigenvalues $\lambda_{i}$ and $\mu_{j}$, respectively. From Sec.
2.2. in [16], one can apply unfolding technique to convert a random Hermitian
tensor into a Hermitian matrix. If we have Assumption 3.1 in [19], we are able
to obtain the $i$-th eigenvalue distribution of a positive definite Hermitian
tensor, see Corollary 3.3 in [19].

######  Corollary 1 (arithmetic mean)

Under conidtions provided by Theorem 1, if the function $\psi$ has the
following form:

$\displaystyle\psi(x,y)=\frac{x+y}{2},$ (19)

we have

$\displaystyle\mathrm{Pr}\left(\left\|T_{\mathcal{A},\mathcal{B},\psi=\frac{x+y}{2}}(\mathcal{X})\right\|_{\rho}\geq\theta\right)\leq\frac{\left(\mathbb{I}_{1}^{N}\right)^{2}\left\|\mathcal{X}\right\|_{\rho}}{2\theta}\sum\limits_{i=1}^{\mathbb{I}_{1}^{N}}\sum\limits_{j=1}^{\mathbb{I}_{1}^{N}}\left[\mathbb{E}\left(\left|\lambda_{i}\right|\right)+\mathbb{E}\left(\left|\mu_{j}\right|\right)\right],$
(20)

where we have

$\displaystyle\mathbb{E}\left(\left|\lambda_{i}\right|\right)$
$\displaystyle=$
$\displaystyle\int_{0}^{\infty}t\left[p_{\lambda_{i}}(t)+p_{\lambda_{i}}(-t)\right]dt,$
(21)

and

$\displaystyle\mathbb{E}\left(\left|\mu_{j}\right|\right)$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}t\left[p_{\mu_{j}}(t)+p_{\mu_{j}}(-t)\right]dt.$
(22)

Proof: The key is to evaluate
$\mathbb{E}\left(\left|\psi(\lambda_{i},\mu_{j})\right|\right)$, we have

$\displaystyle\mathbb{E}\left(\left|\psi(\lambda_{i},\mu_{j})\right|\right)$
$\displaystyle=$
$\displaystyle\mathbb{E}\left(\left|\frac{\lambda_{i}+\mu_{j}}{2}\right|\right)$
(23) $\displaystyle\leq$
$\displaystyle\frac{1}{2}\mathbb{E}\left(\left|\lambda_{i}\right|\right)+\frac{1}{2}\mathbb{E}\left(\left|\mu_{j}\right|\right).$

This corollary is proved by the following fact for a random variable
$Y=\left|X\right|$:

$\displaystyle F_{Y}(y)=\begin{cases}F_{X}(y)-F_{X}(-y)&y\geq 0,\\\
0&y<0,\end{cases}$ (24)

where $F_{Y}$ and $F_{X}$ are CDFs of random variables $Y$ and $X$.
$\hfill\Box$

######  Corollary 2 (geometric mean)

Under conidtions provided by Theorem 1 with the assumption that $\mathcal{A}$
and $\mathcal{B}$ are random positive definite tensors, if the function $\psi$
has the following form:

$\displaystyle\psi(x,y)=\sqrt{xy},$ (25)

we have

$\displaystyle\mathrm{Pr}\left(\left\|T_{\mathcal{A},\mathcal{B},\psi=\sqrt{xy}}(\mathcal{X})\right\|_{\rho}\geq\theta\right)\leq\frac{\left(\mathbb{I}_{1}^{N}\right)^{2}\left\|\mathcal{X}\right\|_{\rho}}{\theta}\sum\limits_{i=1}^{\mathbb{I}_{1}^{N}}\sum\limits_{j=1}^{\mathbb{I}_{1}^{N}}\left[\mathbb{E}\left(\sqrt{\lambda_{i}}\right)\mathbb{E}\left(\sqrt{\mu_{j}}\right)\right],$
(26)

where we have

$\displaystyle\mathbb{E}\left(\sqrt{\lambda_{i}}\right)$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\sqrt{t}p_{\lambda_{i}}(t)dt,$ (27)

and

$\displaystyle\mathbb{E}\left(\sqrt{\mu_{j}}\right)$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\sqrt{t}p_{\mu_{j}}(t)dt.$ (28)

Proof: To evaluate $\mathbb{E}\left(\sqrt{\lambda_{i}\mu_{j}}\right)$, we have

$\displaystyle\mathbb{E}\left(\left|\psi(\lambda_{i},\mu_{j})\right|\right)$
$\displaystyle=$
$\displaystyle\mathbb{E}\left(\sqrt{\lambda_{i}\mu_{j}}\right)$ (29)
$\displaystyle=$
$\displaystyle\mathbb{E}\left(\sqrt{\lambda_{i}}\right)\mathbb{E}\left(\sqrt{\mu_{j}}\right)$

$\hfill\Box$

######  Corollary 3 (harmonic mean)

Under conidtions provided by Theorem 1 with the assumption that $\mathcal{A}$
and $\mathcal{B}$ are random positive definite tensors, if the function $\psi$
has the following form:

$\displaystyle\psi(x,y)=\frac{2}{x^{-1}+y^{-1}},$ (30)

we have

$\displaystyle\mathrm{Pr}\left(\left\|T_{\mathcal{A},\mathcal{B},\psi=\frac{2}{x^{-1}+y^{-1}}}(\mathcal{X})\right\|_{\rho}\geq\theta\right)$
$\displaystyle\leq$
$\displaystyle\frac{2\left(\mathbb{I}_{1}^{N}\right)^{2}\left\|\mathcal{X}\right\|_{\rho}}{\theta}\sum\limits_{i=1}^{\mathbb{I}_{1}^{N}}\sum\limits_{j=1}^{\mathbb{I}_{1}^{N}}\int_{0}^{\infty}\mathrm{H}_{\lambda_{i}}(t)\mathrm{H}_{\mu_{j}}(t)dt,$
(31)

where we have

$\displaystyle\mathrm{H}_{\lambda_{i}}(t)$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\frac{e^{-tw}}{w^{2}}p_{\lambda_{i}}\left(1/w\right)dw,$
(32)

and

$\displaystyle\mathrm{H}_{\mu_{j}}(t)$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\frac{e^{-tz}}{z^{2}}p_{\mu_{j}}\left(1/z\right)dz.$
(33)

Proof: Let $X$ and $Y$ are two positive random variables with distribution
fnuctions $f_{X}(x)$ and $f_{Y}(y)$. Then, we have

$\displaystyle\mathbb{E}\left(\frac{2}{X^{-1}+Y^{-1}}\right)$
$\displaystyle=_{1}$ $\displaystyle\mathbb{E}\left(\frac{2}{W+Z}\right)$ (34)
$\displaystyle=$ $\displaystyle
2\mathbb{E}\left(\int_{0}^{\infty}\exp\left(-t(W+Z)\right)dt\right)$
$\displaystyle=$ $\displaystyle
2\int_{0}^{\infty}\mathbb{E}\left(\exp\left(-t(W+Z)\right)\right)dt$
$\displaystyle=_{2}$ $\displaystyle
2\int_{0}^{\infty}\mathbb{E}\left(\exp\left(-tW\right)\right)\mathbb{E}\left(\exp\left(-tZ\right)\right)dt$

where we set $W=X^{-1}$ and $Z=Y^{-1}$ at $=_{1}$ and we use indepedent
assumptions of random variables $W=X^{-1}$ and $Z=Y^{-1}$ at $=_{2}$. This
corollary is proved since we have following distribution functions for random
variables $W$ and $Z$, which will correspond to random variables $\lambda_{i}$
and $\mu_{j}$, expressed as

$\displaystyle f(w)=\frac{1}{w^{2}}p_{\lambda_{i}}\left(1/w\right),$ (35)

and

$\displaystyle f(z)=\frac{1}{z^{2}}p_{\mu_{j}}\left(1/z\right).$ (36)

$\hfill\Box$

We have to prepare a lemma about the expectation of ratio between two depedent
random variables before presenting the next corollary.

###### Lemma 2

Given two random variables $X$ and $Y$ such that $Y\neq 0$ always, we have

$\displaystyle\mathbb{E}\left(\frac{X}{Y}\right)=\frac{\mathbb{E}\left(X\right)}{\mathbb{E}\left(Y\right)}+\lim\limits_{\epsilon\rightarrow
0}\sum\limits_{i=1}^{\infty}\frac{(-1)^{i}\left[\mathbb{E}(X)\mathbb{E}(\acute{Y}^{i})+\mathbb{E}(\acute{X}\acute{Y}^{i})\right]}{\prod\limits_{j=0}^{i}(\mathbb{E}(Y)+j\epsilon)},$
(37)

where we define the following random variables:

$\displaystyle\acute{X}$ $\displaystyle\stackrel{{\scriptstyle\mbox{\tiny
def}}}{{=}}$ $\displaystyle X-\mathbb{E}(X),$ $\displaystyle\acute{Y}$
$\displaystyle\stackrel{{\scriptstyle\mbox{\tiny def}}}{{=}}$ $\displaystyle
Y-\mathbb{E}(Y).$

Proof:

We have the following expression about $\mathbb{E}(\frac{X}{Y})$, it is

$\displaystyle\mathbb{E}(\frac{X}{Y})$ $\displaystyle=$
$\displaystyle\mathbb{E}\left(\frac{\mathbb{E}(X)}{\mathbb{E}(Y)}\times\frac{1+\frac{\acute{X}}{\mathbb{E}(X)}}{1+\frac{\acute{Y}}{\mathbb{E}(Y)}}\right)$
(39) $\displaystyle=$
$\displaystyle\frac{\mathbb{E}(X)}{\mathbb{E}(Y)}\mathbb{E}\left(\left(1+\frac{\acute{X}}{\mathbb{E}(X)}\right)\left(1+\frac{\acute{Y}}{\mathbb{E}(Y)}\right)^{-1}\right)$
$\displaystyle=$
$\displaystyle\frac{\mathbb{E}(X)}{\mathbb{E}(Y)}\mathbb{E}\left(\left(1+\frac{\acute{Y}}{\mathbb{E}(Y)}\right)^{-1}\right)+\frac{1}{\mathbb{E}(Y)}\mathbb{E}\left(\acute{X}\left(1+\frac{\acute{Y}}{\mathbb{E}(Y)}\right)^{-1}\right)$

Given a real function $g$, we have the following approximation form:

$\displaystyle g(a+x)=\lim\limits_{\epsilon\rightarrow
0}\sum\limits_{i=0}^{\infty}\frac{x^{i}\Delta_{\epsilon}^{i}g(a)}{\epsilon^{i}i!},$
(40)

where $\Delta_{\epsilon}^{i}$ is the finite difference operator of degree $i$
and the step size $\epsilon$ is defined as

$\displaystyle\Delta_{\epsilon}^{i}g(a)=\sum\limits_{j=0}^{i}(-1)^{j}{i\choose
j}g(a+(i-j)\epsilon).$ (41)

If we apply Eq. (40) to the following function at $\acute{Y}=0$

$\displaystyle
g(\acute{Y})=\left(1+\frac{\acute{Y}}{\mathbb{E}(Y)}\right)^{-1},$ (42)

we will get

$\displaystyle\left(1+\frac{\acute{Y}}{\mathbb{E}(Y)}\right)^{-1}$
$\displaystyle=$ $\displaystyle\lim\limits_{\epsilon\rightarrow
0}\sum\limits_{i=0}^{\infty}\frac{(-1)^{i}(\acute{Y})^{i}\mathbb{E}(Y)}{\prod\limits_{j=0}^{i}(\mathbb{E}(Y)+j\epsilon)}.$
(43)

This lemma is proved by applying Eq (43) to Eq. (39). $\hfill\Box$

Following corollary is about the unitarily invariant norm tail bound for the
two tensors mean general format. Note that we have logarithmic mean if
$\alpha\rightarrow 1$ in Eq. (4).

######  Corollary 4 (general mean)

Under conidtions provided by Theorem 1 with the assumption that $\mathcal{A}$
and $\mathcal{B}$ are random positive definite tensors, if the function $\psi$
has the following form:

$\displaystyle\psi(x,y)=\frac{\alpha-1}{\alpha}\frac{x^{\alpha}-y^{\alpha}}{x^{\alpha-1}-y^{\alpha-1}},$
(44)

where $\alpha\in\mathbb{R}$ and $x\neq y$ 111If $x=y$, this situation has
measure zero.. Then we have

$\displaystyle\mathrm{Pr}\left(\left\|T_{\mathcal{A},\mathcal{B},\psi=\frac{\alpha-1}{\alpha}\frac{x^{\alpha}-y^{\alpha}}{x^{\alpha-1}-y^{\alpha-1}}}(\mathcal{X})\right\|_{\rho}\geq\theta\right)\leq\frac{\left(\mathbb{I}_{1}^{N}\right)^{2}\left\|\mathcal{X}\right\|_{\rho}(\alpha-1)}{\theta\alpha}\times\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ $
$\displaystyle\sum\limits_{i=1}^{\mathbb{I}_{1}^{N}}\sum\limits_{j=1}^{\mathbb{I}_{1}^{N}}\left\\{\frac{\mathbb{E}\left(X_{\lambda_{i},\mu_{j}}\right)}{\mathbb{E}\left(Y_{\lambda_{i},\mu_{j}}\right)}+\lim\limits_{\epsilon\rightarrow
0}\sum\limits_{\i=1}^{\infty}\frac{(-1)^{\i}\left[\mathbb{E}(X_{\lambda_{i},\mu_{j}})\mathbb{E}((\acute{Y}_{\lambda_{i},\mu_{j}})^{\i})+\mathbb{E}(\acute{X}_{\lambda_{i},\mu_{j}}(\acute{Y}_{\lambda_{i},\mu_{j}})^{\i})\right]}{\prod\limits_{\j=0}^{\i}(\mathbb{E}(Y_{\lambda_{i},\mu_{j}})+\j\epsilon)}\right\\},$
(45)

where we have random variables $X_{\lambda_{i},\mu_{j}}$ and
$Y_{\lambda_{i},\mu_{j}}$ defined by

$\displaystyle X_{\lambda_{i},\mu_{j}}$ $\displaystyle=$
$\displaystyle\lambda^{\alpha}_{i}-\mu^{\alpha}_{j},$ (46)

and

$\displaystyle Y_{\lambda_{i},\mu_{j}}$ $\displaystyle=$
$\displaystyle\lambda^{\alpha-1}_{i}-\mu^{\alpha-1}_{j}.$ (47)

Proof:

Since $\left|\psi(\lambda_{i},\mu_{j})\right|$ is

$\displaystyle\psi(\lambda_{i},\mu_{j})$ $\displaystyle=$
$\displaystyle\frac{\alpha-1}{\alpha}\frac{\lambda^{\alpha}_{i}-\mu^{\alpha}_{j}}{\lambda^{\alpha-1}_{i}-\mu^{\alpha-1}_{j}}=\frac{\alpha-1}{\alpha}\frac{X_{\lambda_{i},\mu_{j}}}{Y_{\lambda_{i},\mu_{j}}},$
(48)

this corollary is proved by applying Lemma 2 to the expectation of Eq. (48).
$\hfill\Box$

Note that each of the following terms
$\mathbb{E}\left(X_{\lambda_{i},\mu_{j}}\right),\mathbb{E}\left(Y_{\lambda_{i},\mu_{j}}\right),\mathbb{E}((\acute{Y}_{\lambda_{i},\mu_{j}})^{\i})$
and
$\mathbb{E}(\acute{X}_{\lambda_{i},\mu_{j}}(\acute{Y}_{\lambda_{i},\mu_{j}})^{\i})$
can be evaluated exactly since we know all density distributions
$p_{\lambda_{i}}(\leavevmode\nobreak\ )$ and $p_{\mu_{j}}(\leavevmode\nobreak\
)$.

##  4 Tail Bounds for Random Lipschitz Estimates

In this section, we will try to provide tail bounds for the Lipschitz estimate
for the unitarily invariant norm for a given function, which is the main
result of this section. We will begin with the perturbation lemma. The vanilla
case is discussed in Section 4.1. The case about considering quasi-commutator
is provided by Section 4.2

###  4.1 Vanilla Case

We will begin by providing a perturbation formula for DTI.

###### Lemma 3

Let $\mathcal{A},\mathcal{B}\in\mathbb{C}^{I_{1}\times\cdots\times I_{N}\times
I_{1}\times\cdots\times I_{N}}$ be Hermitian tensors with
$\mathcal{A}=\sum\limits_{i=1}^{\mathbb{I}_{1}^{N}}\lambda_{i}\mathcal{U}_{i}\star_{1}\mathcal{U}^{H}_{i}$
and
$\mathcal{B}=\sum\limits_{j=1}^{\mathbb{I}_{1}^{N}}\mu_{j}\mathcal{U}_{j}\star_{1}\mathcal{U}^{H}_{j}$.
Also, let $f:\mathbb{R}\rightarrow\mathbb{R}$ such that $f^{\prime}(x)$
exists, we define the bivariate function $f^{[1]}$ as

$\displaystyle f^{[1]}(x,y)=\begin{cases}\frac{f(x)-f(y)}{x-y},&x\neq y;\\\
f^{\prime}(x),&x=y.\end{cases}$ (49)

Then, we have

$\displaystyle
f(\mathcal{A})-f(\mathcal{B})=T_{\mathcal{A},\mathcal{B},f^{[1]}}\left(\mathcal{A}-\mathcal{B}\right).$
(50)

Proof:

Since we have

$\displaystyle\mathcal{A}=\sum\limits_{i=1}^{\mathbb{I}_{1}^{N}}\lambda_{i}\mathcal{U}_{i}\star_{1}\mathcal{U}^{H}_{i},$
(51)

and

$\displaystyle\mathcal{B}=\sum\limits_{j=1}^{\mathbb{I}_{1}^{N}}\mu_{j}\mathcal{U}_{j}\star_{1}\mathcal{U}^{H}_{j},$
(52)

then, we will obtain the following

$\displaystyle
f(\mathcal{A})=\sum\limits_{i=1}^{\mathbb{I}_{1}^{N}}f(\lambda_{i})\mathcal{U}_{i}\star_{1}\mathcal{U}^{H}_{i},$
(53)

and

$\displaystyle
f(\mathcal{B})=\sum\limits_{j=1}^{\mathbb{I}_{1}^{N}}f(\mu_{j})\mathcal{U}_{j}\star_{1}\mathcal{U}^{H}_{j}.$
(54)

Because, we also have

$\displaystyle
f(\mathcal{A})\star_{N}\mathcal{U}_{i}=f(\lambda_{i})\star_{N}\mathcal{U}_{i},$
(55)

and

$\displaystyle
f(\mathcal{B})\star_{N}\mathcal{V}_{j}=f(\mu_{j})\star_{N}\mathcal{V}_{j},$
(56)

then, we have

$\displaystyle\langle\left(f(\mathcal{A})-f(\mathcal{B})\star_{N}\mathcal{V}_{j},\mathcal{U}_{i}\right)\rangle$
$\displaystyle=$ $\displaystyle
f^{[1]}(\lambda_{i},\mu_{j})\langle\left(\mathcal{A}-\mathcal{B}\right)\star_{N}\mathcal{V}_{j},\mathcal{U}_{i}\rangle.$
(57)

By applying Eq. (6) definition and Eq. (11) to Eq. (57), this Lemma is proved.
$\hfill\Box$

From perturbation formula given by Lemma 3, we can have the following theoem
about the tail bounds of the unitarily invariant norm for the Lipschitz
estimate of tensor-valued function with random tensors as inputs.

###### Theorem 2

Let $\mathcal{A},\mathcal{B}\in\mathbb{C}^{I_{1}\times\cdots\times I_{N}\times
I_{1}\times\cdots\times I_{N}}$ be independent random Hermitian tensors with
$\mathcal{A}=\sum\limits_{i=1}^{\mathbb{I}_{1}^{N}}\lambda_{i}\mathcal{U}_{i}\star_{1}\mathcal{U}^{H}_{i}$
and
$\mathcal{B}=\sum\limits_{j=1}^{\mathbb{I}_{1}^{N}}\mu_{j}\mathcal{U}_{j}\star_{1}\mathcal{U}^{H}_{j}$,
moreover, we are given a real valud funtion $f(x)$ for $x\in\mathbb{R}$ such
that $\left|f^{[1]}(x,y)\right|$ is bounded by a positive number denoted as
$\Omega_{f^{[1]}}$. Then, for any $\theta>0$, we have

$\displaystyle\mathrm{Pr}\left(\left\|f(\mathcal{A})-f(\mathcal{B})\right\|_{\rho}\geq\theta\right)\leq\frac{\left(\mathbb{I}_{1}^{N}\right)^{2}\Omega_{f^{[1]}}}{\theta}\mathbb{E}\left(\left\|\mathcal{A}-\mathcal{B}\right\|_{\rho}\right).$
(58)

Proof:

We have

$\displaystyle\left\|f(\mathcal{A})-f(\mathcal{B})\right\|_{\rho}$
$\displaystyle=$
$\displaystyle\left\|T_{\mathcal{A},\mathcal{B},f^{[1]}}\left(\mathcal{A}-\mathcal{B}\right)\right\|_{\rho}$
(59) $\displaystyle=$
$\displaystyle\left\|\sum\limits_{i=1}^{\mathbb{I}_{1}^{N}}\sum\limits_{j=1}^{\mathbb{I}_{1}^{N}}f^{[1]}(\lambda_{i},\mu_{j})\mathcal{P}_{\mathcal{U}_{i}}\star_{N}\left(\mathcal{A}-\mathcal{B}\right)\star_{N}\mathcal{P}_{\mathcal{V}_{j}}\right\|_{\rho}$
$\displaystyle\leq_{1}$
$\displaystyle\sum\limits_{i=1}^{\mathbb{I}_{1}^{N}}\sum\limits_{j=1}^{\mathbb{I}_{1}^{N}}\left|f^{[1]}(\lambda_{i},\mu_{j})\right|\left\|\mathcal{P}_{\mathcal{U}_{i}}\star_{N}\left(\mathcal{A}-\mathcal{B}\right)\star_{N}\mathcal{P}_{\mathcal{V}_{j}}\right\|_{\rho}$
$\displaystyle=_{2}$
$\displaystyle\sum\limits_{i=1}^{\mathbb{I}_{1}^{N}}\sum\limits_{j=1}^{\mathbb{I}_{1}^{N}}\left|f^{[1]}(\lambda_{i},\mu_{j})\right|\left\|\mathcal{A}-\mathcal{B}\right\|_{\rho}$
$\displaystyle\leq_{3}$
$\displaystyle\left(\mathbb{I}_{1}^{N}\right)^{2}\Omega_{f^{[1]}}\left\|\mathcal{A}-\mathcal{B}\right\|_{\rho},$

where $\leq_{1}$ comes from triangle inequality of the unitarily invariant
norm, $=_{2}$ comes from the definition of the unitarily invariant norm, and
$\leq_{3}$ is due to that
$\left|f^{[1]}(\lambda_{i},\mu_{j})\right|\leq\Omega_{f^{[1]}}$ (mean value
theorem of divide difference).

Then, we have the following bound for
$\mathrm{Pr}\left(\left\|f(\mathcal{A})-f(\mathcal{B})\right\|_{\rho}\geq\theta\right)$

$\displaystyle\mathrm{Pr}\left(\left\|f(\mathcal{A})-f(\mathcal{B})\right\|_{\rho}\geq\theta\right)$
$\displaystyle\leq_{1}$
$\displaystyle\mathrm{Pr}\left(\left(\mathbb{I}_{1}^{N}\right)^{2}\Omega_{f^{[1]}}\left\|\mathcal{A}-\mathcal{B}\right\|_{\rho}\geq\theta\right)$
(60) $\displaystyle\leq_{2}$
$\displaystyle\frac{\left(\mathbb{I}_{1}^{N}\right)^{2}\Omega_{f^{[1]}}}{\theta}\mathbb{E}\left(\left|\mathcal{A}-\mathcal{B}\right|_{\rho}\right),$

where $\leq_{1}$ comes from the inequality obtained by Eq. (59), and the
$\leq_{2}$ is based on Markov inequality. $\hfill\Box$

Since the upper bound $\Omega_{f^{[1]}}$ depends on $f(x)$, we will consider
the following two corollaries about special types of the function $f(x)$.

######  Corollary 5 (Lipschitz estimate for polynomial fnuctions)

Let $\mathcal{A},\mathcal{B}\in\mathbb{C}^{I_{1}\times\cdots\times I_{N}\times
I_{1}\times\cdots\times I_{N}}$ be independent random Hermitian tensors with
$\mathcal{A}=\sum\limits_{i=1}^{\mathbb{I}_{1}^{N}}\lambda_{i}\mathcal{U}_{i}\star_{1}\mathcal{U}^{H}_{i}$
and
$\mathcal{B}=\sum\limits_{j=1}^{\mathbb{I}_{1}^{N}}\mu_{j}\mathcal{U}_{j}\star_{1}\mathcal{U}^{H}_{j}$.
Besides, we are given a real polynomial funtion $f(x)$ with degree $m$ over an
interval $[a,b]$. For any $\theta>0$, we have

$\displaystyle\mathrm{Pr}\left(\left\|f(\mathcal{A})-f(\mathcal{B})\right\|_{\rho}\geq\theta\right)\leq\frac{\left(\mathbb{I}_{1}^{N}\right)^{2}\sum\limits_{k=1}^{m}\frac{\left|f^{(k)}(x^{*}_{(k)})\right|(b-a)^{k}}{k!}}{\theta}\mathbb{E}\left(\left\|\mathcal{A}-\mathcal{B}\right\|_{\rho}\right),$
(61)

where $x^{*}_{(k)}$ is the maximizer to reach the maximum value for the
function $\left|f^{(k)}(x)\right|$, i.e., the absolute value of the $k$-th
derivative, in the interval $[a,b]$.

Proof:

If we perform Tayler expansion for the function $f$ at $x$, we have

$\displaystyle
f(y)=f(x)+f^{\prime}(x)(y-x)+f^{\prime\prime}(x)\frac{(y-x)^{2}}{2!}+f^{\prime\prime\prime}(x)\frac{(y-x)^{3}}{3!}+\ldots,$
(62)

which is equivalent to have

$\displaystyle\frac{f(y)-f(x)}{y-x}=f^{\prime}(x)+f^{\prime\prime}(x)\frac{(y-x)}{2!}+f^{\prime\prime\prime}(x)\frac{(y-x)^{2}}{3!}+\ldots.$
(63)

Because the polynomial function $f$ has degree $m$, from Eq. (62), we have the
following bound from triangle inequality:

$\displaystyle\left|\frac{f(y)-f(x)}{y-x}\right|$ $\displaystyle=$
$\displaystyle\left|\sum\limits_{k=1}^{m}\frac{f^{(k)}(x)(y-x)^{k}}{k!}\right|$
(64) $\displaystyle\leq$
$\displaystyle\sum\limits_{k=1}^{m}\frac{\left|f^{(k)}(x^{*}_{(k)})\right|(b-a)^{k}}{k!},$

where $x^{*}_{(k)}$ is the maximizer to reach the maximum value for the
function $\left|f^{(k)}(x^{*}_{(k)})\right|$ in the interval $[a,b]$. This
corrollary is proved by Theorem 2. $\hfill\Box$

Following corollary is about Lipschitz estimate for polygamma functions.
Recall that a digamma function $\omega(x)$ is defined as

$\displaystyle\omega(x)=\frac{\Gamma^{\prime}(x)}{\Gamma(x)},$ (65)

where $\Gamma(x)=\int_{0}^{\infty}e^{-t}t^{x-1}dt$ [20]. Then polygamma
functions are defined as the $k$-th derivative $\omega^{(k)}(x)$ for any
$k\in\mathbb{N}$.

######  Corollary 6 (Lipschitz estimate for polygamma functions)

Let $\mathcal{A},\mathcal{B}\in\mathbb{C}^{I_{1}\times\cdots\times I_{N}\times
I_{1}\times\cdots\times I_{N}}$ be independent random positive definite
tensors with
$\mathcal{A}=\sum\limits_{i=1}^{\mathbb{I}_{1}^{N}}\lambda_{i}\mathcal{U}_{i}\star_{1}\mathcal{U}^{H}_{i}$
and
$\mathcal{B}=\sum\limits_{j=1}^{\mathbb{I}_{1}^{N}}\mu_{j}\mathcal{U}_{j}\star_{1}\mathcal{U}^{H}_{j}$.
Besides, we are given a polygamma funtion $\omega^{(k)}(x)$ for any
$k\in\mathbb{N}$ and $x>0$. For any $\theta>0$, we have

$\displaystyle\mathrm{Pr}\left(\left\|\omega^{(k)}(\mathcal{A})-\omega^{(k)}(\mathcal{B})\right\|_{\rho}\geq\theta\right)\leq\frac{\left(\mathbb{I}_{1}^{N}\right)^{2}\Omega_{\omega^{[k+1]}}}{\theta}\mathbb{E}\left(\left\|\mathcal{A}-\mathcal{B}\right\|_{\rho}\right),$
(66)

where

$\displaystyle\Omega_{\omega^{[k+1]}}$ $\displaystyle=$
$\displaystyle\max\left\\{\omega^{(k+1)}\left(\frac{1}{e}\left(\frac{(b^{*})^{b^{*}}}{(a^{*})^{a^{*}}}\right)^{\frac{1}{b^{*}-a^{*}}}\right),\omega^{(k+1)}(x^{*})\right\\}.$
(67)

The values $a^{*}$ and $b^{*}$ are the maximizers of the function
$\omega^{(k+1)}\left(\frac{1}{e}\left(\frac{b^{b}}{a^{a}}\right)^{\frac{1}{b-a}}\right)$
given $a\neq b$ and $a,b>0$. The value $x^{*}$ is the maximizer of the
function $\omega^{(k+1)}(x)$ given $x>0$.

Proof:

This corrollary is proved according to Theorem 2 by applying Theorem 1 from
[20] to Eq. (65). $\hfill\Box$

###  4.2 Quasi-Commutator Case

In this section, we will extend the perturbation formula provided by Lemma 4
to the quasi-commutator
$\mathcal{D}\star\mathcal{A}-\mathcal{B}\star\mathcal{D}$. Tail bouns based on
the quasi-commutator $\mathcal{D}\star\mathcal{A}-\mathcal{B}\star\mathcal{D}$
will be given in this section.

###### Lemma 4

Let $\mathcal{A},\mathcal{B}\in\mathbb{C}^{I_{1}\times\cdots\times I_{N}\times
I_{1}\times\cdots\times I_{N}}$ be random Hermitian tensors with
$\mathcal{A}=\sum\limits_{i=1}^{\mathbb{I}_{1}^{N}}\lambda_{i}\mathcal{U}_{i}\star_{1}\mathcal{U}^{H}_{i}$
and
$\mathcal{B}=\sum\limits_{j=1}^{\mathbb{I}_{1}^{N}}\mu_{j}\mathcal{V}_{j}\star_{1}\mathcal{V}^{H}_{j}$.
Besides, we have the function $f:\mathbb{R}\rightarrow\mathbb{R}$ such that
$f^{\prime}(x)$ exists, then, we have

$\displaystyle\mathcal{D}\star_{N}f(\mathcal{A})-f(\mathcal{B})\star_{N}\mathcal{D}=T_{\mathcal{A},\mathcal{B},f^{[1]}}\left(\mathcal{D}\star_{N}\mathcal{A}-\mathcal{B}\star_{N}\mathcal{D}\right).$
(68)

where $f^{[1]}$ has been defined by Eq. (49).

Proof:

Since we have

$\displaystyle\mathcal{A}=\sum\limits_{i=1}^{\mathbb{I}_{1}^{N}}\lambda_{i}\mathcal{U}_{i}\star_{1}\mathcal{U}^{H}_{i},$
(69)

and

$\displaystyle\mathcal{B}=\sum\limits_{j=1}^{\mathbb{I}_{1}^{N}}\mu_{j}\mathcal{V}_{j}\star_{1}\mathcal{V}^{H}_{j},$
(70)

then, we will obtain the following

$\displaystyle\mathcal{D}\star_{N}\mathcal{A}=T_{\mathcal{A},\mathcal{I},\lambda_{i}}\star_{N}\mathcal{D},$
(71)

and

$\displaystyle\mathcal{B}\star_{N}\mathcal{D}=T_{\mathcal{I},\mathcal{B},\mu_{j}}\star_{N}\mathcal{D}.$
(72)

By applying Lemma 1, we have

$\displaystyle
T_{\mathcal{A},\mathcal{B},f^{[1]}}\left(\mathcal{D}\star_{N}\mathcal{A}-\mathcal{B}\star_{N}\mathcal{D}\right)$
$\displaystyle=$ $\displaystyle
T_{\mathcal{A},\mathcal{B},f^{[1]}}\left(T_{\mathcal{A},\mathcal{I},\lambda_{i}}\star_{N}\mathcal{D}-T_{\mathcal{I},\mathcal{B},\mu_{j}}\star_{N}\mathcal{D}\right)$
(73) $\displaystyle=$ $\displaystyle
T_{\mathcal{A},\mathcal{B},f^{[1]}(\lambda_{i},\mu_{j})}\left(\mathcal{D}\right)$
$\displaystyle=$ $\displaystyle
T_{\mathcal{A},\mathcal{B},f(\lambda_{i})-f(\mu_{j})}\left(\mathcal{D}\right)$
$\displaystyle=$
$\displaystyle\mathcal{D}\star_{N}f(\mathcal{A})-f(\mathcal{B})\star_{N}\mathcal{D}.$

$\hfill\Box$

Following theorem is the tail bound for
$\mathcal{D}\star_{N}f(\mathcal{A})-f(\mathcal{B})\star_{N}\mathcal{D}$.

###### Theorem 3

Let $\mathcal{A},\mathcal{B}\in\mathbb{C}^{I_{1}\times\cdots\times I_{N}\times
I_{1}\times\cdots\times I_{N}}$ be independent random Hermitian tensors with
$\mathcal{A}=\sum\limits_{i=1}^{\mathbb{I}_{1}^{N}}\lambda_{i}\mathcal{U}_{i}\star_{1}\mathcal{U}^{H}_{i}$
and
$\mathcal{B}=\sum\limits_{j=1}^{\mathbb{I}_{1}^{N}}\mu_{j}\mathcal{V}_{j}\star_{1}\mathcal{V}^{H}_{j}$,
moreover, we are given a real valud funtion $f(x)$ for $x\in\mathbb{R}$ such
that $f^{\prime}(x)$ exists and $\left|f^{\prime}(x)\right|$ is bounded by a
positive number denoted as $\Omega_{f^{[1]}}$. Then, for any $\theta>0$, we
have

$\displaystyle\mathrm{Pr}\left(\left\|\mathcal{D}\star_{N}f(\mathcal{A})-f(\mathcal{B})\star_{N}\mathcal{D}\right\|_{\rho}\geq\theta\right)\leq\frac{\left(\mathbb{I}_{1}^{N}\right)^{2}\Omega_{f^{[1]}}}{\theta}\mathbb{E}\left(\left\|\mathcal{D}\star_{N}\mathcal{A}-\mathcal{B}\star_{N}\mathcal{D}\right\|_{\rho}\right).$
(74)

Proof:

From Lemma 4, we have

$\displaystyle\left\|\mathcal{D}\star_{N}f(\mathcal{A})-f(\mathcal{B})\star_{N}\mathcal{D}\right\|_{\rho}$
$\displaystyle=$
$\displaystyle\left\|T_{\mathcal{A},\mathcal{B},f^{[1]}}\left(\mathcal{D}\star_{N}\mathcal{A}-\mathcal{B}\star_{N}\mathcal{D}\right)\right\|_{\rho}$
(75) $\displaystyle=$
$\displaystyle\left\|\sum\limits_{i=1}^{\mathbb{I}_{1}^{N}}\sum\limits_{j=1}^{\mathbb{I}_{1}^{N}}f^{[1]}(\lambda_{i},\mu_{j})\mathcal{P}_{\mathcal{U}_{i}}\star_{N}\left(\mathcal{D}\star_{N}\mathcal{A}-\mathcal{B}\star_{N}\mathcal{D}\right)\star_{N}\mathcal{P}_{\mathcal{V}_{j}}\right\|_{\rho}$
$\displaystyle\leq_{1}$
$\displaystyle\sum\limits_{i=1}^{\mathbb{I}_{1}^{N}}\sum\limits_{j=1}^{\mathbb{I}_{1}^{N}}\left|f^{[1]}(\lambda_{i},\mu_{j})\right|\left\|\mathcal{P}_{\mathcal{U}_{i}}\star_{N}\left(\mathcal{D}\star_{N}\mathcal{A}-\mathcal{B}\star_{N}\mathcal{D}\right)\star_{N}\mathcal{P}_{\mathcal{V}_{j}}\right\|_{\rho}$
$\displaystyle=_{2}$
$\displaystyle\sum\limits_{i=1}^{\mathbb{I}_{1}^{N}}\sum\limits_{j=1}^{\mathbb{I}_{1}^{N}}\left|f^{[1]}(\lambda_{i},\mu_{j})\right|\left\|\mathcal{D}\star_{N}\mathcal{A}-\mathcal{B}\star_{N}\mathcal{D}\right\|_{\rho}$
$\displaystyle\leq_{3}$
$\displaystyle\left(\mathbb{I}_{1}^{N}\right)^{2}\Omega_{f^{[1]}}\left\|\mathcal{D}\star_{N}\mathcal{A}-\mathcal{B}\star_{N}\mathcal{D}\right\|_{\rho},$

where $\leq_{1}$ comes from triangle inequality of the unitarily invariant
norm, $=_{2}$ comes from the definition of the unitarily invariant norm, and
$\leq_{3}$ is due to that
$\left|f^{[1]}(\lambda_{i},\mu_{j})\right|\leq\Omega_{f^{[1]}}$ (mean value
theorem of divide difference).

Then, we have the following bound for
$\mathrm{Pr}\left(\left\|\mathcal{D}\star_{N}f(\mathcal{A})-f(\mathcal{B})\star_{N}\mathcal{D}\right\|_{\rho}\geq\theta\right)$

$\displaystyle\mathrm{Pr}\left(\left\|\mathcal{D}\star_{N}f(\mathcal{A})-f(\mathcal{B})\star_{N}\mathcal{D}\right\|_{\rho}\geq\theta\right)$
$\displaystyle\leq_{1}$
$\displaystyle\mathrm{Pr}\left(\left(\mathbb{I}_{1}^{N}\right)^{2}\Omega_{f^{[1]}}\left\|\mathcal{D}\star_{N}\mathcal{A}-\mathcal{B}\star_{N}\mathcal{D}\right\|_{\rho}\geq\theta\right)$
(76) $\displaystyle\leq_{2}$
$\displaystyle\frac{\left(\mathbb{I}_{1}^{N}\right)^{2}\Omega_{f^{[1]}}}{\theta}\mathbb{E}\left(\left\|\mathcal{D}\star_{N}\mathcal{A}-\mathcal{B}\star_{N}\mathcal{D}\right\|_{\rho}\right),$

where $\leq_{1}$ comes from the inequality obtained by Eq. (75), and the
$\leq_{2}$ is based on Markov inequality. $\hfill\Box$

Similar to Corollaries 5 and 6, we have following two corollaries for
Lipschitz estimate for polynomial and polygamma under quasi-commutator case.

######  Corollary 7 (Lipschitz estimate for polynomial fnuctions, quasi-
commutator case)

Let $\mathcal{A},\mathcal{B}\in$
$\mathbb{C}^{I_{1}\times\cdots\times I_{N}\times I_{1}\times\cdots\times
I_{N}}$ be independent random Hermitian tensors with
$\mathcal{A}=\sum\limits_{i=1}^{\mathbb{I}_{1}^{N}}\lambda_{i}\mathcal{U}_{i}\star_{1}\mathcal{U}^{H}_{i}$
and
$\mathcal{B}=\sum\limits_{j=1}^{\mathbb{I}_{1}^{N}}\mu_{j}\mathcal{V}_{j}\star_{1}\mathcal{V}^{H}_{j}$.
Besides, we are given a real polynomial funtion $f(x)$ with degree $m$ over an
interval $[a,b]$. For any $\theta>0$, we have

$\displaystyle\mathrm{Pr}\left(\left\|\mathcal{D}\star_{N}f(\mathcal{A})-f(\mathcal{B})\star_{N}\mathcal{D}\right\|_{\rho}\geq\theta\right)\leq\frac{\left(\mathbb{I}_{1}^{N}\right)^{2}\sum\limits_{k=1}^{m}\frac{\left|f^{(k)}(x^{*}_{(k)})\right|(b-a)^{k}}{k!}}{\theta}\mathbb{E}\left(\left\|\mathcal{D}\star_{N}\mathcal{A}-\mathcal{B}\star_{N}\mathcal{D}\right\|_{\rho}\right),$
(77)

where $x^{*}_{(k)}$ is the maximizer to reach the maximum value for the
function $\left|f^{(k)}(x^{*}_{(k)})\right|$ in the interval $[a,b]$.

Proof:

If we perform Tayler expansion for the function $f$ at $x$, we have

$\displaystyle
f(y)=f(x)+f^{\prime}(x)(y-x)+f^{\prime\prime}(x)\frac{(y-x)^{2}}{2!}+f^{\prime\prime\prime}(x)\frac{(y-x)^{3}}{3!}+\ldots,$
(78)

which is equivalent to have

$\displaystyle\frac{f(y)-f(x)}{y-x}=f^{\prime}(x)+f^{\prime\prime}(x)\frac{(y-x)}{2!}+f^{\prime\prime\prime}(x)\frac{(y-x)^{2}}{3!}+\ldots.$
(79)

Because the polynomial function $f$ has degree $m$, from Eq. (62), we have the
following bound from triangle inequality:

$\displaystyle\left|\frac{f(y)-f(x)}{y-x}\right|$ $\displaystyle=$
$\displaystyle\left|\sum\limits_{k=1}^{m}\frac{f^{(k)}(x)(y-x)^{k}}{k!}\right|$
(80) $\displaystyle\leq$
$\displaystyle\sum\limits_{k=1}^{m}\frac{\left|f^{(k)}(x^{*}_{(k)})\right|(b-a)^{k}}{k!},$

where $x^{*}_{(k)}$ is the maximizer to reach the maximum value for the
function $\left|f^{(k)}(x^{*}_{(k)})\right|$ in the interval $[a,b]$. This
corrollary is proved by Theorem 3. $\hfill\Box$

######  Corollary 8 (Lipschitz estimate for polygamma functions: quasi-
commutator case)

Let $\mathcal{A},\mathcal{B}\in$
$\mathbb{C}^{I_{1}\times\cdots\times I_{N}\times I_{1}\times\cdots\times
I_{N}}$ be independent random positive definite tensors with
$\mathcal{A}=\sum\limits_{i=1}^{\mathbb{I}_{1}^{N}}\lambda_{i}\mathcal{U}_{i}\star_{1}\mathcal{U}^{H}_{i}$
and
$\mathcal{B}=\sum\limits_{j=1}^{\mathbb{I}_{1}^{N}}\mu_{j}\mathcal{V}_{j}\star_{1}\mathcal{V}^{H}_{j}$.
Besides, we are given a polygamma funtion $\omega^{(k)}(x)$ for any
$k\in\mathbb{N}$ and $x>0$. For any $\theta>0$, we have

$\displaystyle\mathrm{Pr}\left(\left\|\mathcal{D}\star_{N}\omega^{(k)}(\mathcal{A})-\omega^{(k)}(\mathcal{B})\star_{N}\mathcal{D}\right\|_{\rho}\geq\theta\right)\leq\frac{\left(\mathbb{I}_{1}^{N}\right)^{2}\Omega_{\omega^{[k+1]}}}{\theta}\mathbb{E}\left(\left\|\mathcal{D}\star_{N}\mathcal{A}-\mathcal{B}\star_{N}\mathcal{D}\right\|_{\rho}\right),$
(81)

where

$\displaystyle\Omega_{\omega^{[k+1]}}$ $\displaystyle=$
$\displaystyle\max\left\\{\omega^{(k+1)}\left(\frac{1}{e}\left(\frac{(b^{*})^{b^{*}}}{(a^{*})^{a^{*}}}\right)^{\frac{1}{b^{*}-a^{*}}}\right),\omega^{(k+1)}(x^{*})\right\\}.$
(82)

The values $a^{*}$ and $b^{*}$ are the maximizers of the function
$\omega^{(k+1)}\left(\frac{1}{e}\left(\frac{b^{b}}{a^{a}}\right)^{\frac{1}{b-a}}\right)$
given $a\neq b$ and $a,b>0$. The value $x^{*}$ is the maximizer of the
function $\omega^{(k+1)}(x)$ given $x>0$.

Proof:

This corrollary is proved by applying Theorem 1 from [20] to Eq. (65) and
Theorem 3. $\hfill\Box$

##  5 Continuity of Random Tensor Integral

In this section, we will establish continuity of DTI. We need the following
definition to define the convergence in mean for random tensors.

###### Definition 1

We say that a sequence of random tensor $\mathcal{X}_{n}$ converges in the
$r$-th mean towards the random tensor $\mathcal{X}$ with respect to the tensor
norm $\left\|\cdot\right\|_{\rho}$, if we have

$\displaystyle\mathbb{E}\left(\left\|\mathcal{X}_{n}\right\|^{r}_{\rho}\right)\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \mbox{exists,}$ (83)

and

$\displaystyle\mathbb{E}\left(\left\|\mathcal{X}\right\|^{r}_{\rho}\right)\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \mbox{exists,}$ (84)

and

$\displaystyle\lim\limits_{n\rightarrow\infty}\mathbb{E}\left(\left\|\mathcal{X}_{n}-\mathcal{X}\right\|^{r}_{\rho}\right)=0.$
(85)

We adopt the notatation $\mathcal{X}_{n}\xrightarrow[]{r}\mathcal{X}$ to
represent that random tensors $\mathcal{X}_{n}$ converges in the $r$-th mean
to the random tensor $\mathcal{X}$ with respect to the tensor norm
$\left\|\cdot\right\|_{\rho}$.

Besides random tensor convergence definition, we also need to define triple
tensor integral and second-order divide difference.

We define _triple tensor integrals_ (TTI) with respect to Hermitian tensors
$\mathcal{A},\mathcal{B},\mathcal{C}\in\mathbb{C}^{I_{1}\times\cdots\times
I_{N}\times I_{1}\times\cdots\times I_{N}}$ such that
$\mathcal{A}=\sum\limits_{i=1}^{\mathbb{I}_{1}^{N}}\lambda_{i}\mathcal{U}_{i}\star_{1}\mathcal{U}^{H}_{i}\stackrel{{\scriptstyle\mbox{\tiny
def}}}{{=}}\sum\limits_{i=1}^{\mathbb{I}_{1}^{N}}\lambda_{i}\mathcal{P}_{\mathcal{U}_{i}}$,
$\mathcal{B}=\sum\limits_{j=1}^{\mathbb{I}_{1}^{N}}\mu_{j}\mathcal{V}_{j}\star_{1}\mathcal{V}^{H}_{j}\stackrel{{\scriptstyle\mbox{\tiny
def}}}{{=}}\sum\limits_{j=1}^{\mathbb{I}_{1}^{N}}\mu_{j}\mathcal{P}_{\mathcal{V}_{j}}$
and
$\mathcal{C}=\sum\limits_{k=1}^{\mathbb{I}_{1}^{N}}\nu_{k}\mathcal{W}_{j}\star_{1}\mathcal{W}^{H}_{j}\stackrel{{\scriptstyle\mbox{\tiny
def}}}{{=}}\sum\limits_{k=1}^{\mathbb{I}_{1}^{N}}\nu_{k}\mathcal{P}_{\mathcal{W}_{k}}$.
Given the function $\varphi:\mathbb{R}^{3}\rightarrow\mathbb{R}$, the TTI
associated with tensors $\mathcal{A},\mathcal{B},\mathcal{C}$ and the function
$\varphi$, denoted as
$T_{\mathcal{A},\mathcal{B},\mathcal{C},\varphi}(\mathcal{X},\mathcal{Y})$,
can be expressed as

$\displaystyle
T_{\mathcal{A},\mathcal{B},\mathcal{C},\varphi}(\mathcal{X},\mathcal{Y})$
$\displaystyle=$
$\displaystyle\sum\limits_{i=1}^{\mathbb{I}_{1}^{N}}\sum\limits_{j=1}^{\mathbb{I}_{1}^{N}}\sum\limits_{k=1}^{\mathbb{I}_{1}^{N}}\varphi(\lambda_{i},\mu_{j},\nu_{k})\mathcal{P}_{\mathcal{U}_{i}}\star_{N}\mathcal{X}\star_{N}\mathcal{P}_{\mathcal{V}_{j}}\star_{N}\mathcal{Y}\star_{N}\mathcal{P}_{\mathcal{W}_{j}},$
(86)

where $\mathcal{X},\mathcal{Y}\in\mathbb{C}^{I_{1}\times\cdots\times
I_{N}\times I_{1}\times\cdots\times I_{N}}$.

The second-order divide difference for a function $f(x)$, denoted as
$f^{[2]}(x,y,z)$, can be defined as

$\displaystyle f^{[2]}(x,y,z)\stackrel{{\scriptstyle\mbox{\tiny
def}}}{{=}}\frac{f^{[1]}(y,z)-f^{[1]}(x,y)}{z-x}$ (87)

###### Lemma 5

Given three Hermitian tensors
$\mathcal{A},\mathcal{B},\mathcal{C}\in\mathbb{C}^{I_{1}\times\cdots\times
I_{N}\times I_{1}\times\cdots\times I_{N}}$ such that
$\mathcal{A}=\sum\limits_{i=1}^{\mathbb{I}_{1}^{N}}\lambda_{i}\mathcal{U}_{i}\star_{1}\mathcal{U}^{H}_{i}\stackrel{{\scriptstyle\mbox{\tiny
def}}}{{=}}\sum\limits_{i=1}^{\mathbb{I}_{1}^{N}}\lambda_{i}\mathcal{P}_{\mathcal{U}_{i}}$,
$\mathcal{B}=\sum\limits_{j=1}^{\mathbb{I}_{1}^{N}}\mu_{j}\mathcal{V}_{j}\star_{1}\mathcal{V}^{H}_{j}\stackrel{{\scriptstyle\mbox{\tiny
def}}}{{=}}\sum\limits_{j=1}^{\mathbb{I}_{1}^{N}}\mu_{j}\mathcal{P}_{\mathcal{V}_{j}}$
and
$\mathcal{C}=\sum\limits_{k=1}^{\mathbb{I}_{1}^{N}}\nu_{k}\mathcal{W}_{j}\star_{1}\mathcal{W}^{H}_{j}\stackrel{{\scriptstyle\mbox{\tiny
def}}}{{=}}\sum\limits_{k=1}^{\mathbb{I}_{1}^{N}}\nu_{k}\mathcal{P}_{\mathcal{W}_{k}}$
and the function $f(x)$ with $f^{[2]}(x)$ bounded by $\Omega_{f^{[2]}}$, we
then have the following norm estimate for
$T_{\mathcal{A},\mathcal{B},\mathcal{C},\varphi}(\mathcal{X},\mathcal{Y})$:

$\displaystyle\left\|T_{\mathcal{A},\mathcal{B},\mathcal{C},\varphi}(\mathcal{X},\mathcal{Y})\right\|_{\rho}$
$\displaystyle\leq$
$\displaystyle\left(\mathbb{I}_{1}^{N}\right)^{3}\Omega_{f^{[2]}}\left\|\mathcal{X}\right\|_{\rho}\cdot\left\|\mathcal{Y}\right\|_{\rho}.$
(88)

Proof:

Since we have

$\displaystyle\left\|T_{\mathcal{A},\mathcal{B},\mathcal{C},\varphi}(\mathcal{X},\mathcal{Y})\right\|_{\rho}$
$\displaystyle=$
$\displaystyle\left\|\sum\limits_{i=1}^{\mathbb{I}_{1}^{N}}\sum\limits_{j=1}^{\mathbb{I}_{1}^{N}}\sum\limits_{k=1}^{\mathbb{I}_{1}^{N}}\varphi(\lambda_{i},\mu_{j},\nu_{k})\mathcal{P}_{\mathcal{U}_{i}}\star_{N}\mathcal{X}\star_{N}\mathcal{P}_{\mathcal{V}_{j}}\star_{N}\mathcal{Y}\star_{N}\mathcal{P}_{\mathcal{W}_{j}}\right\|_{\rho}$
(89) $\displaystyle\leq_{1}$
$\displaystyle\sum\limits_{i=1}^{\mathbb{I}_{1}^{N}}\sum\limits_{j=1}^{\mathbb{I}_{1}^{N}}\sum\limits_{k=1}^{\mathbb{I}_{1}^{N}}\left|\varphi(\lambda_{i},\mu_{j},\nu_{k})\right|\left\|\mathcal{P}_{\mathcal{U}_{i}}\star_{N}\mathcal{X}\star_{N}\mathcal{P}_{\mathcal{V}_{j}}\star_{N}\mathcal{Y}\star_{N}\mathcal{P}_{\mathcal{W}_{j}}\right\|_{\rho}$
$\displaystyle\leq_{2}$
$\displaystyle\sum\limits_{i=1}^{\mathbb{I}_{1}^{N}}\sum\limits_{j=1}^{\mathbb{I}_{1}^{N}}\sum\limits_{k=1}^{\mathbb{I}_{1}^{N}}\left|\varphi(\lambda_{i},\mu_{j},\nu_{k})\right|\left\|\mathcal{X}\right\|_{\rho}\cdot\left\|\mathcal{Y}\right\|_{\rho}$
$\displaystyle\leq$
$\displaystyle\left(\mathbb{I}_{1}^{N}\right)^{3}\Omega_{f^{[2]}}\left\|\mathcal{X}\right\|_{\rho}\cdot\left\|\mathcal{Y}\right\|_{\rho}$

where $\leq_{1}$ comes from triangle inequality of the unitarily invariant
norm and $\leq_{2}$ comes from the definition of the unitarily invariant norm
and submultiplicative property of any unitarily invariant norm. $\hfill\Box$

Following theorem is about the continuity of a tensor integral.

###### Theorem 4

Let
$\mathcal{A}_{n},\mathcal{A},\mathcal{B}_{n},\mathcal{B}\in\mathbb{C}^{I_{1}\times\cdots\times
I_{N}\times I_{1}\times\cdots\times I_{N}}$ be random Hermitian tensors such
that

$\displaystyle\mathcal{A}_{n}\xrightarrow[]{r}\mathcal{A},\leavevmode\nobreak\
\mbox{and}\leavevmode\nobreak\ \mathcal{B}_{n}\xrightarrow[]{r}\mathcal{B},$
(90)

where $1\leq r<\infty$. Moreover, a real values function $f(x)$ for
$x\in\mathbb{R}$ such that $f^{[2]}(x)$ exists and bounded by
$\Omega_{f^{[2]}}$, respectively. Then, we have

$\displaystyle
T_{\mathcal{A}_{n},\mathcal{B}_{n},f^{[1]}}(\mathcal{X})\xrightarrow[]{r}T_{\mathcal{A},\mathcal{B},f^{[1]}}(\mathcal{X}),$
(91)

where $\mathcal{X}\in\mathbb{C}^{I_{1}\times\cdots\times I_{N}\times
I_{1}\times\cdots\times I_{N}}$ is a fixed tensor.

Proof:

From Lemma 3 and telescoping summation, we have the following:

$\displaystyle\left\|T_{\mathcal{A}_{n},\mathcal{B}_{n},f^{[1]}}(\mathcal{X})-T_{\mathcal{A},\mathcal{B},f^{[1]}}(\mathcal{X})\right\|_{\rho}$
$\displaystyle=$
$\displaystyle\left\|T_{\mathcal{A}_{n},\mathcal{B}_{n},f^{[1]}}(\mathcal{X})-T_{\mathcal{A},\mathcal{B}_{n},f^{[1]}}(\mathcal{X})\right.$
(92)
$\displaystyle\left.+T_{\mathcal{A},\mathcal{B}_{n},f^{[1]}}(\mathcal{X})-T_{\mathcal{A},\mathcal{B},f^{[1]}}(\mathcal{X})\right\|_{\rho}$
$\displaystyle=_{1}$
$\displaystyle\left\|T_{\mathcal{A}_{n},\mathcal{A},\mathcal{B}_{n},f^{[2]}}(\mathcal{A}_{n}-\mathcal{A},\mathcal{X})+T_{\mathcal{A},\mathcal{B}_{n},\mathcal{B},f^{[2]}}(X,\mathcal{B}_{n}-\mathcal{B})\right\|_{\rho}$
$\displaystyle\leq$
$\displaystyle\Omega_{f^{[2]}}\left\|\mathcal{A}_{n}-\mathcal{A}\right\|_{\rho}\left\|\mathcal{X}\right\|_{\rho}+\Omega_{f^{[2]}}\left\|\mathcal{B}_{n}-\mathcal{B}\right\|_{\rho}\left\|\mathcal{X}\right\|_{\rho},$

where $=_{1}$ is obtained from the $f^{[2]}$ definition given by Eq. (87), and
the $\leq$ comes from Lemma 5 and triangle inequality.

By raising the power $r$ and taking the expectation at the both sides of the
inequality provided by Eq. (92), we have proved this theorem by conditions
given by Eq. (90) and the following inequality:

$\displaystyle(a+b)^{r}\leq 2^{r}(a^{r}+b^{r})\leavevmode\nobreak\
\leavevmode\nobreak\ \mbox{given $a,b\geq 0$}.$ (93)

$\hfill\Box$

##  6 Applications of Tensor Integral

In this section, we will apply Theorem 4 and perturbation formulas provided by
Lemma 3 and Lemma 4 to bound the tail probability of the derivative of tensor-
valued function norm.

Given a fixed perturbation tensor
$\mathcal{X}\in\mathbb{C}^{I_{1}\times\cdots\times I_{N}\times
I_{1}\times\cdots\times I_{N}}$ with respect to the random Hermitian tensor
$\mathcal{A}\in\mathbb{C}^{I_{1}\times\cdots\times I_{N}\times
I_{1}\times\cdots\times I_{N}}$, and a tensor-valued function $f(x)$, we
define the derivative of $f(x)$ at $\mathcal{A}$ with respect to the
perturbation $\mathcal{X}$, represented by
$f^{\prime}_{\mathcal{X}}(\mathcal{A})$, as

$\displaystyle
f^{\prime}_{\mathcal{X}}(\mathcal{A})\stackrel{{\scriptstyle\mbox{\tiny
def}}}{{=}}\lim\limits_{t\rightarrow
0}\frac{f(\mathcal{A}+t\mathcal{X})-f(\mathcal{A})}{t}$ (94)

Following theorem is about the tail bound for the norm of
$f^{\prime}_{\mathcal{X}}(\mathcal{A})$.

###### Theorem 5

Given a fixed perturbation tensor
$\mathcal{X}\in\mathbb{C}^{I_{1}\times\cdots\times I_{N}\times
I_{1}\times\cdots\times I_{N}}$ and a random Hermitian tensor
$\mathcal{A}\in\mathbb{C}^{I_{1}\times\cdots\times I_{N}\times
I_{1}\times\cdots\times I_{N}}$ with a tensor-valued function $f(x)$. Suppose
we have
$\mathcal{A}=\sum\limits_{i=1}^{\mathbb{I}_{1}^{N}}\lambda_{i}\mathcal{P}_{\mathcal{U}_{i}}$.Then,
we have the tail bound for
$\left\|f^{\prime}_{\mathcal{X}}(\mathcal{A})\right\|_{\rho}$ as

$\displaystyle\mathrm{Pr}\left(\left\|f^{\prime}_{\mathcal{X}}(\mathcal{A})\right\|_{\rho}\geq\theta\right)$
$\displaystyle\leq$
$\displaystyle\frac{\left(\mathbb{I}_{1}^{N}\right)^{2}\left\|\mathcal{X}\right\|_{\rho}}{\theta}\sum\limits_{i=1}^{\mathbb{I}_{1}^{N}}\sum\limits_{j=1}^{\mathbb{I}_{1}^{N}}\mathbb{E}\left(\left|f^{[1]}(\lambda_{i},\lambda_{j})\right|\right).$
(95)

Proof:

From Lemma 3, we have

$\displaystyle f(\mathcal{A}+t\mathcal{X})-f(\mathcal{A})$ $\displaystyle=$
$\displaystyle
T_{\mathcal{A}+t\mathcal{X},\mathcal{A},f^{[1]}}\left(\mathcal{A}+t\mathcal{X}-\mathcal{A}\right).$
(96)

Then, we have

$\displaystyle\frac{f(\mathcal{A}+t\mathcal{X})-f(\mathcal{A})}{t}$
$\displaystyle=$ $\displaystyle
T_{\mathcal{A}+t\mathcal{X},\mathcal{A},f^{[1]}}\left(\mathcal{X}\right).$
(97)

From Theorem 4, we have

$\displaystyle f^{\prime}_{\mathcal{X}}(\mathcal{A})$ $\displaystyle=$
$\displaystyle\lim\limits_{t\rightarrow
0}\frac{f(\mathcal{A}+t\mathcal{X})-f(\mathcal{A})}{t}$ (98)
$\displaystyle\xrightarrow[]{r}$ $\displaystyle
T_{\mathcal{A},\mathcal{A},f^{[1]}}\left(\mathcal{X}\right).$

Since the convergence in the $r$-th mean implies the convergence in
probability, we then have

$\displaystyle\mathrm{Pr}\left(\left\|f^{\prime}_{\mathcal{X}}(\mathcal{A})\right\|_{\rho}\geq\theta\right)$
$\displaystyle=$
$\displaystyle\mathrm{Pr}\left(\left\|T_{\mathcal{A},\mathcal{A},f^{[1]}}\left(\mathcal{X}\right)\right\|_{\rho}\geq\theta\right)$
(99) $\displaystyle\leq$
$\displaystyle\frac{\left(\mathbb{I}_{1}^{N}\right)^{2}\left\|\mathcal{X}\right\|_{\rho}}{\theta}\sum\limits_{i=1}^{\mathbb{I}_{1}^{N}}\sum\limits_{j=1}^{\mathbb{I}_{1}^{N}}\mathbb{E}\left(\left|f^{[1]}(\lambda_{i},\lambda_{j})\right|\right),$

where the inequality comes from Theorem 1 $\hfill\Box$

The commutator of a tensor $\mathcal{A}\in\mathbb{C}^{I_{1}\times\cdots\times
I_{N}\times I_{1}\times\cdots\times I_{N}}$ with respect to the tensor
$\mathcal{B}\in\mathbb{C}^{I_{1}\times\cdots\times I_{N}\times
I_{1}\times\cdots\times I_{N}}$ if we have
$\mathcal{A}\star_{N}\mathcal{B}=\mathcal{B}\star_{N}\mathcal{A}$.

Given a fixed perturbation tensor
$\mathcal{X}\in\mathbb{C}^{I_{1}\times\cdots\times I_{N}\times
I_{1}\times\cdots\times I_{N}}$ with respect to the random Hermitian tensor
$\mathcal{A}\in\mathbb{C}^{I_{1}\times\cdots\times I_{N}\times
I_{1}\times\cdots\times I_{N}}$, and a tensor-valued function $f(x)$, we
define the derivative of $f(x)$ at $\mathcal{A}$ with respect to the
perturbation $\mathcal{X}$ and the commutator tensor
$\mathcal{D}\in\mathbb{C}^{I_{1}\times\cdots\times I_{N}\times
I_{1}\times\cdots\times I_{N}}$ of $\mathcal{A}$, represented by
$f^{\prime}_{\mathcal{X}|\mathcal{D}}(\mathcal{A})$, as

$\displaystyle
f^{\prime}_{\mathcal{X}|\mathcal{D}}(\mathcal{A})\stackrel{{\scriptstyle\mbox{\tiny
def}}}{{=}}\lim\limits_{t\rightarrow
0}\frac{\mathcal{D}\star_{N}f(\mathcal{A}+t\mathcal{X})-f(\mathcal{A})\star_{N}\mathcal{D}}{t}$
(100)

Following theorem is about the tail bound for the norm of
$f^{\prime}_{\mathcal{X}|\mathcal{D}}(\mathcal{A})$.

###### Theorem 6

Given a fixed perturbation tensor
$\mathcal{X}\in\mathbb{C}^{I_{1}\times\cdots\times I_{N}\times
I_{1}\times\cdots\times I_{N}}$ and a random Hermitian tensor
$\mathcal{A}\in\mathbb{C}^{I_{1}\times\cdots\times I_{N}\times
I_{1}\times\cdots\times I_{N}}$ with a tensor-valued function $f(x)$. Suppose
we have
$\mathcal{A}=\sum\limits_{i=1}^{\mathbb{I}_{1}^{N}}\lambda_{i}\mathcal{P}_{\mathcal{U}_{i}}$
and the tensor $\mathcal{D}\in\mathbb{C}^{I_{1}\times\cdots\times I_{N}\times
I_{1}\times\cdots\times I_{N}}$ is the commutator of the tensor $\mathcal{A}$.
Then, we have the tail bound for
$\left\|f^{\prime}_{\mathcal{X}|\mathcal{D}}(\mathcal{A})\right\|_{\rho}$ as

$\displaystyle\mathrm{Pr}\left(\left\|f^{\prime}_{\mathcal{X}|\mathcal{D}}(\mathcal{A})\right\|_{\rho}\geq\theta\right)$
$\displaystyle\leq$
$\displaystyle\frac{\left(\mathbb{I}_{1}^{N}\right)^{2}\left\|\mathcal{D}\star_{N}\mathcal{X}\right\|_{\rho}}{\theta}\sum\limits_{i=1}^{\mathbb{I}_{1}^{N}}\sum\limits_{j=1}^{\mathbb{I}_{1}^{N}}\mathbb{E}\left(\left|f^{[1]}(\lambda_{i},\lambda_{j})\right|\right).$
(101)

Proof:

From Lemma 4, we have

$\displaystyle\mathcal{D}\star_{N}f(\mathcal{A}+t\mathcal{X})-f(\mathcal{A})\star_{N}\mathcal{D}$
$\displaystyle=$ $\displaystyle
T_{\mathcal{A}+t\mathcal{X},\mathcal{A},f^{[1]}}\left(\mathcal{D}\star_{N}\left(\mathcal{A}+t\mathcal{X}\right)-\mathcal{A}\star_{N}\mathcal{D}\right).$
(102)

Also, from $\mathcal{D}\star_{N}\mathcal{A}=\mathcal{A}\star_{N}\mathcal{D}$,
we also have

$\displaystyle\frac{\mathcal{D}\star_{N}f(\mathcal{A}+t\mathcal{X})-f(\mathcal{A})\star_{N}\mathcal{D}}{t}$
$\displaystyle=$ $\displaystyle
T_{\mathcal{A}+t\mathcal{X},\mathcal{A},f^{[1]}}\left(\mathcal{D}\star_{N}\mathcal{X}\right).$
(103)

From Theorem 4, we have

$\displaystyle f^{\prime}_{\mathcal{X}|\mathcal{D}}(\mathcal{A})$
$\displaystyle=$ $\displaystyle\lim\limits_{t\rightarrow
0}\frac{\mathcal{D}\star_{N}f(\mathcal{A}+t\mathcal{X})-f(\mathcal{A})\star_{N}\mathcal{D}}{t}$
(104) $\displaystyle\xrightarrow[]{r}$ $\displaystyle
T_{\mathcal{A},\mathcal{A},f^{[1]}}\left(\mathcal{D}\star_{N}\mathcal{X}\right).$

Since the convergence in the $r$-th mean implies the convergence in
probability, we then have

$\displaystyle\mathrm{Pr}\left(\left\|f^{\prime}_{\mathcal{X}|\mathcal{D}}(\mathcal{A})\right\|_{\rho}\geq\theta\right)$
$\displaystyle=$
$\displaystyle\mathrm{Pr}\left(\left\|T_{\mathcal{A},\mathcal{A},f^{[1]}}\left(\mathcal{D}\star_{N}\mathcal{X}\right)\right\|_{\rho}\geq\theta\right)$
(105) $\displaystyle\leq$
$\displaystyle\frac{\left(\mathbb{I}_{1}^{N}\right)^{2}\left\|\mathcal{D}\star_{N}\mathcal{X}\right\|_{\rho}}{\theta}\sum\limits_{i=1}^{\mathbb{I}_{1}^{N}}\sum\limits_{j=1}^{\mathbb{I}_{1}^{N}}\mathbb{E}\left(\left|f^{[1]}(\lambda_{i},\lambda_{j})\right|\right),$

where the inequality comes from Theorem 1 again. $\hfill\Box$

##  7 Conclusions

We first define what is the random DTI and derive the tail bound of the
unitarily invariant norm for a random DTI. This bound assists us to establish
tail bounds of the unitarily invariant norm for various types of dual tensor
means, e.g., arithmetic mean, geometric mean, harmonic mean, and general mean.
The random DTI is also being applied to build the random Lipschitz estimate in
contexts of random tensors. Finally, we derive the continuity property for
random DTI in the sense of convergence in the random tensor mean, and apply
this fact to obtain the tail bound of the unitarily invariant norm for the
derivative of the tensor-valued function. Possible future works will be to
extend DTI to multiple tensor integrals.

## References

  * [1] L. Qi and Z. Luo, _Tensor analysis: spectral theory and special tensors_. SIAM, 2017.
  * [2] N. D. Sidiropoulos, L. De Lathauwer, X. Fu, K. Huang, E. E. Papalexakis, and C. Faloutsos, “Tensor decomposition for signal processing and machine learning,” _IEEE Transactions on Signal Processing_ , vol. 65, no. 13, pp. 3551–3582, Jul. 2017.
  * [3] R. Gurau, _Random tensors_. Oxford University Press, Oxford, 2017.
  * [4] I. R. Klebanov and G. Tarnopolsky, “Uncolored random tensors, melon diagrams, and the Sachdev-Ye-Kitaev models,” _Phys. Rev. D_ , vol. 95, no. 4, pp. 046 004, 13, 2017. [Online]. Available: https://doi.org/10.1103/physrevd.95.046004
  * [5] R. Vershynin, “Concentration inequalities for random tensors,” _Bernoulli_ , vol. 26, no. 4, pp. 3139–3162, 2020. [Online]. Available: https://doi.org/10.3150/20-BEJ1218
  * [6] S. Y. Chang and W.-W. Lin, “Convenient tail bounds for sums of random tensors,” _Taiwanese Journal of Mathematics_ , vol. 1, no. 1, pp. 1–36, 2021\.
  * [7] S. Y. Chang, “General tail bounds for random tensors summation: Majorization approach,” 2021.
  * [8] ——, “Hanson-wright inequality for random tensors under einstein product,” _arXiv preprint arXiv:2111.12169_ , 2021.
  * [9] ——, “Generalized hanson-wright inequality for random tensors,” _arXiv preprint arXiv:2203.00659_ , 2022.
  * [10] M. E. Kilmer, C. D. Martin, and L. Perrone, “A third-order generalization of the matrix svd as a product of third-order tensors,” _Tufts University, Department of Computer Science, Tech. Rep. TR-2008-4_ , 2008.
  * [11] S. Y. Chang and Y. Wei, “T-square tensors—Part I: inequalities,” _Computational and Applied Mathematics_ , vol. 41. [Online]. Available: https://doi.org/10.1007/s40314-022-01770-0
  * [12] ——, “T-square tensors—Part II: tail bounds for sums of random t-product tensors,” _Computational and Applied Mathematics_ , vol. 41. [Online]. Available: https://doi.org/10.1007/s40314-022-01811-8
  * [13] F. Hiai, “Matrix analysis: matrix monotone functions, matrix means, and majorization,” _Interdisciplinary Information Sciences_ , vol. 16, no. 2, pp. 139–248, 2010.
  * [14] A. Skripka and A. Tomskova, _Multilinear Operator Integrals Theory and Applications Introduction_. Springer, 2019, vol. 2250.
  * [15] M. Liang and B. Zheng, “Further results on Moore-Penrose inverses of tensors with application to tensor nearness problems,” _Comput. Math. Appl._ , vol. 77, no. 5, pp. 1282–1293, 2019. [Online]. Available: https://doi.org/10.1016/j.camwa.2018.11.001
  * [16] ——, “Further results on moore–penrose inverses of tensors with application to tensor nearness problems,” _Computers & Mathematics with Applications_, vol. 77, no. 5, pp. 1282–1293, Mar. 2019.
  * [17] D. Huang, “Generalizing lieb’s concavity theorem via operator interpolation,” _Advances in Mathematics_ , vol. 369, p. 107208, 2020.
  * [18] K. Fan and A. J. Hoffman, “Some metric inequalities in the space of matrices,” _Proceedings of the American Mathematical Society_ , vol. 6, no. 1, pp. 111–116, 1955.
  * [19] L. G. Ordóñez, D. P. Palomar, and J. R. Fonollosa, “Ordered eigenvalues of a general class of hermitian random matrices with application to the performance analysis of MIMO systems,” _IEEE Transactions on Signal Processing_ , vol. 57, no. 2, pp. 672–689, 2008.
  * [20] F. Qi, S. Guo, and S.-X. Chen, “A new upper bound in the second kershaw’s double inequality and its generalizations,” _Journal of computational and applied mathematics_ , vol. 220, no. 1-2, pp. 111–118, 2008.


# Double occupancies in confined attractive fermions on optical lattices

Ji-Hong Hu  Department of Physics, Zhejiang Normal University, Jinhua 340012,
China  Gao Xianlong<EMAIL_ADDRESS>Department of Physics, Zhejiang Normal
University, Jinhua 340012, China

###### Abstract

We perform a numerical study of a one-dimensional Fermion-Hubbard model in
harmonic traps within the Thomas-Fermi approximation based on the exact Bethe-
ansatz solution. The $\rho-U/t$ phase diagram is shown for the systems of
attractive interactions ($\rho$ is the characteristic density and $U/t$ the
interaction strength scaled in units of the hopping parameter.). We study the
double occupancy, the local central density and their derivatives. Their roles
are discussed in details in detecting the composite phases induced by the
trapping potential.

###### pacs:

05.30.Fk,03.75.Ss,71.10.Pm,71.15.Pd

The recent remarkable experimental progress in cooling and manipulating
ultracold atomic gases in optical lattices provides us an alternative way to
investigate the many-body effects in condensed matter system review paper .
For the optical lattices loaded with atomic gases, the on-site interaction
between different species can be tuned to be very strong by means of a
technique named Feshbach resonance Feshbach  or by increasing the lattice
depth to decrease the hopping between neighbor sites. Thus the outstanding
controllability offers the possibility to use these systems to ”quantum
simulate” quantum many-body physics in a well-controlled way. An example is
the experimentally realized Tonks-Girardeau gas in a one-dimensional (1D)
optical lattice Paredes .

In the experiments, the trapping potential used to confine atoms induces
inhomogeneity and complexity in characterizing the quantum phases. Though a
few techniques exist in experiments observing the new quantum phases and
probing various properties in ultracold gases in optical lattices, such as the
time-of-flight imaging of the momentum distribution, noise correlations and
Bragg spectroscopy review paper , there are difficulties in observing the
composite phases induced by the inherent trapping potential, for example, the
incompressible bulk Mott- and band-insulating phases surrounded by the
compressible fluid phases at the edges in the 1D Fermi-Hubbard model under
harmonic confinement. Scarola et al. proposed to use the core compressibility
by filtering out the edge effects, which offers a direct probe of
incompressible phases independent of inhomogeneity Scarola . Duan et al.
suggested to use the Fourier sampling of time-of-flight images to reveal new
correlation functions Duan . Kollath et al. used the double occupancy (DO)
induced by periodic lattice modulations via the creation of molecules to
identify the pair-binding energy and the spin-ordering in a Fermi gas in an
optical lattice Kollath , which was recently measured in the experiments by
molecule formation and used to identify an incompressible Mott-insulator phase
of two hyperfine states of strongly repulsive fermionic 40K atoms Jordens .
The compressibility of the quantum gas in the trap was used recently in
identifying the incompressible Mott-insulating phase at strong interactions by
independent control of external confinement and lattice depth Schneider . Here
in this paper we discuss how to characterize the different composite phases of
1D confined fermions on optical lattices, by investigating observables like
the double occupation sites, the local central density and their derivatives.

We consider a two-component Fermi gas in a tube with $N_{f}$ atoms and $N_{s}$
lattice sites, which can be realized by a strong confinement in the transverse
directions moritz_2005  with an additional periodic potential applied along
the tube and can be described by a one-band inhomogeneous Fermi-Hubbard model
jaksch_98 ,

$\displaystyle\hat{H}_{s}$ $\displaystyle=$
$\displaystyle-t\sum_{i=1,\sigma}^{N_{s}}({\hat{c}}^{\dagger}_{i\sigma}{\hat{c}}_{i+1\sigma}+{H}.{c}.)+U\sum_{i=1}^{N_{s}}\,{\hat{D}}_{i}$
(1) $\displaystyle+V_{2}\sum_{i=1}^{N_{s}}(i-N_{s}/2)^{2}{\hat{n}}_{i}$

where the spin degrees of freedom $\sigma=\uparrow,\downarrow$ are
pseudospin-$1/2$ labels for two internal hyperfine states.
${\hat{c}}_{i\sigma}$ (${\hat{c}}^{\dagger}_{i\sigma}$) are fermionic
operators annihilating (creating) particles with spin $\sigma$ in a Wannier
state at site $i$.
${\hat{n}}_{i}=\sum_{\sigma}{\hat{n}}_{i\sigma}=\sum_{\sigma}{\hat{c}}^{\dagger}_{i\sigma}{\hat{c}}_{i\sigma}$
is the total site occupation operator, $t$ is the tunneling between nearest
neighbors, $U$ is the strength of the on-site attractive interaction,
${\hat{D}}_{i}={\hat{n}}_{i\uparrow}{\hat{n}}_{i\downarrow}$ is the double
occupancy operator at site $i$, and $V_{i}$ describes the strength of the
trapping potential.

The homogeneous 1D Fermi-Hubbard model belongs to the universality class of
Luttinger liquids. At zero temperature, the properties of this model in the
thermodynamic limit ($N_{\sigma},N_{s}\rightarrow\infty$, but with finite
$N_{\sigma}/N_{s}$) are determined by the fillings
$n_{\sigma}=N_{\sigma}/N_{s}$ and by the dimensionless coupling constant
$u=U/t$. According to Lieb and Wu LiebWu , the ground state properties for
different fillings in the thermodynamic limit are described by the coupled
integral equations (details see Refs. KocharianPRB, ; gaoprb78, ).

Figure 1: (Color online) Phase diagram for the one-dimensional inhomogeneous
fermion-Hubbard model as a function of attractive interaction strength $u=U/t$
and characteristic particle density $\rho=N_{f}\sqrt{V_{2}/t}$. For fixed $u$,
the shape of the density profile is completely determined by $\rho$. The band
insulator (BI) phase is determined with the local occupation in the trap
center satisfying $|n_{i}-2|<0.005$. Below the BI phase, it is the Luther-
Emery (LE) liquid phase. The boundary between BI and LE phases is plotted as
open squares ($\square$) with solid lines as a guide for the eye. The solid
circles ($\bullet$) represent the phase boundary determined from the
derivative of the double occupancy, and the solid upper triangles
($\blacktriangle$) represent the one determined from the derivative of the
local central density (See details in the text and in Fig. 4). All of them can
be used as identifying the onset of the BI phase.

Figure 2: (Color online) The ground-state concentration of doubly occupied
sites in the homogeneous system $D^{\rm hom}$ (a) as a function of interaction
strength $u$ for various $n$ and (b) as a function of particle density $n$ for
various interacting strengths $u$. For $0\leq n\leq 2$, $D^{\rm hom}$ is in
the range $0\leq D^{\rm hom}\leq 1$. For fixed $n$, $D^{\rm hom}$ increases
monotonically with increasing $|u|$ and for fixed $u$ it increases
monotonically with increasing $n$. At $U=0$, $D^{\rm hom}=n^{2}/4$ and at
strongly attractive coupling limit $D^{\rm hom}=n/2$. The attractive
interaction favors the formation of pairs between opposite spins. In the inset
of (b), the derivative of the double occupancy with respect to the particle
density $n$ is plotted, which is used to calculate the double occupancy for
the inhomogeneous system in the spirit of local density approximation.

The inhomogeneity of 1D Fermi-Hubbard model caused by the boundary conditions
or the external potential, normally invalidates a reliable analytical method
that usually used in the homogeneous system. Many numerically accurate
schemes, such as the density-matrix renormalization group Machida1 ; gaoprb78
; Molina , quantum-Monto-carlo (QMC) Rigol ; Rigol2004 ; Astrakhardik ; Casula
, exact diagnolization Machida1 ; exact diagnolization ; Nikkarila  and
density-functional theory based on the exact Bethe-ansatz solution gaoprb78 ,
are used in obtaining the phase diagram and the collective oscillations of the
atomic mass density. The Thomas-Fermi approximation based on the exact Bethe-
ansatz solution has also been successfully used in characterizing the phase
diagrams Liu . The rich phases induced by the trapping potential or the phase
transitions driven by the system parameters can be identified by the site
occupation, the variance of the local density, the local compressibility
Rigol2004 , or the thermodynamic stiffness (i.e., the inverse of the global
compressibility) gaoprb78 .

Here we study the inhomogeneous Fermi-Hubbard model described by Eq. (1),
which allows for the spatial coexistence of different quantum states Rigol ;
Liu . For attractive interactions, $U<0$, there are two different composite
phases, well classified in the phase diagram defined by a scaled dimensionless
characteristic density $\rho=N_{f}\sqrt{V_{2}/t}$ and $u=U/t$ Rigol : the
Luther-Emery (LE) liquid phase characterized by fermion pairs and atomic
density wave, and band insulator (BI) phase with a plateau of $n_{i}=2$ in the
trap center of tightly bound spin-singlet dimers surrounded by Luther-Emery
layers gaoprb78 .

All these phases have been analyzed in the $\rho-u$ Rigol2004 ; Heiselberg ,
$\mu-u$ ($\mu$ is the chemical potential) Campo , and $n_{0}-\nu$ Liu
($n_{0}$ is the central density and $\nu=2\sqrt{N_{f}}/\pi$ is the
characteristic filling factor) phase diagrams.

Figure 3: (Color online) (a) The ground-state concentration of doubly occupied
sites $D$ as a function of characteristic density $\rho$ in the system under
the presence of trapping potential for various interacting strength $u$. (b)
The local central density $n_{0}$ as a function of characteristic density
$\rho$ for various interacting strength $u$. For comparison, the double
occupancy and the local central density for the noninteracting case ($u=0$)
are also shown in (a) and (b), respectively.

Figure 4: (Color online) (a) The derivative of the ground-state concentration
of doubly occupied sites $D$ as a function of characteristic density $\rho$ in
the system under the presence of trapping potential. In the inset, the local
inverse compressibility $\kappa^{-1}_{i}$ (in units of $t$) is shown as a
function of $i$ with the system parameters determined from the critical
characteristic density on the peaks. (b) The derivative of the local central
density $n_{0}$ as a function of characteristic density $\rho$ for different
attractive interactions. For comparison, the derivatives of the double
occupancy and the local central density for the noninteracting case ($u=0$)
are also shown in (a) and (b), respectively.

The harmonic trapping potential in Eq. (1) induces a position dependence of
the chemical potentials which can be evaluated within the local density
approximation based on the Bethe-ansatz equations. This amounts to determining
the global chemical potential $\mu$ of the system from the local equilibrium
condition,

$\mu=\left.\mu_{\rm hom}(n,u)\right|_{n\rightarrow
n_{i}}+V_{2}(i-N_{s}/2)^{2},$ (2)

while imposing the normalization condition $\sum_{i}n_{i}=N_{f}$. Here
$\mu_{\rm hom}(n,u)$ is the chemical potential of the homogeneous system
calculated from the Bethe-ansatz coupled equations gao_long_2005 ; gaoprb78 .

In Fig. 1, we show the above-mentioned composite phases through the $\rho-u$
phase diagram, which is insensitive to the individual values of $N_{f}$ and
$V_{2}$, that is, for fixed $u$, the shape of the density profile in this
system is completely determined by $\rho$ Rigol . Below the BI phase, it is
the LE liquid phase. We would like to mention that the similar phase diagram
has been obtained by Rigol et al. for the repulsive interaction Rigol  and
Heiselberg but with a fitting function for $\mu_{\rm hom}(n,u)$ Heiselberg .
Here we solve the local equilibrium condition with the exact numerical data
for $\mu_{\rm hom}(n,u)$ calculated from the coupled Bethe-ansatz equations.
We conclude that the phase diagram based on the Thomas-Fermi approximation
gives qualitatively correct results when comparing to QMC.

To detect the composite phases induced by the trapping potential, we
investigate firstly the concentration of the doubly occupied sites, which is
defined as the derivative of the ground state energy with respect to $U$ Leo ,

$D=\frac{2}{N_{f}}\sum^{N_{s}}_{i=1}\langle\hat{D}_{i}\rangle=\frac{2}{N_{f}}\frac{\partial
E_{0}}{\partial U}\,.$

Here $E_{0}$ is the ground state energy per site of the system. The DO in the
homogeneous system, $D^{\rm hom}$, can be calculated exactly from the Bethe-
ansatz coupled equations KocharianPRB , which is shown in Fig. 2. Notice, at
$U=0$, we have $D^{\rm hom}=n^{2}/4$ and at strongly attractive coupling limit
$D^{\rm hom}=n/2$. In the presence of trapping potential, we calculate the DO
in the spirit of local density approximation, $D=\frac{2}{N_{f}}\sum_{i}D^{\rm
hom}(n)|_{n\rightarrow n_{i}}$. From Fig. 2, it becomes obvious that the
attractive interaction favors the formation of pairs between the two spin
species while the repulsive interaction suppresses them.

Now we show in Fig. 3 (a) the DO and in Fig. 3 (b) the local central density
$n_{0}$ for the inhomogeneous system described in Eq. (1) as a function of the
characteristic density $\rho$. For large $\rho$, the DO and the central
density saturate at $D\approx 1$ and $n_{0}\approx 2$, respectively,
corresponding to the BI regime, where an incompressible phase of band-
insulating type is realized in the core of the trap. The more attractive of
the interaction is, the sooner for $D$ to reach the value of $1$ and for
$n_{0}$ to reach $2$. These manifest that the attractive interaction favors
pairing, and thus increases DO, different from the case of repulsive
interaction where it suppresses DO and supports Mott-insulating phase.

In the following, we show the derivative of the DO in Fig. 4 (a) and the
derivative of the local central density in Fig. 4 (b), respectively, as probes
for the composite quantum phase transitions and especially as signs for band-
insulating phases. The derivatives of the DO and the local central density are
performed numerically using the central finite difference method. We find that
the singularities in these two quantities provide a further support and clear
signatures of the composite quantum phase transitions accompanying Fig. 3. In
Fig. 4 (a), the phase transition from LE to BI is manifested by an obvious
peak. The onset of the peak gives a critical value of the characteristic
density, which determines the critical particle number for a given $V_{2}/t$.
Based on these parameters, we plot in the inset of Fig. 4 (a) the local
inverse compressibility defined by $\kappa^{-1}_{i}=\partial^{2}E_{0}/\partial
n^{2}|_{n\rightarrow n_{i}}$ as a function of site position. The
incompressible nature of the BI phase is manifested by the disappearance of
$\kappa^{-1}_{i}$, which is compatible with the onset of the peak in
$dD/d\rho$. In Fig. 4 (b), the vanishing of the derivative of the local
central density gives very direct signatures of the incompressible nature of
BI in the center of the trap. The phase boundaries of the onset of the BI
phase are plotted in Fig. 1, which are determined by the appearance of the
peak in the derivative of the DO and the vanishing of the derivative of the
local central density. We find that the phase boundary determined from the DO
slightly overestimates the onset of the BI phase, while the central density
slightly underestimates it.

In conclusion, motivated by earlier theoretical work and by ongoing
experimental efforts, we present a fully numerical study for the ground state
properties of a one-dimensional attractive Fermion-Hubbard model in harmonic
traps within the Thomas-Fermi approximation based on the exact Bethe-ansatz
solution. We investigate the double occupancy, the local central density and
their derivatives. We found they are good indicators of the formation of an
incompressible band-insulating phase in the center of the trap, which can be
detected in the near-future in the experiments in the 1D cold fermionic atomic
gases.

This work was supported by NSF of China under Grant No. 10704066, 10974181,
Qianjiang River Fellow Fund 2008R10029, and Program for Innovative Research
Team in Zhejiang Normal University.

## References

  * (1) M. Lewenstein, A. Sanpera, V. Ahufinger, B. Damski, A. Sen (De), and U. Sen, Adv. Phys. 56, 243 (2007); I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 80, 885 (2008). S. Giorgini, L. P. Pitaevskii, and S. Stringari, Rev. Mod. Phys. 80, 1215 (2008).
  * (2) H. Feshbach, Ann. Phys. 5, 357 (1958); C. A. Regal, M. Greiner, and D. S. Jin, Phys. Rev. Lett. 92, 040403 (2004).
  * (3) B. Paredes, A. Widera, V. Murg, O. Mandel, S. Fölling, I. Cirac, G. V. Shlyapnikov, T. W. Hänsch, and I. Bloch, Nature 429, 277 (2004); T. Kinoshita, T. Wenger, and D. S. Weiss, Science 305, 1125 (2004).
  * (4) V. W. Scarola, L. Pollet, J. Oitmaa, and M. Troyer, Phys. Rev. Lett. 102, 135302 (2009).
  * (5) B. Wang and L. M. Duan, New J. Phys. 10, 073007 (2008).
  * (6) C. Kollath, A. Iucci, I. P. McCulloch, and T. Giamarchi, Phys. Rev. A 74, 041604(R) (2006).
  * (7) R. Jördens, N. Strohmaier, K. Günter, H. Moritz, and T. Esslinger, Nature (London) 455, 204 (2008). U. Schneider, L. Hackermüller, S. Will, Th. Best, I. Bloch, T. A. Costi, R. W. Helmes, D. Rasch, and A. Rosch, Science 322, 1520 (2008).
  * (8) U. Schneider, L. Hackermuller, S. Will, Th. Best, I. Bloch, T. A. Costi, R. W. Helmes, D. Rasch, and A. Rosch, Science 322, 1520 (2008).
  * (9) H. Moritz, T. Stöferle, K. Günter, M. Köhl, and T. Esslinger, Phys. Rev. Lett. 94, 210401 (2005).
  * (10) D. Jaksch, C. Bruder, J. I. Cirac, C. W. Gardiner, and P. Zoller, Phys. Rev. Lett. 81, 3108 (1998); W. Hofstetter, J. I. Cirac, P. Zoller, E. Demler, and M. D. Lukin, Phys. Rev. Lett. 89, 220407 (2002).
  * (11) E. H. Lieb and F. Y. Wu, Phys. Rev. Lett 20, 1445 (1968).
  * (12) G. Xianlong, M. Rizzi, M. Polini, R. Fazio, M. P. Tosi, V. L. Campo, Jr., and K. Capelle, Phys. Rev. Lett 98, 030404 (2007); G. Xianlong, Phys. Rev. B78, 085108 (2008).
  * (13) A. N. Kocharian, C. Yang, and Y. L. Chiang, Phys. Rev. B 59, 7458 (1999).
  * (14) M. Machida, S. Yamada, Y. Ohashi, and H. Matsumoto, Phys. Rev. A 74, 053621 (2006); M. Machida, S. Yamada, M. Okumura, Y. Ohashi, and H. Matsumoto, Phys. Rev. A 77, 053614 (2008).
  * (15) R. A. Molina, J. Dukelsky, and P. Schmitteckert, Phys. Rev. Lett. 99, 080404 (2007).
  * (16) M. Casula, D. M. Ceperley, and E. J. Mueller, Phys. Rev. A 78, 033607 (2008).
  * (17) G. E. Astrakharchik, D. Blume, S. Giorgini, and L. P. Pitaevskii, Phys. Rev. Lett. 93, 050402 (2004).
  * (18) M. Rigol, A. Muramatsu, G. G. Batrouni, and R. T. Scalettar, Phys. Rev. Lett. 91, 130403 (2003).
  * (19) M. Rigol and A. Muramatsu, Phys. Rev. A 69, 053612 (2004); Opt. Commun. 243, 33 (2004).
  * (20) T. Husslein, W. Fettes, and I. Morgenstern, Int. J. Mod. Phys. C 8, 397 (1997).
  * (21) J.-P. Nikkarila, M. Koskinen, S. M. Reimann, and M. Manninen, New J. Phys. 10, 063013 (2008); J.-P. Nikkarila, M. Koskinen, and M. Manninen, Eur. Phys. J. B 64, 95 (2008).
  * (22) X.-J. Liu, P. D. Drummond, and H. Hu, Phys. Rev. Lett. 94, 136406 (2005).
  * (23) H. Heiselberg, Phys. Rev. A 74, 033608 (2006).
  * (24) V. L. Campo, Jr. and K. Capelle, Phys. Rev. A 72, 061602(R) (2005).
  * (25) G. Xianlong, M. Polini, M.P. Tosi, V.L. Campo, Jr., K. Capelle, and M. Rigol, Phys. Rev. B 73, 165120 (2006).
  * (26) L. De Leo, C. Kollath, A. Georges, M. Ferrero, and O. Parcollet, Phys. Rev. Lett. 101, 210403 (2008).