# Extremes of generalized inversions and descents
on permutation groups

Philip Dörr  Department of Mathematics, Magdeburg…<EMAIL_ADDRESS>
###### Abstract.

The classical inversion statistic on symmetric groups is the sum of all
indicators $\textbf{1}\\{\pi(i)>\pi(j)\\}$ for a random permutation
$\pi=(\pi(1),\ldots,\pi(n))$ and the pairs $(i,j)$ with $1\leq i<j\leq n$. The
descent statistic counts all $i\in\\{1,\ldots,n-1\\}$ with $\pi(i)>\pi(i+1)$.
The number of inversions can be generalized by restricting the indicators to
pairs $(i,j)$ with $i<j$ and $|i-j|\leq d$ for some $d\in\\{1,\ldots,n-1\\}$.
Likewise, the number of descents can be generalized by counting all
$i\in\\{1,\ldots,n-d\\}$ with $\pi(i)>\pi(i+d)$. These generalized statistics
can be further extended to the signed and even-signed permutation groups. The
bandwidth index $d$ can be chosen in dependence of $n$, and the magnitude of
$d$ is significant for asymptotic considerations. It is known that each of
these statistics is asymptotically normal for suitable choices of $d$. In this
paper we prove the bivariate asymptotic normality and determine the extreme
value asymptotics of generalized inversions and descents.

###### Key words and phrases:

extreme values, CLT, permutation statistics, Coxeter group

###### 2010 Mathematics Subject Classification:

Primary: 60G70, 05A16; Secondary: 20F55

Philip Dörr is supported by the DFG (314838170, GRK 2297, “MathCoRe”).

##  1\. Introduction

The numbers of inversions and descents are two important quantities of
permutations. For a permutation
$\pi\negmedspace:\\{1,\ldots,n\\}\rightarrow\\{1,\ldots,n\\},$ an inversion is
a pair $(i,j)$ with $i<j,$ but $\pi(i)>\pi(j)$. A descent is an index $i$ with
$\pi(i)>\pi(i+1)$, i.e., descents correspond to adjacent inversions. We now
equip the symmetric group $S_{n}$ with the discrete uniform probability
measure induced by the point masses $\mathbb{P}(\\{\pi\\})=1/n!$
$\forall\pi\in S_{n}$. In this context, we denote the random numbers of
inversions and descents as $X_{\mathrm{inv}}$ and $X_{\mathrm{des}}$. These
random variables can also be represented with help of i.i.d. standard uniform
variables $Z_{1},Z_{2},\ldots,\sim U(0,1)$, namely,

(1) $\displaystyle X_{\mathrm{inv}}$ $\displaystyle=\sum_{1\leq i<j\leq
n}\textbf{1}\\{Z_{i}>Z_{j}\\}\,,$ $\displaystyle X_{\mathrm{des}}$
$\displaystyle=\sum_{i=1}^{n-1}\textbf{1}\\{Z_{i}>Z_{i+1}\\}\,.$

In (1), we have to sum the indicators $\textbf{1}\\{Z_{i}>Z_{j}\\}$ over all
pairs $(i,j)$ with $i<j$. A class of generalized inversion statistics can be
constructed by restricting the sum in (1) to pairs $(i,j)$ with $1\leq j-i\leq
d$, for some $d\in\\{1,\ldots,n-1\\}$. Writing
$\mathfrak{N}_{n,d}:=\\{(i,j)\in\\{1,\ldots,n\\}^{2}\mid 1<j-i\leq d\\}$, we
can state that

(2)
$X_{\mathrm{inv}}^{(d)}:=\sum_{(i,j)\in\mathfrak{N}_{n,d}}\textbf{1}\\{Z_{i}>Z_{j}\\}$

counts the so-called $d$-inversions. The term generalized inversions is an
umbrella term for all $d$-inversions. By choosing $d=n-1$ or $d=1,$ it is seen
that this class includes common inversions and descents. Moreover,
$d$-descents or generalized descents are given by

(3) $X_{\mathrm{des}}^{(d)}:=\sum_{i=1}^{n-d}\textbf{1}\\{Z_{i}>Z_{i+d}\\}\,.$

Generalized inversions were first introduced by de Mari & Shayman [5] to solve
a problem in algebraic geometry. The first stochastic considerations for the
random variables $X_{\mathrm{inv}}^{(d)}$ are due to Bona [2] and Pike [10],
who computed the mean and variance and proved a central limit theorem (CLT).

Symmetric groups belong to the class of finite irreducible Coxeter groups,
which have been classified by [4]. The concept of inversions and descents can
be transferred to these groups, see [1, Section 1.4]. Two other important
subfamilies of the finite irreducible Coxeter groups are the signed
permutation groups $B_{n}$ and the even-signed permutation groups $D_{n}$. The
group $B_{n}$ consists of all maps
$\pi\negthickspace:\\{1,\ldots,n\\}\longrightarrow\\{1,\ldots,n\\}\cup\\{-1,\ldots,-n\\}$
for which $|\pi|$ is a permutation. Such a map is called a signed permutation.
The groups $D_{n}$ are the subgroups of $B_{n}$ containing all signed
permutations with an even number of negative signs. The random number of
inversions on these groups can be represented with help of i.i.d. variables
$Z_{1},\ldots,Z_{n}\sim U(-1,1)$ as follows:

(4a) $\displaystyle X_{\mathrm{inv}}^{B}$ $\displaystyle=\sum_{1\leq i<j\leq
n}\textbf{1}\\{Z_{i}>Z_{j}\\}+\sum_{1\leq i<j\leq
n}\textbf{1}\\{-Z_{i}>Z_{j}\\}+\sum_{i=1}^{n}\textbf{1}\\{Z_{i}<0\\}\,,$ (4b)
$\displaystyle X_{\mathrm{inv}}^{D}$ $\displaystyle=\sum_{1\leq i<j\leq
n}\textbf{1}\\{Z_{i}>Z_{j}\\}+\sum_{1\leq i<j\leq
n}\textbf{1}\\{-Z_{i}>Z_{j}\\}\,.$

We refer to the three families $S_{n},B_{n},$ and $D_{n}$ as the classical
Weyl groups. Recently, Meier & Stump [9] extended the concept of generalized
inversions and descents to the groups $B_{n}$ and $D_{n}$, again proving a
CLT.

The investigation of extreme value asymptotics for common inversions and
descents was initiated in [7, 6], so we aim to extend this knowledge to
generalized inversions. We aim to prove Gumbel attraction for both the stand-
alone statistics $X_{\mathrm{inv}}^{(d_{1})},X_{\mathrm{des}}^{(d_{2})}$ and
the joint statistic
$(X_{\mathrm{inv}}^{(d_{1})},X_{\mathrm{des}}^{(d_{2})})^{\top}$. Here,
$d_{1}$ and $d_{2}$ can be either fixed or dependent on $n$. We will
investigate the impact of the choice of $d$ for these results.

This paper is structured as follows. Section 2 gives basic properties on
generalized inversions and descents, and introduces Hájek projections that are
used as independent sum approximations. Section 3 deduces the bivariate CLT
and the extreme value theorems for generalized inversions and descents.
Section 4 gives the proofs of two auxiliary lemmas determining the choices of
$d$ for which the Hájek projection serves as a working approximation. We use
typical Landau notation for positive sequences $a_{n},b_{n}$ as follows:

  * •

$a_{n}=O(b_{n})$ means that $\limsup_{n\rightarrow\infty}a_{n}/b_{n}<\infty$.

  * •

$a_{n}=o(b_{n})$ means that $\lim_{n\rightarrow\infty}a_{n}/b_{n}=0$. This is
also written as $a_{n}\ll b_{n}$ or $b_{n}\gg a_{n}$.

  * •

$a_{n}=\Theta(b_{n})$ means that $a_{n}$ and $b_{n}$ have the same order of
magnitude, i.e., both $a_{n}=O(b_{n})$ and $b_{n}=O(a_{n})$ hold.

  * •

$a_{n}=b_{n}+o_{\mathbb{P}}(1)$ means that $a_{n},b_{n}$ are sequences of
random variables with $a_{n}-b_{n}\overset{\mathbb{P}}{\longrightarrow}0$.

##  2\. Basic definitions

######  Definition 2.1.

Let $S_{n}$ be a symmetric group and let $\pi\in S_{n}$. For
$d\in\\{1,\ldots,n-1\\}$, $d$-inversions are all pairs $(i,j)$ in
$\mathfrak{N}_{n,d}$ with $\pi(i)>\pi(j)$. In this sense, common descents
equal $1$-inversions, and common inversions equal $(n-1)$-inversions.
Moreover, $d$-descents are all numbers $i\in\\{1,\ldots,n-d\\}$ with
$\pi(i)>\pi(i+d)$. We write $X_{\mathrm{inv}}^{(d)}$ for the random number of
$d$-inversions and $X_{\mathrm{des}}^{(d)}$ for the number of $d$-descents.
Probabilistic representations are given in (2), (3). 111In the literature,
there are different terminologies, e.g., in [2, 10], $d$-inversions are called
$d$-descents. However, we use the terms and notation provided in [9].

######  Remark 2.2.

Obviously, every number $k\in\\{1,\ldots,n\\}$ can appear in at most $2d$
$d$-inversions. This bound is redundant if $d>n/2$. In fact, it is an
important case distinction whether $d\leq n/2$ or $d>n/2$, e.g., when
calculating the mean and variance of $X_{\mathrm{inv}}^{(d)}$ and
$X_{\mathrm{des}}^{(d)}$. In the case of $d\leq n/2,$ we can split
$\\{1,\ldots,n\\}$ into the regions

$K_{1}:=\\{1,\ldots,d\\},\qquad K_{2}:=\\{d+1,\ldots,n-d\\},\qquad
K_{3}:=\\{n-d+1,\ldots,n\\}.$

For each $k\in K_{2},$ all larger indices $k+1,\ldots,k+d$ and all smaller
indices $k-1,\ldots,k-d$ allow to form $d$-inversions. For any $k\notin
K_{2},$ there are less than $d$ indices available in one direction. If $k\in
K_{1}$, then $k-1$ smaller indices are available for $d$-inversions, which is
less than $d$ numbers. If $k>n-d$ is large, then there are only $n-k<d$ larger
indices available.

On the contrary, if $d>n/2,$ then $n-d<d$ and the above partition into three
regions is now written as

$K_{1}:=\\{1,\ldots,n-d\\},\qquad K_{2}:=\\{n-d+1,\ldots,d\\},\qquad
K_{3}:=\\{d+1,\ldots,n\\}.$

Now, if $k\in K_{2},$ then there are only $k-1$ smaller indices and $n-k$
larger indices available to form $d$-inversions.

The mean and variance of $X_{\mathrm{inv}}^{(d)}$ have been extensively
computed by Pike [10]. It is easy to verify that the special cases $d=1$ and
$d=n-1$ are consistent with [8, Corollaries 3.2 and 4.2]. The proof of [10,
Theorem 1] is reviewed in [9, Theorem A.1], where the variance of $d$-descents
is provided as well.

######  Theorem 2.3 (see [10], Theorem 1 and [9], Theorem A.1).

For all $d=1,\ldots,n-1,$ it holds that

$\displaystyle\mathbb{E}\left(X_{\mathrm{inv}}^{(d)}\right)$
$\displaystyle=\frac{2nd-d^{2}-d}{4},$
$\displaystyle\mathbb{E}\left(X_{\mathrm{des}}^{(d)}\right)$
$\displaystyle=\frac{n-d}{2}.$

Moreover, if $d\leq n/2,$ then

$\displaystyle\mathrm{Var}\left(X_{\mathrm{inv}}^{(d)}\right)$
$\displaystyle=\frac{6nd+4d^{3}+3d^{2}-d}{72},$
$\displaystyle\mathrm{Var}\left(X_{\mathrm{des}}^{(d)}\right)$
$\displaystyle=\frac{n+d}{12}.$

If $d>n/2,$ then

$\mathrm{Var}\left(X_{\mathrm{inv}}^{(d)}\right)=-\frac{1}{6}d^{3}+\left(\frac{1}{3}n-\frac{7}{24}\right)d^{2}-\left(\frac{1}{6}n^{2}-\frac{5}{12}n+\frac{1}{8}\right)d+\frac{1}{36}n^{3}-\frac{1}{12}n^{2}+\frac{1}{18}n$

and $\mathrm{Var}\left(X_{\mathrm{des}}^{(d)}\right)=(n-d)/4$.

Meier & Stump [9] introduced an extension of generalized inversions and
descents from symmetric groups to the other classical Weyl groups $B_{n}$ and
$D_{n}$. This extension is based on the root poset of a classical Weyl group.
We refer to [9, Section 2] for the details. On the symmetric group $S_{n}$,
the ordered pairs of indices $(i,j)$ correspond to the positive roots
$[ij]:=e_{i}-e_{j}$, where $e_{i},e_{j}$ are unit vectors in $\mathbb{R}^{n}$,
and the height of $[ij]$ within the root poset is $\text{ht}([ij])=j-i$.

On the signed permutation group $B_{n},$ we also have to consider the positive
roots $[\widetilde{ij}]:=e_{i}+e_{j}$ and $[i]:=e_{i}$ for $1\leq i<j\leq n$.
The heights of these additional roots are $\text{ht}([\widetilde{ij}])=i+j$
and $\text{ht}([i])=i$. The root $[\widetilde{ij}]$ corresponds to the
indicator $\textbf{1}\\{-\pi(i)>\pi(j)\\}$ appearing in (4a), while $[i]$
corresponds to $\textbf{1}\\{\pi(i)<0\\}$. See [9, Example 2.2] for an
illustration of the root poset of $B_{n}$. On the even-signed permutation
group $D_{n},$ the roots $[i]$ are disregarded, and $[\widetilde{ij}]$ has
height $i+j-2$.

######  Definition 2.4.

For any classical Weyl group, $d$-inversions are determined by roots of height
at most $d$, and $d$-descents are determined by roots of height exactly $d$,
see [9, Definition 2.4]. For symmetric groups, this coincides with Definition
2.7. In addition to $\mathfrak{N}_{n,d}$, we introduce

$\widetilde{\mathfrak{N}}_{n,d}:=\\{(i,j)\in\\{1,\ldots,n\\}^{2}\mid
i<j,j+i\leq d\\}.$

Then, on the signed and even-signed permutation groups,
$X_{\mathrm{inv}}^{(d)}$ and $X_{\mathrm{des}}^{(d)}$ can be expressed as
follows:

(5a) $\displaystyle X_{\mathrm{inv}}^{B,(d)}$
$\displaystyle=\sum_{(i,j)\in\mathfrak{N}_{n,d}}\textbf{1}\\{Z_{i}>Z_{j}\\}+\sum_{(i,j)\in\widetilde{\mathfrak{N}}_{n,d}}\textbf{1}\\{-Z_{i}>Z_{j}\\}+\sum_{i=1}^{d}\textbf{1}\\{Z_{i}<0\\}\,,$
(5b) $\displaystyle X_{\mathrm{inv}}^{D,(d)}$
$\displaystyle=\sum_{(i,j)\in\mathfrak{N}_{n,d}}\textbf{1}\\{Z_{i}>Z_{j}\\}+\sum_{(i,j)\in\widetilde{\mathfrak{N}}_{n,d+2}}\textbf{1}\\{-Z_{i}>Z_{j}\\}\,,$
(5c) $\displaystyle X_{\mathrm{des}}^{B,(d)}$
$\displaystyle=\sum_{i=1}^{n-d}\textbf{1}\\{Z_{i}>Z_{i+d}\\}+\sum_{i=1}^{\lceil
d/2\rceil-1}\textbf{1}\\{-Z_{i}>Z_{d-i}\\}+\textbf{1}\\{\pi(i)>d\\}\,,$ (5d)
$\displaystyle X_{\mathrm{des}}^{D,(d)}$
$\displaystyle=\sum_{i=1}^{n-d}\textbf{1}\\{Z_{i}>Z_{i+d}\\}+\sum_{i=1}^{\lceil
d/2\rceil}\textbf{1}\\{-Z_{i}>Z_{d+2-i}\\}\,.$

######  Remark 2.5.

The largest possible choice of $d$ equals the total height of the root poset,
namely, $d_{\max}-1$, where $d_{\max}$ denotes the largest degree of the
underlying classical Weyl group. In particular, $d_{\max}=n$ for $S_{n},$
$d_{\max}=2n$ for $B_{n}$ and $d_{\max}=2n-2$ for $D_{n}$. To precisely
compute the variance of $X_{\mathrm{inv}}^{(d)}$ and $X_{\mathrm{des}}^{(d)}$
on the groups $B_{n}$ and $D_{n},$ one needs to distinguish eight cases, as
seen in [9, Theorems A.4 and A.13]. However, many of these cases give the same
asymptotic quantification, which can be stated as follows:

######  Lemma 2.6 (cf. [9], Theorems A.4 and A.13).

For the generalized inversions and descents on both the groups $B_{n}$ and
$D_{n},$ it holds that

$\displaystyle\mathrm{Var}\left(X_{\mathrm{inv}}^{(d)}\right)$
$\displaystyle=\begin{cases}\frac{1}{36}d^{3}+\frac{1}{12}nd+O(d^{2}),&d\leq
n/2\\\ \frac{1}{36}d^{3}+O(d^{2}),&n/2\leq d<n\\\
-\frac{1}{12}d^{3}+\frac{1}{3}nd^{2}-\frac{1}{3}n^{2}d+\frac{1}{9}n^{3}+O(d^{2}),&d\geq
n\end{cases},$ $\displaystyle\mathrm{Var}\left(X_{\mathrm{des}}^{(d)}\right)$
$\displaystyle=\begin{cases}\frac{1}{24}d+\frac{1}{12}n+O(1),&d<n\\\
-\frac{1}{8}d+\frac{1}{4}n+O(1),&d\geq n\end{cases}.$

An evident issue seen in (2), (5a), and (5b) is the mutual dependence of the
indicators $\textbf{1}\\{Z_{i}>Z_{j}\\}$ and $\textbf{1}\\{-Z_{i}>Z_{j}\\}$.
Therefore, we aim to approximate $X_{\mathrm{inv}}^{(d)}$ with a sum of
independent variables. Hájek projections are a common tool for this.

######  Definition 2.7.

Let $Z_{1},\ldots,Z_{n}$ be independent random variables, and let $X$ be
another random variable. Then, the Hájek projection of $X$ with respect to
$Z_{1},\ldots,Z_{n}$ is given by

$\hat{X}:=\sum_{k=1}^{n}\mathbb{E}(X\mid Z_{k})-(n-1)\mathbb{E}(X).$

As every $\mathbb{E}(X\mid Z_{k})$ is a measurable function only in $Z_{k}$,
the Hájek projection is a sum of independent random variables. If
$(X_{n})_{n\in\mathbb{N}}$ is a sequence of random variables with each $X_{n}$
depending on $Z_{1},\ldots,Z_{n}$, then the Hájek projections $\hat{X}_{n}$
should give a sufficiently accurate approximation to $X_{n}$. To ensure this,
the following criterion is useful.

######  Theorem 2.8 (cf. [11], Theorem 11.2).

Consider a sequence $(X_{n})_{n\geq 1}$ of random variables and their
associated Hájek projections $(\hat{X}_{n})_{n\geq 1}$. If
$\mathrm{Var}(\hat{X}_{n})\sim\mathrm{Var}(X_{n})$ as $n\rightarrow\infty,$
then

$\frac{X_{n}-\mathbb{E}(X_{n})}{\mathrm{Var}(X_{n})^{1/2}}=\frac{\hat{X}_{n}-\mathbb{E}(\hat{X}_{n})}{\mathrm{Var}(\hat{X}_{n})^{1/2}}+o_{\mathbb{P}}(1).$

##  3\. Asymptotic results

### Asymptotic normality

For the univariate statistics $X_{\mathrm{inv}}^{(d)}$ and
$X_{\mathrm{des}}^{(d)}$, the asymptotic normality has been proven in [9]. For
the joint statistic
$\left(X_{\mathrm{inv}}^{(d)},X_{\mathrm{des}}^{(d)}\right)^{\top}$, the
simplest case is when $d$ remains fixed. Then, both $X_{\mathrm{inv}}^{(d)}$
and $X_{\mathrm{des}}^{(d)}$ are $m$-dependent $(m=d)$ and it is not necessary
to apply the Hájek approximation. From the CLT for $m$-dependent random
vectors, it follows that:

######  Theorem 3.1.

For any two fixed numbers $d_{1},d_{2},$ the joint distribution
$\left(X_{\mathrm{inv}}^{(d_{1})},X_{\mathrm{des}}^{(d_{2})}\right)^{\top}$
satisfies the CLT.

An even stronger statement than the CLT is the uniform Gaussian approximation
over all hyperrectangles for $m$-dependent random vectors given in [3, Theorem
2]. It also applies for random vectors with a sparse dependency graph. Let
$\left(X_{j}^{(n)}\right)_{j=1,\ldots,n}$ be a triangular array of (w.l.o.g.,
centered) $p$-dimensional random vectors and let

$\displaystyle X^{(n)}$ $\displaystyle:=\sum_{t=1}^{n}X_{t}^{(n)},$
$\displaystyle\Sigma^{(n)}$ $\displaystyle:=\mathrm{Var}(X^{(n)})\,,$
$\displaystyle\mathcal{N}$
$\displaystyle\sim\mathrm{N}_{p}(0,\Sigma^{(n)})\,.$

Let $\mathcal{A}^{\mathrm{re}}$ be the system of all $p$-dimensional
hyperrectangles, including infinite bounds, i.e.,
$\mathcal{A}^{\mathrm{re}}:=\bigl{\\{}\\{\textbf{w}\in\mathbb{R}^{p}\negmedspace:\textbf{a}\leq\textbf{w}\leq\textbf{b}\\}\mid\textbf{a},\textbf{b}\in[-\infty,\infty]^{p}\bigr{\\}}\,.$
Let

$r_{n}(\mathcal{A}^{\mathrm{re}}):=\sup_{A\in\mathcal{A}^{\mathrm{re}}}|\mathbb{P}(X^{(n)}\in
A)-\mathbb{P}(\mathcal{N}\in A)|\,.$

In particular, $r_{n}(\mathcal{A}^{\mathrm{re}})$ is an upper bound for

$r_{n}(\mathcal{A}^{\mathrm{CLT}}):=\sup_{u\in\mathbb{R}^{p}}|\mathbb{P}(X^{(n)}\leq
u)-\mathbb{P}(\mathcal{N}_{n}\leq u)|\,,$

which is relevant for the CLT, and for

$r_{n}(\mathcal{A}^{\mathrm{ext}}):=\sup_{u\in\mathbb{R}^{p}}|\mathbb{P}(X^{(n)}>u)-\mathbb{P}(\mathcal{N}_{n}>u)|\,,$

which is relevant for the asymptotics of extreme values. Now, let $G_{n}$ be a
dependency graph for $X_{1}^{(n)},\ldots,X_{n}^{(n)},$ which consists of all
edges $(i,j)$ for which $X_{i}^{(n)}$ and $X_{j}^{(n)}$ are dependent. Let
$\Delta_{n}$ be the maximum degree of $G_{n}$ and let $\Delta_{n}^{*}$ be the
maximum degree of the 2-reachability graph of $G_{n}$. If the graphs $G_{n}$
are not too dense as $n\rightarrow\infty,$ then
$r_{n}(\mathcal{A}^{\mathrm{re}})$ can be bounded as follows:

######  Theorem 3.2 (cf. [3], Theorem 2).

Let $(X_{j}^{(n)})_{j=1,\ldots,n},$ $r_{n}(\mathcal{A}^{\mathrm{re}}),$
$\Delta_{n},\Delta_{n}^{*}$ be as above. Under some conditions which are
satisfied by bounded variables, it holds that

$r_{n}(\mathcal{A}^{\mathrm{re}})=O\left(\frac{(\Delta_{n}\Delta_{n}^{*})^{1/3}\log(p)^{7/6}}{n^{1/6}}\right).$

In particular, if $(X_{j}^{(n)})_{j=1,\ldots,n}$ is $m$-dependent for some
global constant $m\in\mathbb{N},$ then
$r_{n}(\mathcal{A}^{\mathrm{re}})=O\bigl{(}n^{-1/6}\log(p)^{7/6}\bigr{)}$.

This theorem was stated by Chang et al. for high dimensions, but it also works
with fixed $p$ by artificially repeating the components of a random vector in
fixed dimension.

From now, we consider the univariate statistics
$X_{\mathrm{inv}}^{(d_{1})},X_{\mathrm{des}}^{(d_{2})}$ and the joint
statistic
$\left(X_{\mathrm{inv}}^{(d_{1})},X_{\mathrm{des}}^{(d_{2})}\right)^{\top}$
for two sequences $d_{1}=d_{1}(n),$ $d_{2}=d_{2}(n)$ with
$d_{1}(n),d_{2}(n)\leq d_{\max}-1$ $\forall n\in\mathbb{N}$. Both
$X_{\mathrm{inv}}^{(d_{1})}$ and $X_{\mathrm{des}}^{(d_{2})}$ are based on a
sequence of classical Weyl groups $(W_{n})_{n\in\mathbb{N}},$ where each
$W_{n}$ is one of $S_{n},B_{n},$ or $D_{n}$.

If $d_{2}=d_{2}(n)$ grows as $n\rightarrow\infty$, then there exists no
constant $m\in\mathbb{N}$ for which all $X_{\mathrm{des}}^{(d_{2})}$ are
$m$-dependent. However, the dependency structure of
$X_{\mathrm{des}}^{(d_{2})}$ is still sparse. Recall that
$X_{\mathrm{inv}}^{(d_{1})}$ and $X_{\mathrm{des}}^{(d_{2})}$ are constructed
from i.i.d. random variables $Z_{1},\ldots,Z_{n}$. According to (3), (5c),
(5d), we can represent $X_{\mathrm{des}}^{(d_{2})}$ as a sum of indicator
variables, each of which depends on at most three others. Therefore, the
first-order and second-order maximum degrees $\Delta_{n},\Delta_{n}^{*}$ are
bounded in the way of $\Delta_{n}\leq 3$ and $\Delta_{n}^{*}\leq 9$. In
conclusion, there are no issues with applying Theorem 3.2 if $d_{2}$ grows.

If $d_{1}=d_{1}(n)$ grows as well, then the maximum degrees
$\Delta_{n},\Delta_{n}^{*}$ of
$\left(X_{\mathrm{inv}}^{(d_{1})},X_{\mathrm{des}}^{(d_{2})}\right)^{\top}$
are bounded in the way of $\Delta_{n}\leq 4d_{1}$ and $\Delta_{n}^{*}\leq
8d_{1}$. Moreover, we have to take into account that by (2), (5a), (5b),
$X_{\mathrm{inv}}^{(d_{1})}$ is based on $\Theta(nd_{1})$ summands, so we have
to replace $n$ with $nd$ in Theorem 3.2. Therefore, Theorem 3.2 gives an
$o(1)$ bound of $r_{n}(\mathcal{A}^{\mathrm{re}})$ if and only if
$(d_{1})^{2/3}=o((nd_{1})^{1/6})\Longleftrightarrow d_{1}=o(n^{1/3})$. For any
faster growth rate of $d_{1},$ the dependency structure of
$\left(X_{\mathrm{inv}}^{(d_{1})},X_{\mathrm{des}}^{(d_{2})}\right)^{\top}$ is
too complex to apply Theorem 3.2. In that case, we have to replace
$X_{\mathrm{inv}}^{(d_{1})}$ with the Hájek projection
$\hat{X}_{\mathrm{inv}}^{(d_{1})}$. The growth rate of $d_{1}$ determines
whether the condition in Theorem 2.8 is fulfilled or not. To classify the
validity of this condition, we distinguish the short case $d\leq d_{\max}/2$
and the long case $d>d_{\max}/2$.

######  Lemma 3.3.

Given the short case, we have for symmetric groups that

$\mathrm{Var}\left(\hat{X}_{\mathrm{inv}}^{(d)}\right)=\frac{1}{72}\left(4d^{3}+6d^{2}+2d\right)\,,$

and we have
$\mathrm{Var}\left(\hat{X}_{\mathrm{inv}}^{(d)}\right)=d^{3}/36+O(d^{2})$ for
the groups $B_{n}$ and $D_{n}$. So, if $d_{\max}/2>d\gg\sqrt{n},$ then
$\mathrm{Var}\left(\hat{X}_{\mathrm{inv}}^{(d)}\right)\sim\mathrm{Var}\left(X_{\mathrm{inv}}^{(d)}\right)$.

######  Lemma 3.4.

Given the long case,
$\mathrm{Var}\left(\hat{X}_{\mathrm{inv}}^{(d)}\right)\sim\mathrm{Var}\left(X_{\mathrm{inv}}^{(d)}\right)$
holds on symmetric groups if and only if $d\sim d_{\max}$. On the other
classical Weyl groups,
$\mathrm{Var}\left(\hat{X}_{\mathrm{inv}}^{(d)}\right)\sim\mathrm{Var}\left(X_{\mathrm{inv}}^{(d)}\right)$
holds if and only if $d\sim d_{\max}$ or $d\sim d_{\max}/2$.

From Theorem 3.2 and Slutsky’s lemma, we can conclude:

######  Corollary 3.5.

$\left(X_{\mathrm{inv}}^{(d_{1})},X_{\mathrm{des}}^{(d_{2})}\right)^{\top}$
satisfies the CLT if $d_{1}=o(n^{1/3})$ or if $d_{1}$ satisfies the conditions
of either Lemma _3.3_ or Lemma _3.4_.

### Extreme values of the univariate statistics $X_{\mathrm{inv}}^{(d)}$ and
$X_{\mathrm{des}}^{(d)}$

We now postulate the univariate and bivariate extreme value limit theorems
(EVLTs) for $X_{\mathrm{inv}}^{(d)}$ and $X_{\mathrm{des}}^{(d)}$ with the
help of the methods used in [6]. The proof of [6, Theorem 4.1] consists of two
parts. The first part is to show the Gumbel attraction of the statistic that
replaces $X_{\mathrm{inv}}^{(d)}$ with $\hat{X}_{\mathrm{inv}}^{(d)}$. This is
done by an application of [3, Theorem 2], providing Gaussian approximation for
$m$-dependent random vectors. The second part is to show that the error
resulting from this replacement vanishes in probability. Recall that by use of

$\displaystyle\alpha_{n}$ $\displaystyle=\frac{1}{\sqrt{2\log n}}\,,$
$\displaystyle\beta_{n}$
$\displaystyle=\frac{1}{\alpha_{n}}-\frac{1}{2}\alpha_{n}\bigl{(}\log\log
n+\log(4\pi)\bigr{)}\,,$

the maximum of $n$ i.i.d. variables $N_{1},\ldots,N_{n}\sim\mathrm{N}(0,1)$ is
attracted to the Gumbel distribution, that is, $\forall x\in\mathbb{R}$:

$\mathbb{P}\left(\max_{i=1,\ldots,n}N_{i}\leq\alpha_{n}x+\beta_{n}\right)\longrightarrow\exp(-\exp(-x))=:\Lambda(x).$

Likewise, the two-dimensional standard normal distribution is attracted to the
two-dimensional Gumbel distribution with independent marginals. Taking
$\boldsymbol{\upalpha}_{n}:=(\alpha_{n},\alpha_{n}),$
$\boldsymbol{\upbeta}_{n}:=(\beta_{n},\beta_{n})\in\mathbb{R}^{2}$ and writing
”$*$” for component-wise multiplication, we get that for i.i.d. bivariate
standard normal $\mathbf{N}_{1},\ldots,\mathbf{N}_{n}$ and
$\mathbf{x}=(x_{1},x_{2})\in\mathbb{R}^{2}$:

$\mathbb{P}\left(\max_{i=1,\ldots,n}\mathbf{N}_{i}\leq\boldsymbol{\upalpha}_{n}*\mathbf{x}+\boldsymbol{\upbeta}_{n}\right)\longrightarrow\exp\bigl{(}-\exp(-x_{1})-\exp(-x_{2})\bigr{)}=:\Lambda_{2}(x).$

Since the random numbers of generalized inversions and descents are discrete
distributions, we are interested in the extremes of a triangular array
$(X_{n1},\ldots,X_{nk_{n}}),$ where for each $n\in\mathbb{N},$ the block
$X_{n1},\ldots,X_{nk_{n}}$ consists of i.i.d. samples taken from the $n$-th
symmetric group $S_{n}$. Each of the following EVLTs imposes an upper bound on
the number of samples $k_{n}$.

For a univariate triangular array consisting of generalized descents, it is
not necessary to use the Hájek projection. It has already been argued in [6,
Remark 4.2] that a subexponential bound of $k_{n}$ can be obtained if the
Hájek projection is not needed. So, in the univariate EVLT for generalized
descents, we can use the bound of $k_{n}=\exp(o(n^{1/7}))$ stated in [6,
Remark 4.2]:

######  Theorem 3.6.

Let $(X_{nj})_{j=1,\ldots,k_{n}}$ be a row-wise i.i.d. triangular array with
$X_{n1}\overset{\mathcal{D}}{=}X_{\mathrm{des}}^{(d)}$ for any sequence
$d=d(n)$ with $1\leq d\leq d_{\max}-1,$ and let
$M_{n}:=\max\\{X_{n1},\ldots,X_{nk_{n}}\\}$. Let
$a_{n}:=\sigma\bigl{(}X_{\mathrm{des}}^{(d)}\bigr{)}\alpha_{k_{n}}$ and let
$b_{n}:=\sigma\bigl{(}X_{\mathrm{des}}^{(d)}\bigr{)}\beta_{k_{n}}+\mu\bigl{(}X_{\mathrm{inv}}^{(d)}\bigr{)}.$
If $k_{n}=\exp(o(n^{1/7})),$ then

$\forall x\in\mathbb{R}\negmedspace:\quad\mathbb{P}(M_{n}\leq
a_{n}x+b_{n})\longrightarrow\Lambda(x)\,.$

A similar statement applies to $X_{\mathrm{inv}}^{(d)}$ if $d$ grows slow
enough to permit the application of Theorem 3.2.

######  Theorem 3.7.

Let $(X_{n1},\ldots,X_{nk_{n}})$ be a row-wise i.i.d. triangular array with
$X_{n1}\overset{\mathcal{D}}{=}X_{\mathrm{inv}}^{(d)}$, where
$d=d(n)=o(n^{1/3})$. Let $M_{n},a_{n},b_{n}$ be given in analogy to Theorem
_3.6_. If $k_{n}=\exp\bigl{(}o(n^{1/7}/d^{3/7})\bigr{)},$ then

$\forall x\in\mathbb{R}\negmedspace:\quad\mathbb{P}(M_{n}\leq
a_{n}x+b_{n})\longrightarrow\Lambda(x)\,.$

###### Proof.

According to the above considerations, the maximum degrees
$\Delta_{n},\Delta_{n}^{*}$ in the dependency graph of the representations
(2), (5a), and (5b) are bounded in the way of $\Delta_{n}\leq 4d$ and
$\Delta_{n}^{*}\leq 8d$. Since these representations are based on $\Theta(nd)$
summands, we need to replace $n$ with $nd$ when applying Theorem 3.2. An
application of Theorem 3.2 with $\max\\{n,k_{n}\\}$ i.i.d. iterations of
$X_{\mathrm{inv}}^{(d)}$ yields

$|\mathbb{P}(M_{n}\leq
a_{n}x+b_{n})-\mathbb{P}(\mathcal{M}_{n}\leq\alpha_{n}x+\beta_{n})|=O(n^{-1/6}d^{1/2})\log(k_{n})^{7/6}=o(1)\,,$

where $\mathcal{M}_{n}$ is the maximum of $n$ i.i.d. copies of the standard
normal distribution. The claim follows. ∎

For any other growth rate of $d,$ we can state an EVLT only for the cases
covered by Lemmas 3.3 and 3.4, and due to the use of Hájek’s projection, we
can only impose a strongly reduced asymptotic bound for $k_{n}$.

######  Theorem 3.8.

Let $(X_{n1},\ldots,X_{nk_{n}})$ be a row-wise i.i.d. triangular array with
$X_{n1}\overset{\mathcal{D}}{=}X_{\mathrm{inv}}^{(d)}$, where $d=d(n)$
satisfies the conditions in either Lemma _3.3_ or _3.4_. Assume
$k_{n}\log(k_{n})=o(n^{-1}d^{2})$ in the short case and
$k_{n}\log(k_{n})=o\left(\min\left\\{n,\displaystyle{\frac{n}{n-d}}\right\\}\right)$
in the long case. Let $M_{n}:=\max\\{X_{n1},\ldots,X_{nk_{n}}\\}$ be the row-
wise maximum. Then,

$\forall x\in\mathbb{R}\negthickspace:\quad\mathbb{P}(M_{n}\leq
a_{n}x+b_{n})\longrightarrow\Lambda(x)\,,$

with $a_{n},b_{n}$ given in analogy to Theorem _3.6_.

###### Proof.

The conditions on $d$ ensure that
$1-\mathrm{Var}(X_{\mathrm{inv}}^{(d)})/\mathrm{Var}(\hat{X}_{\mathrm{inv}}^{(d)})=o(1),$
but the rate of convergence determines the bound for $k_{n}$ by means of [6,
Eq. (9), (10)]. We compute this rate for symmetric groups, since the same
conclusions can be obtained on the other classical Weyl groups. In the short
case, we have

$\displaystyle
1-\frac{\mathrm{Var}(X_{\mathrm{inv}}^{(d)})}{\mathrm{Var}(\hat{X}_{\mathrm{inv}}^{(d)})}$
$\displaystyle=1-\frac{4d^{3}+6nd+3d^{2}-d}{4d^{3}+6d^{2}+2d}=\frac{4d^{2}+6n+3d-1}{4d^{2}+6d+2}$
$\displaystyle=1-\frac{4d^{2}}{4d^{2}+6d+2}-\frac{3d+1}{4d^{2}+6d+2}-\frac{6n}{4d^{2}+6d+2}$
$\displaystyle=\Theta\left(\frac{1}{d}\right)-\frac{6n}{4d^{2}+6d+2}.$

Apparently, $6n/(4d^{2}+6d+2)$ always dominates $1/d$. In conclusion,

$1-\frac{\mathrm{Var}(X_{\mathrm{inv}}^{(d)})}{\mathrm{Var}(\hat{X}_{\mathrm{inv}}^{(d)})}=\Theta\left(\frac{n}{d^{2}}\right),$

giving the condition of $k_{n}\log(k_{n})=o(d^{2}/n)$ according to the
arguments in the proof of [6, Theorem 4.1]. The remaining steps of the proof
are the same.

In the long case, we write $\delta:=n-d$ and note that

$\displaystyle\mathrm{Var}(\hat{X}_{\mathrm{inv}})$
$\displaystyle=-\frac{1}{6}(n-\delta)^{3}+\frac{1}{3}(n-\delta)^{2}n-\frac{5}{36}n^{3}+\frac{16}{36}n^{2}\delta+O(n^{2}),$
$\displaystyle\mathrm{Var}(X_{\mathrm{inv}})$
$\displaystyle=-\frac{1}{6}(n-\delta)^{3}+\frac{1}{3}(n-\delta)^{2}n-\frac{5}{36}n^{3}+\frac{6}{36}n^{2}\delta+O(n^{2}).$

Therefore,

$1-\frac{\mathrm{Var}(X_{\mathrm{inv}})}{\mathrm{Var}(\hat{X}_{\mathrm{inv}})}=\frac{10n^{2}\delta+O(n^{2})}{\Theta(n^{3})}=\Theta\left(\max\left\\{\frac{\delta}{n},\frac{1}{n}\right\\}\right),$

giving the requirement of

$k_{n}\log(k_{n})\Theta\left(\max\left\\{\frac{\delta}{n},\frac{1}{n}\right\\}\right)=o(1)\Longleftrightarrow
k_{n}\log(k_{n})=o\left(\min\left\\{n,\frac{n}{\delta}\right\\}\right).$

Again, the proof now follows the same steps as in [6, Theorem 4.1]. ∎

######  Remark 3.9.

This result is only useful if $d\gg n^{3/4}$. Otherwise, we can obtain the
EVLT with a bound of $k_{n}=o(n^{1/2})$ from the universal EVLT [7, Theorem
5.1].

### Extreme values of the joint distribution of $X_{\mathrm{inv}}^{(d)}$ and
$X_{\mathrm{des}}^{(d)}$

The bivariate EVLT for
$\left(X_{\mathrm{inv}}^{(d_{1})},X_{\mathrm{des}}^{(d_{2})}\right)^{\top}$
follows the same way, since the descent component does not interfere with any
of the previously used arguments.

######  Theorem 3.10.

Let $(X_{n1},\ldots,X_{nk_{n}})$ be a row-wise i.i.d. triangular array with
$X_{n1}\overset{\mathcal{D}}{=}\left(X_{\mathrm{inv}}^{(d_{1})},X_{\mathrm{des}}^{(d_{2})}\right)^{\top}$
on $S_{n}.$

  1. _(a)_

If $d_{1}=o(n^{1/3}),$ then assume
$k_{n}=\exp\bigl{(}o(n^{1/7}/d^{3/7})\bigr{)}.$

  2. _(b)_

If $d_{1}$ satisfies the conditions in either Lemma _3.3_ or _3.4_ , then
assume $k_{n}\log(k_{n})=o(n^{-1}d_{1}^{2})$ in the short case and
$k_{n}\log(k_{n})=o\left(\min\left\\{n,\displaystyle{\frac{n}{n-d_{1}}}\right\\}\right)$
in the long case.

Let $M_{n}:=\max\\{X_{n1},\ldots,X_{nk_{n}}\\}$ be the row-wise maximum. Let
$\mu_{n}:=\mathbb{E}(X_{n1})$ and let
$s_{n}:=\bigl{(}\sigma(X_{\mathrm{inv}}^{(d)}),\sigma(X_{\mathrm{des}})\bigr{)}$.
Let $a_{n}:=s_{n}*\boldsymbol{\upalpha}_{k_{n}}$ and let
$b_{n}:=s_{n}*\boldsymbol{\upbeta}_{k_{n}}+\mu_{n}$. Then,

$\forall\mathbf{x}\in\mathbb{R}^{2}\negthickspace:\quad\mathbb{P}(M_{n}\leq
a_{n}*\mathbf{x}+b_{n})\longrightarrow\Lambda_{2}(\mathbf{x})\,.$

##  4\. Proofs of Lemmas 3.3 and 3.4

We first consider $\hat{X}_{\mathrm{inv}}^{(d)}$ on symmetric groups and then
provide the analogous observations on the other classical Weyl groups.

###### Proof of Lemma 3.3 for symmetric groups.

By (2), we have

$\mathbb{E}(X_{\mathrm{inv}}^{(d)}\mid
Z_{k})=\sum_{(i,j)\in\mathfrak{N}_{n,d}}\mathbb{P}(Z_{i}>Z_{j}\mid
Z_{k})=\sum_{(i,j)\in\mathfrak{N}_{n,d}}\begin{cases}1/2,&k\notin\\{i,j\\}\\\
Z_{k},&k=i\\\ 1-Z_{k},&k=j\end{cases}.$

Only the pairs $(i,j)$ with $k\in\\{i,j\\}$ contribute to
$\mathrm{Var}\left(\mathbb{E}\left(X_{\mathrm{inv}}^{(d)}\mid
Z_{k}\right)\right)$. These contributions are called the non-trivial parts for
simplicity. The number of these pairs depends on whether $k$ belongs to
$K_{1},K_{2},$ or $K_{3}$, as already seen in Remark 2.2. Figure 1 visualizes
this case distinction for the exemplary choice of $n=15$ and $d=4$.

Figure 1. Overview of relevant pairs $(i,j)$ for computing the variance of
$\mathbb{E}\left(X_{\mathrm{inv}}^{(d)}\mid Z_{k}\right)$, where $n=15$ and
$d=4$. For each of the regions $K_{1},K_{2},K_{3},$ an exemplary index $k$ is
chosen, and the pairs that give a non-trivial contribution are highlighted in
red for $k\in K_{1}$, in blue for $k\in K_{2}$, and in green for $k\in K_{3}$.

If $k\in K_{2},$ then the non-trivial parts are

$\underbrace{Z_{k}+Z_{k}+\ldots+Z_{k}}_{d~{}\text{times}}+\underbrace{(1-Z_{k})+(1-Z_{k})+\ldots+(1-Z_{k})}_{d~{}\text{times}}=d.$

This means $\mathbb{E}\left(X_{\mathrm{inv}}^{(d)}\mid Z_{k}\right)$ is
constant due to cancellation, and vanishes when computing the variance. So, in
the short case, $\mathrm{Var}\left(\mathbb{E}\left(X_{\mathrm{inv}}^{(d)}\mid
Z_{k}\right)\right)$ originates only from $K_{1}$ and $K_{3}$. If $k\in
K_{1},$ then $\mathbb{E}\left(X_{\mathrm{inv}}^{(d)}\mid
Z_{k}\right)=(k-1)(1-Z_{k})+dZ_{k}+\text{const}=(d+1-k)Z_{k}+\text{const}$. If
$k\in K_{3}$, then $\mathbb{E}\left(X_{\mathrm{inv}}^{(d)}\mid
Z_{k}\right)=(n-k)Z_{k}+d(1-Z_{k})+\text{const}=(n-d-k)Z_{k}+\text{const}$.
So, the overall representation of $\hat{X}_{\mathrm{inv}}^{(d)}$ in the short
case is

$\hat{X}_{\mathrm{inv}}^{(d)}=\sum_{k=1}^{n}\omega_{d}(k)Z_{k}+\text{const}\,,\qquad\text{with}\quad\omega_{d}(k):=\begin{cases}d-k+1\,,&k\in
K_{1}\\\ 0,&k\in K_{2}\\\ n-d-k\,,&k\in K_{3}\end{cases}.$

Therefore,

$\displaystyle\mathrm{Var}\left(\hat{X}_{\mathrm{inv}}^{(d)}\right)$
$\displaystyle=\mathrm{Var}\left(\sum_{k=n-d+1}^{n}(n-k-d)Z_{k}+\sum_{k=1}^{d}(d+1-k)Z_{k}\right)$
$\displaystyle=\sum_{k=n-d+1}^{n}\frac{1}{12}(n-k-d)^{2}+\sum_{k=1}^{d}\frac{1}{12}(d+1-k)^{2}$
$\displaystyle=\frac{1}{12}\sum_{k=1}^{d}k^{2}+\frac{1}{12}\sum_{k=1}^{d}(-k)^{2}=\frac{1}{6}\frac{d(d+1)(2d+1)}{6}$
$\displaystyle=\frac{1}{72}\left(4d^{3}+6d^{2}+2d\right).$

So, the leading term is always $4d^{3}/72$. In light of Theorem 2.3, according
to which $\mathrm{Var}\left(X_{\mathrm{inv}}^{(d)}\right)$ contains the
monomials $6nd/72$ and $4d^{3}/72,$ it must be ensured that $d^{3}$ is
dominant over $nd$. This means $d\gg\sqrt{n}$. ∎

Now, we consider the long case. Due to $d>n/2,$ we now have $n-d<d$, and the
regions $K_{1},K_{2},K_{3}$ are redefined according to Remark 2.2. For the
non-trivial parts, we note that:

  * •

If $k\in K_{1}$ or $k\in K_{3},$ then the non-trivial parts are the same as in
the short case.

  * •

If $k\in K_{2},$ then the non-trivial parts yield
$\mathbb{E}\left(X_{\mathrm{inv}}^{(d)}\mid
Z_{k}\right)=(n-2k+1)Z_{k}+\text{const}$.

So, we again obtain a representation in the way of

(6)
$\hat{X}_{\mathrm{inv}}^{(d)}=\sum_{k=1}^{n}\omega_{d}(k)Z_{k}+\text{const}\,,\qquad\text{with}\quad\omega_{d}(k):=\begin{cases}d-k+1,&k\in
K_{1}\\\ n-2k+1,&k\in K_{2}\\\ n-d-k,&k\in K_{3}\end{cases}.$

The proof of the following lemma shows that if $d\not\sim n,$ then the leading
terms of $\mathrm{Var}\left(X_{\mathrm{inv}}^{(d)}\right)$ and
$\mathrm{Var}\left(\hat{X}_{\mathrm{inv}}^{(d)}\right)$ do not match
perfectly.

###### Proof of Lemma 3.4 for symmetric groups.

From (6), we state that

(6a) $\displaystyle\mathrm{Var}\left(\hat{X}_{\mathrm{inv}}^{(d)}\right)$
$\displaystyle=\frac{1}{12}\sum_{k=1}^{n-d}(d+1-k)^{2}$ (6b)
$\displaystyle+\frac{1}{12}\sum_{k=n-d+1}^{d}(n-2k+1)^{2}$ (6c)
$\displaystyle+\frac{1}{12}\sum_{k=d+1}^{n}(n-k-d)^{2}.$

By appropriate index shifting, we calculate

$\displaystyle\eqref{6a},\eqref{6c}$
$\displaystyle=\frac{1}{12}\left(\sum_{k=1}^{d}k^{2}-\sum_{k=1}^{2d-n}k^{2}\right)$
$\displaystyle\Longrightarrow\eqref{6a}+\eqref{6c}$
$\displaystyle=\frac{1}{6}\left(\sum_{k=1}^{d}k^{2}-\sum_{k=1}^{2d-n}k^{2}\right)$
$\displaystyle=\frac{1}{36}(n-d)\bigl{(}14d^{2}+d(9-10n)+12n^{2}-3n+1\bigr{)},$
$\displaystyle\eqref{6b}$
$\displaystyle=\frac{1}{36}(2d-n)(4d^{2}-4dn+n^{2}-1).$

This gives the total result

$\displaystyle\mathrm{Var}\left(\hat{X}_{\mathrm{inv}}^{(d)}\right)=\frac{1}{36}\left(-6d^{3}+12d^{2}n-9d^{2}-16dn^{2}+12dn-3d+11n^{3}-3n^{2}+2n\right).$

In contrast, by Theorem 2.3,

$\mathrm{Var}\left(X_{\mathrm{inv}}^{(d)}\right)=\frac{1}{36}\left(-6d^{3}+\left(12n-\frac{21}{2}\right)d^{2}-\left(6n^{2}-15n+\frac{9}{2}\right)d+n^{3}-3n^{2}+2n\right).$

Since the long case implies $n/2<d<n$, all monomials of order 3 are leading
terms. The monomials $-16dn^{2}+11n^{3}$ appearing in
$\mathrm{Var}\left(\hat{X}_{\mathrm{inv}}^{(d)}\right)$ and $-6n^{2}d+n^{3}$
appearing in $\mathrm{Var}\left(X_{\mathrm{inv}}^{(d)}\right)$ do not match in
general. Using $\delta:=n-d$, we can write

$\displaystyle\mathrm{Var}(\hat{X}_{\mathrm{inv}}^{(d)})$
$\displaystyle=-\frac{1}{6}(n-\delta)^{3}+\frac{1}{3}(n-\delta)^{2}n-\frac{5}{36}n^{3}+\frac{16}{36}n^{2}\delta+O(n^{2}),$
$\displaystyle\mathrm{Var}(X_{\mathrm{inv}}^{(d)})$
$\displaystyle=-\frac{1}{6}(n-\delta)^{3}+\frac{1}{3}(n-\delta)^{2}n-\frac{5}{36}n^{3}+\frac{6}{36}n^{2}\delta+O(n^{2}).$

This shows that we need $\delta=o(n),$ since otherwise, there are two
different leading coefficients of $n^{3}$. ∎

We now derive the analogous statements for the other classical Weyl groups
$B_{n}$ and $D_{n}$. It is sufficient to prove the lemmas 3.3 and 3.4 for the
groups $B_{n},$ since the difference between $X_{\mathrm{inv}}^{B,(d)}$ and
$X_{\mathrm{inv}}^{D,(d)}$ is asymptotically negligible (cf. (5a) and (5b)).
Recall the asymptotic quantification of $\mathrm{Var}(X_{\mathrm{inv}}^{(d)})$
given in Lemma 2.6.

To compute $\mathrm{Var}(\hat{X}_{\mathrm{inv}}^{(d)}),$ we ignore all
constant parts appearing in $\hat{X}_{\mathrm{inv}}^{(d)}$. By (5a), we can
split $\mathbb{E}(X_{\mathrm{inv}}^{B,(d)}\mid
Z_{k})=\mathbb{E}(X_{\mathrm{inv}}^{B,(d)}\mid
Z_{k})^{+}+\mathbb{E}(X_{\mathrm{inv}}^{B,(d)}\mid
Z_{k})^{-}+\mathbb{E}(X_{\mathrm{inv}}^{B,(d)}\mid Z_{k})^{\circ},$ where

$\displaystyle\mathbb{E}(X_{\mathrm{inv}}^{B,(d)}\mid Z_{k})^{+}$
$\displaystyle:=\sum_{(i,j)\in\mathfrak{N}_{n,d}}\mathbb{P}(Z_{i}>Z_{j}\mid
Z_{k})\,,$ $\displaystyle\mathbb{E}(X_{\mathrm{inv}}^{B,(d)}\mid Z_{k})^{-}$
$\displaystyle:=\sum_{(i,j)\in\widetilde{\mathfrak{N}}_{n,d}}\mathbb{P}(-Z_{i}>Z_{j}\mid
Z_{k})\,,$ $\displaystyle\mathbb{E}(X_{\mathrm{inv}}^{B,(d)}\mid
Z_{k})^{\circ}$ $\displaystyle:=\sum_{i=1}^{n}\mathbb{P}(Z_{i}<0\mid
Z_{k})\,.$

The third sum $\mathbb{E}(X_{\mathrm{inv}}^{B,(d)}\mid Z_{k})^{\circ}$ is
asymptotically negligible. So, to obtain a representation
$\hat{X}_{\mathrm{inv}}^{B,(d)}=\displaystyle{\sum\nolimits_{i=1}^{n}\omega_{d}(k)Z_{k}}+\text{const},$
we write $\omega_{d}(k)=\omega_{d}(k)^{+}+\omega_{d}(k)^{-}$, with
$\omega_{d}(k)^{+}$ stemming from $\mathbb{E}(X_{\mathrm{inv}}^{B,(d)}\mid
Z_{k})^{+}$ and $\omega_{d}(k)^{-}$ stemming from
$\mathbb{E}(X_{\mathrm{inv}}^{B,(d)}\mid Z_{k})^{-},$ respectively.

For $\omega_{d}(k)^{+},$ we can use the previous counting method. However, we
have to take into account that on $B_{n}$,

$\displaystyle\mathbb{P}(Z_{k}>Z_{j}\mid Z_{k})$
$\displaystyle=\frac{Z_{k}+1}{2}\,,$ $\displaystyle\mathbb{P}(Z_{i}>Z_{k}\mid
Z_{k})$ $\displaystyle=\frac{1-Z_{k}}{2}\,,$
$\displaystyle\mathbb{P}(-Z_{k}>Z_{j}\mid Z_{k})$
$\displaystyle=\frac{1-Z_{k}}{2}\,,$ $\displaystyle\mathbb{P}(Z_{i}>Z_{k}\mid
Z_{k})$ $\displaystyle=\frac{1-Z_{k}}{2}\,.$

In conclusion, the coefficients $\omega_{d}(k)^{+}$ on $B_{n}$ are half of the
coefficients $\omega_{d}(k)$ on $S_{n}$. Note that
$\omega_{d}(k)^{+}=(n-2k+1)/2$ for all $d\geq n$. Moreover,
$\omega_{d}(k)^{-}=-\widetilde{N}_{n,d}^{(k)}/2,$ where

$\widetilde{N}_{n,d}^{(k)}:=|\\{(i,j)\in\widetilde{\mathfrak{N}}_{n,d}\mid\text{either}~{}i=k~{}\text{or}~{}j=k\\}|\,.$

Recall that by Remark 2.5, the short case is now $d<n$ and the long case is
$d\geq n$. By counting the number of such $(k,j)$ and $(i,k)$, we can
determine the asymptotic quantification of
$\mathrm{Var}\left(\hat{X}_{\mathrm{inv}}^{B,(d)}\right)$.

###### Proof of Lemma 3.3 for the groups $B_{n}$ and $D_{n}$.

If $d\leq n/2$, then all pairs in $\widetilde{\mathfrak{N}}_{n,d}$ are located
within $K_{1}$, yielding

$\displaystyle\omega_{d}(k)$
$\displaystyle=\begin{cases}(d-k+1)/2-(d-k-1)/2,&k\leq d/2\\\
(d-k+1)/2-(d-k)/2,&d/2<k\leq d\\\ 0,&d<k\leq n-d\\\ (n-d-k)/2,&n-d<k\leq
n\end{cases}$ $\displaystyle=\begin{cases}(n-d-k)/2,&n-d<k\leq n\\\
O(1),&\text{otherwise}\end{cases}\,.$

In conclusion, if $d\leq n/2,$ then

$\displaystyle\mathrm{Var}\left(\hat{X}_{\mathrm{inv}}^{B,(d)}\right)$
$\displaystyle=\frac{1}{3}\sum_{k=1}^{n}\omega_{d}(k)^{2}+O(d^{2})=\frac{1}{12}\sum_{k=n-d+1}^{n}(n-d-k)^{2}+O(d^{2})$
$\displaystyle=\frac{1}{12}\frac{d(d+1)(2d+1)}{6}+O(d^{2})=\frac{1}{36}d^{3}+O(d^{2}).$

Due to
$\mathrm{Var}\left(X_{\mathrm{inv}}^{B,(d)}\right)=d^{3}/36+nd/12+O(d^{2})$
according to Lemma 2.6, we obtain the same condition as in Lemma 3.3, namely,
$d\gg\sqrt{n}$.

If $n/2<d<n$, then the pairs in $\widetilde{\mathfrak{N}}_{n,d}$ also cover
$K_{2}$. For $k\in K_{2},$ there cannot be any pairs $(k,j)$ if $n/2<d\leq
2n/3,$ while this is possible if $d>2n/3$. However, the difference between
these two subcases is only marginal. If $n/2<d\leq 2n/3,$ we obtain

$\omega_{d}(k)=\begin{cases}O(1),&k\in K_{1}\\\
\bigl{(}n-2k+1-(d-k)\bigr{)}/2,&k\in K_{2}\\\ (n-d-k)/2,&k\in
K_{3}\end{cases}\,.$

If $d>2n/3,$ then

$\omega_{d}(k)=\begin{cases}O(1),&k\in K_{1}\\\
\bigl{(}n-2k+1-(d-k-1)\bigr{)}/2,&k\in K_{2},k\leq d/2\\\
\bigl{(}n-2k+1-(d-k)\bigr{)}/2,&k\in K_{2},k>d/2\\\ (n-d-k)/2,&k\in
K_{3}\end{cases}\,.$

In conclusion, if $n/2<d<n,$ then

$\displaystyle\mathrm{Var}\left(\hat{X}_{\mathrm{inv}}^{B,(d)}\right)$
$\displaystyle=\frac{1}{3}\sum_{k=1}^{n}\omega_{d}(k)^{2}+O(d^{2})$
$\displaystyle=\frac{1}{12}\sum_{k=n-d+1}^{d}(n-d-k+1)^{2}+\frac{1}{12}\sum_{k=d+1}^{n}(n-d-k)^{2}+O(d^{2})$
$\displaystyle=\frac{1}{12}\sum_{k=n-d+1}^{d}(n-d-k+1)^{2}+\frac{1}{12}\sum_{k=1}^{2d-n}k^{2}+O(d^{2})$
$\displaystyle=\frac{1}{12}\sum_{k=1}^{2d-n-1}k^{2}+\frac{1}{12}\left(\sum_{k=1}^{d}k^{2}-\sum_{k=1}^{2d-n}k^{2}\right)+O(d^{2})$
$\displaystyle=\frac{1}{36}d^{3}+O(d^{2})\,.\qed$

###### Proof of Lemma 3.4 for the groups $B_{n}$ and $D_{n}$.

In the long case, the main focus is on counting $\widetilde{N}_{n,d}^{(k)}$.
Figure 2 illustrates the positions of pairs
$(i,j)\in\widetilde{\mathfrak{N}}_{n,d}$ in the long case for the exemplary
choice of $n=12,d=16$.

Figure 2. Visualization of pairs $(i,j)$ in $\widetilde{\mathfrak{N}}_{n,d}$
in the long case for $n=12$ and $d=16$. The numbers $d-n$ and $d/2$ are
important case distinction thresholds for counting the pairs $(k,j)$ and
$(i,k),$ respectively.

With help of Figure 2, it is straightforward to count

$\widetilde{N}_{n,d}^{(k)}=\begin{cases}n-1,&1\leq k\leq d-n\\\
d-k-1,&d-n<k\leq d/2\\\ d-k,&d/2<k\leq n\end{cases}\,.$

This result is also illustrated in Figure 3, which displays the number of
pairs $(i,k)$ and the number of pairs $(k,j)$.

Figure 3. Plots of the numbers of pairs $(i,k)$ (red) and $(k,j)$ (blue) in
$\widetilde{\mathfrak{N}}_{n,d}$. The sum of these two numbers is
$\widetilde{N}_{n,d}^{(k)},$ which is displayed by the black crossed points.

Therefore, in the long case, we have

$\omega_{d}(k)=\begin{cases}1-k,&1\leq k\leq d-n\\\ (n+2-d-k)/2,&d-n<k\leq
d/2,\\\ (n+1-d-k)/2,&d/2<k\leq n.\end{cases}$

We compute

$\displaystyle\sum_{k=1}^{n}\omega_{d}(k)^{2}$
$\displaystyle=\sum_{k=1}^{d-n}(k-1)^{2}+\frac{1}{4}\sum_{k=d-n+1}^{n}(n-d+1-k)^{2}+O(n^{2})$
$\displaystyle=\sum_{k=1}^{d-n-1}k^{2}+\frac{1}{4}\sum_{k=d-n+1}^{n}\bigl{(}k-(d-n-1)\bigr{)}^{2}+O(n^{2})$
$\displaystyle=\frac{1}{3}(d-n-1)^{3}+\frac{1}{12}(2n-d-1)^{3}+O(n^{2})$
$\displaystyle=\frac{d^{3}}{4}-\frac{d^{2}n}{2}-\frac{3d^{2}}{4}+dn+\frac{3d}{4}+\frac{n^{3}}{3}-\frac{n}{2}+\frac{1}{4}+O(n^{2})\,.$

Due to $\mathrm{Var}(Z_{k})=1/3,$ we obtain

$\mathrm{Var}(\hat{X}_{\mathrm{inv}}^{B,(d)})=\frac{d^{3}}{12}-\frac{d^{2}n}{6}+\frac{n^{3}}{9}+O(d^{2}).$

On the other hand, by Lemma 2.6,

$\mathrm{Var}(X_{\mathrm{inv}}^{B,(d)})=-\frac{1}{12}d^{3}+\frac{1}{3}nd^{2}-\frac{1}{3}n^{2}d+\frac{1}{9}n^{3}+O(d^{2})\,,$

i.e., the leading terms $d^{3}/12-d^{2}n/6$ and $-d^{3}/12+nd^{2}/3-n^{2}d/3$
need to match. Writing $d=cn(1+o(1))$ for $1\leq c\leq 2,$ this gives the
condition

$\frac{(cn)^{3}}{12}-\frac{c^{2}n^{3}}{6}=-\frac{(cn)^{3}}{12}+\frac{c^{2}n^{3}}{3}-\frac{cn^{3}}{3}\Longleftrightarrow\frac{c^{3}}{6}-\frac{c^{2}}{2}+\frac{c}{3}=0\,,$

which holds precisely for $c=1$ and $c=2$. ∎

## References

  * [1] Anders Björner and Francesco Brenti, _Combinatorics of Coxeter groups_ , vol. 231, Springer Science & Business Media, 2006.
  * [2] Miklós Bóna, _Generalized descents and normality_ , Electronic Journal of Combinatorics 15 (2008), N21.
  * [3] Jinyuan Chang, Xiaohui Chen, and Mingcong Wu, _Central limit theorems for high dimensional dependent data_ , arXiv preprint arXiv:2104.12929 (2021).
  * [4] Harold S. M. Coxeter, _The complete enumeration of finite groups of the form $r_{i}^{2}=(r_{i}r_{j})^{k_{ij}}=1$_, Journal of the London Mathematical Society 1 (1935), no. 1, 21–25.
  * [5] Filippo de Mari and Mark A. Shayman, _Generalized eulerian numbers and the topology of the hessenberg variety of a matrix_ , Acta Applicandae Mathematica 12 (1988), 213–235.
  * [6] Philip Dörr and Johannes Heiny, _Extremes of joint inversions and descents on finite coxeter groups_ , arXiv preprint arXiv:2309.17314 (2023).
  * [7] Philip Dörr and Thomas Kahle, _Extreme values of permutation statistics_ , arXiv preprint arXiv:2205.01426 (2022).
  * [8] Thomas Kahle and Christian Stump, _Counting inversions and descents of random elements in finite Coxeter groups_ , Mathematics of Computation 89 (2020), no. 321, 437–464.
  * [9] Kathrin Meier and Christian Stump, _Central limit theorems for generalized descents and generalized inversions in finite root systems_ , preprint, arXiv:2202.05580 (2022).
  * [10] John Pike, _Convergence rates for generalized descents_ , the electronic journal of combinatorics (2011), P236–P236.
  * [11] Aad W. van der Vaart, _Asymptotic statistics_ , vol. 3, Cambridge university press, 2000.


# An entanglement-aware quantum computer simulation algorithm

Maxime Oliva<EMAIL_ADDRESS>Eviden Quantum Laboratory, Les Clayes-
sous-Bois, France

###### Abstract

The advent of quantum computers promises exponential speed ups in the
execution of various computational tasks. While their capabilities are
hindered by quantum decoherence, they can be exactly simulated on classical
hardware at the cost of an exponential scaling in terms of number of qubits.
To circumvent this, quantum states can be represented as matrix product states
[1, 2] (MPS), a product of tensors separated by so-called bond dimensions.
Limiting bond dimensions growth approximates the state, but also limits its
ability to represent entanglement. Methods based on this representation have
been the most popular tool at simulating large quantum systems. But how to
trust resulting approximate quantum states for such intractable systems sizes
? I propose here a method for inferring the fidelity of an approximate quantum
state without direct comparison to its exact counterpart, and use it to design
an “entanglement-aware” (EA) algorithm for both pure and mixed states. As
opposed to state of the art methods which limit bond dimensions up to an
arbitrary maximum value, this algorithm receives as input a fidelity, and
adapts dynamically its bond dimensions to both local entanglement and noise
such that the final quantum state fidelity at least reaches the input
fidelity. I show that this algorithm far surpasses standard fixed bond
dimension truncation schemes. In particular, a noiseless random circuit of 300
qubits and depth 75 simulated using MPS methods takes one week of computation
time, while EA-MPS only needs 2 hours to reach similar quantum state fidelity.

##  I Introduction

Tensor network methods have been among the most popular avenue at
circumventing the exponential scaling of exact quantum simulations [1, 2, 3,
4, 5]. While lowly entangled pure states can be efficiently simulated as
matrix product states[1, 2], mixed states can be simulated as matrix product
operators[2, 4]. Both approaches have allowed the reach of systems sizes far
beyond what exact computation can achieve.

Any pure state can be represented as a matrix product state (MPS). The state
vector of a quantum state of $N$ qubits is cast into a factorized form of $N$
tensors connected to each other with what are generally called “bond
dimensions”. Physically, bond dimensions can be thought of as the amount of
entanglement a quantum state can encapsulate[1]. Contracting every tensor
along its bond dimensions gives back the quantum state in its standard state
vector form. The usual approach to approximating MPS consists in limiting
these bond dimensions up to a maximum value set arbitrarily prior to the
simulation. The more entanglement, the bigger the truncation errors, which
limits MPS effectiveness to lowly entangled, but possibly very large quantum
states. Such methods suffer two limitations. First, there is no way to assess
how well the resulting quantum state approximates the exact quantum state when
the system size is out of exact simulation reach. Secondly, there is no
general method for guessing which maximum bond dimension is best for a given
problem.

I show that for both pure and mixed states simulations, an approximate quantum
state fidelity with respect to its exact counterpart can be indirectly
computed, that is without ever having to compute the exact quantum state
itself. This computation is done in real-time throughout the simulation, by
efficiently computing the fidelity of every bond dimension truncation.

Based on this result, an entanglement-aware simulation algorithm can be
designed, fully leveraging the presence or absence of both noise and
entanglement. Rather than limiting bond dimensions uniformly as is the case in
standard algorithms, bond dimensions are allowed to increase or decrease
solely based on how this affects the resulting quantum state fidelity. As such
the algorithm does not receive as input a maximum bond dimension, but instead
a fidelity the final quantum state has to at least reach. Decreasing (resp.
increasing) the desired fidelity decides how aggressive (resp. conservative)
the algorithm will be at truncating bond dimensions.

Since bond dimensions are kept as low as possible for every operation, the
algorithm efficiency far surpasses that of regular fixed maximum bond
dimensions algorithms, while simultaneously allowing trustworthy output
quantum states.

##  II The matrix product formalism

Any pure quantum state can be written as a state vector $\ket{\psi}$ defined
by :

$\displaystyle\ket{\psi}=\sum_{\sigma}C_{\sigma_{1},...,\sigma_{N}}\ket{\sigma_{1}}\otimes...\otimes\ket{\sigma_{N}}$
(1)

with $C_{\sigma_{1},...,\sigma_{N}}$ a 1D tensor containing $2^{N}$ complex
values and $\\{\ket{\sigma_{i}}\\}$ forming an orthonormal basis.

As described in Vidal’s article[1], $\ket{\psi}$ can be decomposed as an MPS
via successive singular value decomposition (SVD) of $\ket{\psi}$ in Eq. (1):

$\allowbreak\ket{\psi}=\sum_{\begin{subarray}{c}\sigma,\chi\end{subarray}}A^{[1]\sigma_{1}}_{1,\chi_{1}}A^{[2]\sigma_{2}}_{\chi_{1},\chi_{2}}...A^{[N]\sigma_{N}}_{\chi_{N-1},1}\
\ket{\sigma_{1}}\otimes...\otimes\ket{\sigma_{N}}$ (2)

We obtain a product of $N$ complex-valued tensors $\\{A^{[i]\sigma_{i}}\\}$,
separated by bond dimensions $\\{\chi_{1},\chi_{2},...,\chi_{N-1}\\}$.
Assuming the $\chi$ are bounded by a maximum bond dimension $\chi_{max}$, the
number of values contained in the MPS in Eq. (2) scales in
$\mathcal{O}(N\chi_{max}^{2})$ with $N$ the number of qubits of the system.

Similarly, for mixed state, instead of considering the state vector
$\ket{\psi}$, we consider the density matrix $\hat{\rho}$:

$\allowbreak\hat{\rho}=\sum_{\begin{subarray}{c}\sigma,\sigma^{\prime}\end{subarray}}C_{\sigma_{1}\sigma^{\prime}_{1},...,\sigma_{N}\sigma^{\prime}_{N}}|{\sigma_{1}}\rangle\langle{\sigma^{\prime}_{1}}|\otimes...\otimes|{\sigma_{N}}\rangle\langle{\sigma^{\prime}_{N}}|$
(3)

with $C_{\sigma_{1}\sigma^{\prime}_{1},...,\sigma_{N}\sigma^{\prime}_{N}}$ a
2D tensor containing $2^{N}\times 2^{N}$ complex values and
$\\{\ket{\sigma_{i}}\\}$ forming an orthonormal basis. Again, $\hat{\rho}$ can
be cast into a matrix product operator (MPO) [2] via successive SVD
decompositions:

$\displaystyle\allowbreak\hat{\rho}=\sum_{\chi_{1}...\chi_{N}}A^{[1]\sigma_{1},\sigma^{\prime}_{1}}_{1,\chi_{1}}$
$\displaystyle
A^{[2]\sigma_{2},\sigma^{\prime}_{2}}_{\chi_{1},\chi_{2}}...A^{[N]\sigma_{N},\sigma^{\prime}_{N}}_{\chi_{N-1},1}$
$\displaystyle|{\sigma_{1}}\rangle\langle{\sigma^{\prime}_{1}}|\otimes...\otimes|{\sigma_{N}}\rangle\langle{\sigma^{\prime}_{N}}|$
(4)

Limiting the growth of bond dimensions $\chi\leq\chi_{max}$ approximates the
state, and let the MPO scale in $\mathcal{O}(N^{2}\chi_{max}^{3})$.

When it comes to quantum circuit simulations, such representations can only
simulate circuits with linear nearest neighbour (LNN) topology, meaning gates
have to be applied on neighbouring qubits only. Bond dimensions grow only by
application of operators acting on multiple qubits.

##  III Canonicalization

Casting quantum state vectors or density matrices into their matrix product
representations allow for various simplifications. Most of these
simplifications use the concept of canonicalization[2, 6], which is a
consequence of the SVD or QR operations required for constructing and
truncating quantum states in matrix product representations. Take the SVD of a
matrix $M$, we have $M=U\Lambda V^{\dagger}$, where $\Lambda$ is diagonal,
$U^{\dagger}U=I$ and $(V^{\dagger}V)^{\dagger}=VV^{\dagger}=I$. As such, $U$
is left-normalized, while $V^{\dagger}$ is right-normalized. Similarly, a QR
operation on a matrix $M$ gives $M=QR^{\dagger}$ with $Q$ and $R^{\dagger}$
respectively left and right-normalized matrices. By controlling which tensors
are left- or right-normalized, many operations can be done on a small subset
of tensors, rather than on the entire quantum state.

To illustrate this, let us compute the expectation value of an observable
$\hat{O}$ whose support is on $\ket{\sigma_{k}}$ and $\ket{\sigma_{k+1}}$ for
a quantum state $\ket{\psi}$ in the MPS form defined in Eq. (2):

$\displaystyle\langle\hat{O}\rangle_{\psi}=$
$\displaystyle\sum_{\begin{subarray}{c}\sigma,\sigma^{\prime},\chi\end{subarray}}A^{[0]\sigma^{\prime}_{0}}_{1,\chi_{1}}...A^{[k]\sigma^{\prime}_{k}}_{\chi_{k},\chi_{k+1}}A^{[k+1]\sigma^{\prime}_{k+1}}_{\chi_{k+1},\chi_{k+2}}...$
$\displaystyle
A^{[N]\sigma^{\prime}_{N}}_{\chi_{N-1},1}O^{\sigma^{\prime}_{k},\sigma^{\prime}_{k+1}}_{\sigma_{k},\sigma_{k+1}}A^{[0]\sigma_{0}}_{1,\chi_{1}}...$
$\displaystyle
A^{[k]\sigma_{k}}_{\chi_{k},\chi_{k+1}}A^{[k+1]\sigma_{k+1}}_{\chi_{k+1},\chi_{k+2}}...A^{[N]\sigma_{N}}_{\chi_{N-1},1}\langle\boldmath{\sigma^{\prime}}|\boldmath{\sigma}\rangle$
(5)

For now the MPS $\ket{\psi}$ is in arbitrary form. We apply QR to each tensor
$A^{[i]}$ for $i\in\ [1,k[$ such that $\sum A^{\dagger[i]}A^{[i]}=I$, and
similarly apply a RQ for $i\in\ [N,k+1[$ such that $\sum
A^{[i]}A^{\dagger[i]}=I$. The state is now in “canonical form”. The left-
normalized tensors multiplied with themselves simplify to identity, and
similarly for the right-normalized tensors.

$\hat{O}$ expectation value can then be computed by only multiplying:

$\displaystyle\allowbreak\langle\hat{O}\rangle_{\psi}=\sum_{\begin{subarray}{c}\sigma,\sigma^{\prime},\chi\end{subarray}}$
$\displaystyle
A^{[k]\sigma^{\prime}_{k}}_{\chi_{k},\chi_{k+1}}A^{[k+1]\sigma^{\prime}_{k+1}}_{\chi_{k+1},\chi_{k+2}}$
$\displaystyle
O^{\sigma^{\prime}_{k},\sigma^{\prime}_{k+1}}_{\sigma_{k},\sigma_{k+1}}A^{[k]\sigma_{k}}_{\chi_{k},\chi_{k+1}}A^{[k+1]\sigma_{k+1}}_{\chi_{k+1},\chi_{k+2}}$
(6)

The canonical form is needed for efficiently computing the truncation
fidelities, and is enforced throughout the entanglement-aware simulation.

##  IV Quantum state fidelity

The quantum state fidelity is a measure of the “closeness” of two quantum
states.

For pure states, the fidelity $\mathcal{F}$ between two states $\ket{\psi}$
and $\ket{\phi}$ is defined by the square modulus of their overlap:

$\displaystyle\mathcal{F}(\psi,\phi)=|\langle\psi|\phi\rangle|^{2}$ (7)

Similarly, Josza[7] defined the quantum state fidelity between two mixed
states $\rho$ and $\sigma$ as:

$\displaystyle\mathcal{F}(\rho,\sigma)=Tr\left(\sqrt{\sqrt{\rho}\sigma\sqrt{\rho}}\right)^{2}$
(8)

However, truncations of mixed states in MPO representation can lead to the
loss of their positive semi-definite (PSD) property, which is a necessary
condition for computing the matrix square roots in Eq. (8).

An alternative formulation obeying all four of Josza’s axioms[7] has been
given by X. Wang et al. [8] for mixed states. It is defined by:

$\displaystyle\mathcal{F}(\rho,\sigma)=\frac{|Tr(\rho\sigma)|}{\sqrt{Tr(\rho^{2})Tr(\sigma^{2})}}$
(9)

This definition removes the need for the PSD property of the quantum state,
and is usable for non-normalized density matrices.

For the rest of the article, the quantum state fidelity for pure state will
refer to Eq. (7), while it will refer to Eq. (9) for mixed states.

##  V Truncation fidelity

We call truncation fidelity the quantum state fidelity before and after
truncation. We show truncation fidelities can be efficiently computed using
only the singular values associated with the bond dimension that is to be
truncated.

Let us assume we want to truncate a pure state $\ket{\psi}$ at site $T$, to
obtain the truncated state $|\tilde{\psi}\rangle$. The state $\ket{\psi}$ is
in canonical form. Contracting both tensors neighbouring the bond dimension
$T$, and performing its SVD, we obtain the diagonal matrix $\Lambda$,
containing the singular values in decreasing order. Truncating this matrix
reduces the overall bond dimension at site $T$, and approximates the state.

Since the state is in canonical form, the truncation fidelity of such
truncation is given by:

$\displaystyle
f_{T}(\tilde{\Lambda})=|\langle\psi|\tilde{\psi}\rangle|^{2}=\sum_{i}\Lambda_{ii}^{2}\tilde{\Lambda}_{ii}^{2}$
(10)

with $\Lambda$ and $\tilde{\Lambda}$ the diagonal matrices representing the
singular values obtained at site $T$ _before_ and _after_ truncation
respectively.

Similarly, truncating a mixed state $\rho$ in canonical form at site $T$, we
obtain the truncated mixed state $\tilde{\rho}$. Using Eq. (9), and
simplifying the tensors by using the canonical simplifications previously
described, we obtain:

$\displaystyle
f_{T}(\tilde{\Lambda})=\frac{\sum_{i}\Lambda_{ii}^{2}\tilde{\Lambda}_{ii}^{2}}{\sum_{i}\Lambda_{ii}^{2}\sum_{i}\tilde{\Lambda}_{ii}^{2}}$
(11)

In both cases, the truncation fidelity $f_{T}$ depends only on the actual
truncated singular values in $\tilde{\Lambda}$. Computing the truncation
fidelity is now simply a scalar operation between singular values _before_ and
_after_ truncation.

##  VI Entanglement-aware simulations

Throughout a matrix product-based simulation algorithm, the quantum state
endures successive truncations. The more entangled the state, the bigger the
truncation errors, and thus the higher the bond dimensions need to be. As
shown in Zhou et. al. article[9] and further confirmed in [10], the overall
quantum fidelity $\mathcal{F}$ between an exact quantum state $\ket{\psi}$ and
its truncated counterpart $|\tilde{\psi}\rangle$ can be approximated for
noiseless simulations by:

$\displaystyle\mathcal{F}(n)\approx\prod_{i=1}^{n}f_{i}$ (12)

with $n$ the number of truncations and $f_{i}$ the individual truncation
fidelities obtained at each truncation. This approximation proves remarkably
robust at all regimes studied[9, 10]. In the case of noisy systems, Eq. (12)
underestimates $\mathcal{F}$ since the additional noise reduces the truncation
errors committed. In the general case we have:

$\displaystyle\mathcal{F}(n)\gtrsim\prod_{i=1}^{n}f_{i}$ (13)

By enforcing the canonical form at all times throughout the simulation, we can
efficiently compute truncation fidelities, and thus have direct access to how
close the resulting approximated quantum state is to its exact non-truncated
counterpart. Using Eq. (10) for pure states or Eq. (11) for mixed states, we
know how many singular values have to be truncated in order to reach a
specific target truncation fidelity. Since we know in advance how many 2-qubit
gates are to be applied onto the quantum state, we can revert Eq. (12) and
define the target truncation fidelities $\\{f^{(i)}_{t}\\}$ to reach at every
truncation so that a fidelity $\mathcal{F}_{min}$ is at least obtained. The
bond dimension truncation scheme becomes adaptive.

The simplest way to define the target truncation fidelities
$\\{f^{(i)}_{t}\\}$ is to define them uniformly based on a fidelity
$\mathcal{F}_{min}$ the quantum state has to at least reach:

$\displaystyle f_{t}^{(i)}(n)=\mathcal{F}_{min}^{1/n}\quad\forall i\ \in\
[0,n]$ (14)

Figure 1: Comparison between the target fidelity and the actual quantum state
fidelity obtained for a 25 qubits Haar-random circuit of depth 20. The “naive”
strategy corresponds to no target truncation fidelity update throughout the
simulation. The “nearest” strategy corresponds to the update of only the
upcoming target truncation fidelity, while the “global” strategy updates all
upcoming target truncation fidelities.

Here, $n$ is the number of 2-qubit gates of the quantum circuit to simulate,
i.e. the number of possible truncations. Throughout the simulation, at every
possible truncation, we truncate bond dimensions such that Eq. (10) for pure
states and Eq. (11) for mixed states are as close to but not lower than
$f_{t}$. In practice the actual truncation fidelity will always be higher than
the target truncation fidelity $f_{t}$, resulting in a final quantum state
fidelity higher than the chosen input fidelity $\mathcal{F}_{min}$. To achieve
a quantum state fidelity closer to the initial fidelity desired, one can
update target truncation fidelities dynamically after every truncation so that
the product of truncation fidelities is as close to $\mathcal{F}_{min}$ as
possible, see Fig. 1.

Since bond dimensions now depend on a direct truncation fidelity measure, they
increase or decrease based on how they affect the final quantum state
fidelity. The smaller the $\mathcal{F}_{min}$ value, the lower the target
truncation fidelities, and the more aggressive the truncations will be. As
such, the algorithm is “entanglement-aware” since it is able to detect local
entanglement changes, and adapt bond dimensions accordingly.

Fig. 2 demonstrates this adaptivity on a circuit composed of two subcircuits,
the second subcircuit being the adjoint of the first. The initial state is a
product state requiring only bond dimensions of 1. In the first subcircuit, as
entanglement increases bond dimensions increase, to then decrease in the
second subcircuit. At the end of the simulation, the state goes back to the
initial state with a bond dimensions of 1 everywhere on the tensor.

In Fig. 2, as entanglement dominates, bond dimensions increase. But eventually
noise takes over and reduces the overall state entanglement, leading to
decreased bond dimension needs. The lower the desired fidelity, the faster the
algorithm is at picking up these changes.

Figure 2: (_a_) Evolution of the bond dimensions at each MPS site throughout
the simulation of a 40 qubits circuit composed of a Haar-random circuit of 20
layers followed by its adjoint. All quantum gates are defined by Haar-random
unitary matrices. The bond dimensions increase as the entanglement increases.
When the adjoint circuit is reached, entanglement entropy decreases, and the
bond dimensions follow. (_b_) Simulation of a noisy quantum circuit of 8
qubits with depolarizing noise parameters $\epsilon_{1}=\epsilon_{2}=0.05$ for
both one- and two-qubit gates. Noise increases as the depth of the circuit
increases. This leads to reduced entanglement needs, which is then exploited
by the entanglement-aware algorithm to reduce the state bond dimensions. The
lower the desired fidelity $\mathcal{F}_{min}$, the more aggressive the
algorithm will be at truncating bond dimensions. Inversely, the higher
$\mathcal{F}_{min}$, the more conservative the truncations will be.

This brings new insights on the dynamics of quantum systems. One could for
example expect that the bond dimensions requirements in a Haar-random circuit
would be maximal around the middle junction of the MPS. But this is not always
the case, as shown in Fig. 3. Noisy systems tend to display large bond
dimensions variations that cannot be exploited using fixed maximum bond
dimension truncation schemes, see Fig. 3. This makes adaptive truncation
schemes particularly effective for noisy systems.

Figure 3: (_a_) Simulation of a 40 qubits Haar-random quantum circuit of depth
20 with a target fidelity of $\mathcal{F}_{min}=0.9$. Contrary to intuition,
bond dimensions are not higher at the center of the MPS. (_b_) Noisy
simulation of a 30 qubits Haar-random quantum circuit of depth 10 with a
target fidelity of $\mathcal{F}_{min}=0.9$. Large local entanglement
variations are exploited and induce large bond dimensions variations along the
MPO.

Figure 4:  (_a_) Simulation time of MPS and EA-MPS for a noiseless random
circuit of depth 75. (_b_) Simulation time of MPO and EA-MPO for a noisy
random circuit of depth 25. For both figure and for each point, the MPS (resp.
MPO) maximum bond dimension is taken as the maximum bond dimension encountered
during the EA-MPS (resp. EA-MPO) simulation. This ensures similar final
quantum state fidelity $\mathcal{F}_{min}$. For example, for 300 qubits we
obtain $\mathcal{F}_{min}=0.999$ for both regular MPS and MPO, while we obtain
$\mathcal{F}_{min}=0.998$ for EA-MPS and EA-MPO. (_c_) and (_d_) show the
final state bond dimensions for the 300 qubits data point in (_a_) and (_b_)
respectively. The blue lines show the bond dimensions of the state obtained
through EA simulation, while the red dotted lines are the bond dimensions
obtained through regular fixed bond dimension simulation.

Benchmarks of noiseless and noisy entanglement-aware simulations are presented
in Fig. 4 and Fig. 4, and compared to standard methods. The random circuit
used is defined based on Cheng’s article[5], which alternates randomly chosen
layers of one- and two-qubit gates. The standard MPS simulation requires one
week of simulation time for a 300 qubits random circuit of depth 75, while the
entanglement-aware version only needs 2 hours. For a depth of 25, the standard
MPO simulation requires also one week, but EA-MPO only needs 40 minutes. Fig.
4 and Fig. 4 show the bond dimensions of the resulting approximate quantum
state for 300 qubits. One has to keep in mind that, assuming it is possible to
guess which maximum bond dimension was needed in the first place, standard MPS
and MPO simulation algorithms have to simulate the state using very large bond
dimensions over the entire tensor network for a final quantum state fidelity
gain lower than $0.001$. Note also that while it may be tempting to truncate
some of the high peaks displayed in Fig. 4 or Fig. 4, it is exactly these
truncations which would lead to the largest fidelity losses on the final
quantum state, since it is where the entanglement is highest, and thus where
the truncation errors are maximal.

Both EA-MPS and EA-MPO simulation algorithms vastly outperform regular MPS and
MPO simulations in all cases studied, may it be in terms of quantum state
fidelity, computation time or memory footprint.

It is however difficult to assess the exact scaling improvement over regular
fixed maximum bond dimension methods, as it is problem dependent. Also, as
bond dimensions are allowed to grow indefinitely, a quantum circuit inducing
too much entanglement will remain almost impossible to simulate. For that
reason, it is advisable to cap adaptive bond dimensions up to a maximum value
which, if reached, enforces bond dimension truncations. In that case, the
desired fidelity can no longer be guaranteed.

##  VII Conclusion

In conclusion, I have shown that an adaptive bond dimension method far
outperforms methods based on arbitrarily chosen maximum bond dimensions. The
metric used in this article is the quantum state fidelity because of its
straightforward relationship with truncation fidelities, but other metrics
could be used. This method only adds the computation of the truncation
fidelities, but this additional cost is largely compensated for by the overall
bond dimension reduction it allows. Since bond dimensions growth can still be
constrained by a maximum bond dimension, and because the loss in final quantum
state fidelity is negligible, I expect adaptive truncation methods to push
even further the quantum systems simulation capabilities of tensor network-
based methods on classical hardware.

## References

  * Vidal [2003] G. Vidal, Efficient classical simulation of slightly entangled quantum computations, Physical Review Letters 91, 147902 (2003), arxiv:quant-ph/0301063 .
  * Schollwoeck [2011] U. Schollwoeck, The density-matrix renormalization group in the age of matrix product states, Annals of Physics 326, 96 (2011), arxiv:1008.3477 .
  * White  _et al._ [2018] C. D. White, M. Zaletel, R. S. K. Mong, and G. Refael, Quantum dynamics of thermalizing systems, Physical Review B 97, 035127 (2018), arxiv:1707.01506 .
  * Noh  _et al._ [2020] K. Noh, L. Jiang, and B. Fefferman, Efficient classical simulation of noisy random quantum circuits in one dimension, Quantum 4, 318 (2020), arxiv:2003.13163 .
  * Cheng  _et al._ [2021] S. Cheng, C. Cao, C. Zhang, Y. Liu, S.-Y. Hou, P. Xu, and B. Zeng, Simulating noisy quantum circuits with matrix product density operators, Physical Review Research 3, 023005 (2021).
  * Orus [2014] R. Orus, A Practical Introduction to Tensor Networks: Matrix Product States and Projected Entangled Pair States, Annals of Physics 349, 117 (2014), arxiv:1306.2164 .
  * Jozsa [1994] R. Jozsa, Fidelity for Mixed Quantum States, Journal of Modern Optics 41, 2315 (1994).
  * Wang  _et al._ [2008] X. Wang, C.-S. Yu, and X. X. Yi, An alternative quantum fidelity for mixed states of qudits, Physics Letters A 373, 58 (2008), arxiv:0807.1781 .
  * Zhou  _et al._ [2020] Y. Zhou, E. M. Stoudenmire, and X. Waintal, What limits the simulation of quantum computers?, Physical Review X 10, 041038 (2020), arxiv:2002.07730 .
  * Ayral  _et al._ [2022] T. Ayral, T. Louvet, Y. Zhou, C. Lambert, E. M. Stoudenmire, and X. Waintal, A density-matrix renormalization group algorithm for simulating quantum circuits with a finite fidelity (2022), arxiv:2207.05612 [cond-mat, physics:quant-ph] .