#  Complete Asymptotic Expansions and the High-Dimensional Bingham
Distributions

Armine Bagyan and Donald Richards Department of Statistics, Pennsylvania State
University, University Park, PA 16802, U.S.A. E-mail address:
<EMAIL_ADDRESS>of Statistics, Pennsylvania State University,
University Park, PA 16802, U.S.A. E-mail address<EMAIL_ADDRESS>†Corresponding author.

###### Abstract

For $\displaystyle d\geq 2$, let $\displaystyle X$ be a random vector having a
Bingham distribution on $\displaystyle\mathcal{S}^{d-1}$, the unit sphere
centered at the origin in $\displaystyle\mathbb{R}^{d}$, and let
$\displaystyle\Sigma$ denote the symmetric matrix parameter of the
distribution. Let $\displaystyle\Psi(\Sigma)$ be the normalizing constant of
the distribution and let $\displaystyle\nabla\Psi_{d}(\Sigma)$ be the matrix
of first-order partial derivatives of $\displaystyle\Psi(\Sigma)$ with respect
to the entries of $\displaystyle\Sigma$. We derive complete asymptotic
expansions for $\displaystyle\Psi(\Sigma)$ and
$\displaystyle\nabla\Psi_{d}(\Sigma)$, as $\displaystyle d\to\infty$; these
expansions are obtained subject to the growth condition that
$\displaystyle\|\Sigma\|$, the Frobenius norm of $\displaystyle\Sigma$,
satisfies $\displaystyle\|\Sigma\|\leq\gamma_{0}d^{r/2}$ for all
$\displaystyle d$, where $\displaystyle\gamma_{0}>0$ and $\displaystyle
r\in[0,1)$. Consequently, we obtain for the covariance matrix of
$\displaystyle X$ an asymptotic expansion up to terms of arbitrary degree in
$\displaystyle\Sigma$. Using a range of values of $\displaystyle d$ that have
appeared in a variety of applications of high-dimensional spherical data
analysis we tabulate the bounds on the remainder terms in the expansions of
$\displaystyle\Psi(\Sigma)$ and $\displaystyle\nabla\Psi_{d}(\Sigma)$ and we
demonstrate the rapid convergence of the bounds to zero as $\displaystyle r$
decreases.

Keywords and phrases. Confluent hypergeometric function of matrix argument;
Frobenius norm; Gradient operator; Power sum symmetric function, Zonal
polynomial.

2020 Mathematics Subject Classification. Primary: 60E05, 62H11. Secondary:
62E20, 62R30.

Running head: The High-Dimensional Bingham Distribution.

##  1 Introduction

Let $\displaystyle d$ be a positive integer, where $\displaystyle d\geq 2$.
Denote by $\displaystyle\mathcal{S}^{d-1}$ the unit sphere in $\displaystyle
d$-dimensional Euclidean space, $\displaystyle\mathbb{R}^{d}$, with center at
the origin, and let $\displaystyle\hskip 2.0pt\mathrm{d}x$ be the normalized
surface measure on $\displaystyle\mathcal{S}^{d-1}$. For a symmetric
$\displaystyle d\times d$ matrix $\displaystyle\Sigma=(\sigma_{i,j})$, the
$\displaystyle d\times 1$ random vector $\displaystyle
X=(X_{1},\ldots,X_{d})^{\prime}$ is said to have a Bingham distribution on
$\displaystyle\mathcal{S}^{d-1}$ if the probability density function of
$\displaystyle X$, with respect to the surface measure $\displaystyle\hskip
2.0pt\mathrm{d}x$, is

$\phi(x;{\hskip-3.0pt}\Sigma)=[\Psi_{d}(\Sigma)]^{-1}\exp(x^{\prime}\Sigma
x),$ (1.1)

$\displaystyle x\in\mathcal{S}^{d-1}$, with normalizing constant

$\Psi_{d}(\Sigma)=\int_{\mathcal{S}^{d-1}}\exp(x^{\prime}\Sigma x)\hskip
2.0pt\mathrm{d}x.$ (1.2)

These distributions were introduced by Bingham [4], and they have become
prominent since their appearance. We cite the work of Bagyan [2], Bhattacharya
and Bhattacharya [3], Bingham, Chang, and Richards [5], Chikuse [8], Dryden
[11], Dryden and Mardia [12], Fisher, Lewis, and Embleton [14], Kume and Wood
[19], Mardia and Jupp [20], and Patrangenaru and Ellingson [22] for various
aspects of the theory of these distributions, their applications in random
modulation, spherical data analysis, Procrustes analysis, shape analysis, and
for many references to the literature.

One approach in the study of high-dimensional approximations for the Bingham
distribution is to approximate the distribution of $\displaystyle X$ and then
to deduce approximations for, e.g., $\displaystyle\mathrm{Cov}(X)$, the
covariance matrix of $\displaystyle X$; see, e.g., Chikuse [8], Dryden [11],
and Kume and Wood [19]. In this article, starting from the zonal polynomial
series for $\displaystyle\Psi_{d}(\Sigma)$, we obtain complete asymptotic
expansions, i.e., series expansions in negative powers of $\displaystyle d$,
for $\displaystyle\Psi_{d}(\Sigma)$ and $\displaystyle\mathrm{Cov}(X)$ as
$\displaystyle d\to\infty$. An advantage of our approach is that the resulting
asymptotic expansions are obtained purely in terms of
$\displaystyle\|\Sigma\|=[\hskip 1.0pt\mathrm{tr}\,(\Sigma^{2})]^{1/2}$, the
Frobenius norm of $\displaystyle\Sigma$. Thus we make no assumptions about the
structure of individual matrix entries of $\displaystyle\Sigma$; indeed, our
only constraint on $\displaystyle\Sigma$ is that
$\displaystyle\|\Sigma\|\leq\gamma_{0}d^{r/2}$ for all $\displaystyle d$, with
universal constants $\displaystyle\gamma_{0}>0$ and $\displaystyle r\in[0,1)$.

A key role is played in our derivations by masterful results of Reznick [23],
who obtained bounds on the ratios of power sum symmetric functions in several
variables. Reznick’s results enable us to obtain, for series involving
$\displaystyle\Sigma$, bounds which are valid for all $\displaystyle d$.

In Section 2, we derive a complete asymptotic expansion for
$\displaystyle\Psi_{d}(\Sigma)$. We obtain in Proposition 2.1 a $\displaystyle
O(\cdot)$ bound for the remainder of the zonal polynomial series expansion of
$\displaystyle\Psi_{d}(\Sigma)$ upon truncation after any given number of
terms.

In Section 3, we obtain a complete asymptotic expansion for the matrix

$\displaystyle\Psi_{d}(\Sigma)\mathrm{Cov}(X)=\int_{\mathcal{S}^{d-1}}xx^{\prime}\exp(x^{\prime}\Sigma
x)\hskip 2.0pt\mathrm{d}x.$

We derive in Proposition 3.2 a $\displaystyle O(\cdot)$ bound for the
Frobenius norm of the remainder of that series expansion after truncation at
any number of terms. We obtain the asymptotic expansion for
$\displaystyle\mathrm{Cov}(X)$ by multiplying together the series for
$\displaystyle[\Psi_{d}(\Sigma)]^{-1}$ and
$\displaystyle\Psi_{d}(\Sigma)\mathrm{Cov}(X)$, and we provide in Theorem 3.3
the explicit expansion of $\displaystyle\mathrm{Cov}(X)$ up to terms of
arbitrary degree in $\displaystyle\Sigma$. We remark that our methods can also
be applied to obtain complete asymptotic expansions for higher moments of
$\displaystyle X$.

In Section 4, we describe the range of values of $\displaystyle d$ that have
appeared in the literature in a wide variety of applications of statistical
analyses of spherical data. We compute in Tables 1 and 2 numerous values of
the bounds for the remainder terms in the asymptotic expansions of
$\displaystyle\Psi_{d}(\Sigma)$ and $\displaystyle\nabla\Psi_{d}(\Sigma)$, the
matrix gradient of $\displaystyle\Psi_{d}(\Sigma)$, and we elucidate some
observed features of the tabulated values. Finally we provide all proofs in
the appendix Section A.

##  2 Asymptotic expansion of the normalizing constant

Denote by $\displaystyle\mathbb{R}_{sym}^{d\times d}$ the vector space of
$\displaystyle d\times d$ real symmetric matrices. As with any norm on a
vector space, the Frobenius norm, $\displaystyle\|\Sigma\|=[\hskip
1.0pt\mathrm{tr}\,(\Sigma^{2})]^{1/2}$, satisfies the subadditivity property,
$\displaystyle\|\Sigma_{1}+\Sigma_{2}\|\leq\|\Sigma_{1}\|+\|\Sigma_{2}\|$ for
all $\displaystyle\Sigma_{1},\Sigma_{2}\in\mathbb{R}_{sym}^{d\times d}$. The
Frobenius norm is also submultiplicative:
$\displaystyle\|\Sigma^{j}\|\leq\|\Sigma\|^{j}$ for all
$\displaystyle\Sigma\in\mathbb{R}_{sym}^{d\times d}$ and $\displaystyle
j=1,2,3,\ldots$ [17, p. 341].

Let $\displaystyle g_{1}$ be a real-valued function of $\displaystyle d$ and
$\displaystyle g_{2}$ be a nonnegative function of $\displaystyle d$. We use
the terminology $\displaystyle g_{1}(d)=O(g_{2}(d))$ as $\displaystyle
d\to\infty$ to signify that there exists a constant $\displaystyle\gamma>0$
such that $\displaystyle|g_{1}(d)|\leq\gamma g_{2}(d)$ for all sufficiently
large $\displaystyle d$.

Consider the situation in which $\displaystyle\|\Sigma\|\leq\gamma_{0}$ for
all $\displaystyle d$. By expanding the exponential function in (1.2) in a
Maclaurin expansion and integrating the first two terms of that series, we
find that $\displaystyle\Psi_{d}(\Sigma)\approx 1+d^{-1}\hskip
1.0pt\mathrm{tr}\,(\Sigma)+O(d^{-2})$. However, if $\displaystyle\|\Sigma\|$
is allowed to grow as $\displaystyle d$ increases then a deeper analysis is
required to derive the order of approximation of the remainder term.

Thus, our governing assumption throughout the article is that
$\displaystyle\|\Sigma\|\leq\gamma_{0}d^{r/2}$ for all $\displaystyle d$, with
constants $\displaystyle\gamma_{0}>0$ and $\displaystyle r\in[0,1)$.

For $\displaystyle a\in\mathbb{R}$, define the rising factorial
$\displaystyle(a)_{k}=a(a+1)\cdots(a+k-1)$, $\displaystyle k=0,1,2,\ldots$.
The zonal polynomial corresponding to the partition $\displaystyle(k)$ may be
defined as

$C_{(k)}(\Sigma)=\frac{k!}{(1/2)_{k}}\
\sum_{i_{1}+2i_{2}+3i_{3}+\cdots+ki_{k}=k}\ \prod_{j=1}^{k}\frac{[\hskip
1.0pt\mathrm{tr}\,(\Sigma^{j})]^{i_{j}}}{i_{j}!\,(2j)^{i_{j}}},$ (2.1)

where the multiple sum is over all nonnegative integers $\displaystyle
i_{1},\ldots,i_{k}$ such that $\displaystyle
i_{1}+2i_{2}+3i_{3}+\cdots+ki_{k}=k$ (James [18, p. 493, Eq. (123)]). In
particular,

$C_{(k)}(\Sigma)=\begin{cases}1,&k=0\\\ \hskip
1.0pt\mathrm{tr}\,(\Sigma),&k=1\\\ \frac{1}{3}\big{[}\big{(}\hskip
1.0pt\mathrm{tr}\,(\Sigma)\big{)}^{2}+2\hskip
1.0pt\mathrm{tr}\,(\Sigma^{2})\big{]},&k=2\end{cases}.$ (2.2)

It is well-known that $\displaystyle\Psi_{d}(\Sigma)$, the normalizing
constant in (1.2), can be expressed in terms of a confluent hypergeometric
function of matrix argument. We refer to Bingham [4] for the case
$\displaystyle d=3$ and to Muirhead [21, p. 288, Example 7.8] for the general
case. The confluent hypergeometric function of matrix argument can be
expressed as an infinite series of zonal polynomials, leading to the result,

$\Psi_{d}(\Sigma)={}_{1}F_{1}(1/2;d/2;\Sigma)=\sum_{k=0}^{\infty}\frac{(1/2)_{k}}{(d/2)_{k}}\frac{C_{(k)}(\Sigma)}{k!};$
(2.3)

see James [18, p. 477], Muirhead [21, Chapter 7]. The series (2.3) converges
absolutely for each $\displaystyle\Sigma$ and uniformly on compact regions
(Gross and Richards [15, Theorem 6.3]).

For $\displaystyle m\geq 1$ we will use as an asymptotic approximation to
$\displaystyle\Psi_{d}(\Sigma)$ the sum of the first $\displaystyle m$ terms
of the series (2.3),

$\Psi_{d}(\Sigma)\approx\sum_{k=0}^{m-1}\frac{(1/2)_{k}}{(d/2)_{k}}\frac{C_{(k)}(\Sigma)}{k!}.$
(2.4)

In studying the growth rate of the remainder series,

$R_{m}(\Sigma)=\sum_{k=m}^{\infty}\frac{(1/2)_{k}}{(d/2)_{k}}\frac{C_{(k)}(\Sigma)}{k!},$
(2.5)

we shall encounter the constants

$\gamma_{1}=\tfrac{1}{2}(1+\sqrt{3})\simeq 1.3660,\
\gamma_{2}=\gamma_{0}\gamma_{1},\ \hbox{and}\
\gamma_{3}=2^{3/2}e^{1/2}/\gamma_{1}\simeq 3.4138.$ (2.6)

Then we shall establish the following result.

######  Proposition 2.1.

Suppose that $\displaystyle\|\Sigma\|\leq\gamma_{0}d^{r/2}$, where
$\displaystyle\gamma_{0}>0$ and $\displaystyle 0\leq r<1$. Then, for all
$\displaystyle m\geq 1$ and all $\displaystyle
d\geq(2\gamma_{2}^{2})^{1/(1-r)}$,

$|R_{m}(\Sigma)|\leq\frac{\gamma_{3}\gamma_{2}^{m}}{[(m+1)!]^{1/2}}d^{-m(1-r)/2}.$
(2.7)

In particular, $\displaystyle R_{m}(\Sigma)=O\big{(}d^{-m(1-r)/2}\big{)}$ as
$\displaystyle d\to\infty$.

We note that a consequence of the inequality (2.7) is that $\displaystyle
R_{m}(\Sigma)$ converges to zero, i.e., the accuracy of the approximation
(2.4) increases, as $\displaystyle d$ increases (with $\displaystyle m$ held
fixed) or as $\displaystyle m$ increases (for fixed $\displaystyle d$). The
upper bound (2.7) also increases as $\displaystyle\gamma_{0}$ or
$\displaystyle r$ increases; this is consistent with the fact that the
condition $\displaystyle\|\Sigma\|\leq\gamma_{0}d^{r/2}$ entails an
increasingly larger matrix region as $\displaystyle\gamma_{0}$ or
$\displaystyle r$ increases.

As a consequence of Proposition 2.1, we obtain the following asymptotic
expansion for $\displaystyle[\Psi_{d}(\Sigma)]^{-1}$.

######  Corollary 2.2.

Suppose that $\displaystyle\|\Sigma\|\leq\gamma_{0}d^{r/2}$, where
$\displaystyle\gamma_{0}>0$ and $\displaystyle 0\leq r<1$. Then for all
$\displaystyle l\geq 2$ and all $\displaystyle d>(6\gamma_{2}^{2})^{1/(1-r)}$,

$[\Psi_{d}(\Sigma)]^{-1}=1-\sum_{j=1}^{l-1}\frac{(1/2)_{j}}{(d/2)_{j}}\frac{C_{(j)}(\Sigma)}{j!}+O(d^{-(1-r)}).$
(2.8)

##  3 Asymptotic expansion of the covariance matrix

It is well-known that $\displaystyle\mathbb{R}_{sym}^{d\times d}$, the vector
space of symmetric matrices, is an inner product space with the natural inner
product $\displaystyle\langle H_{1},H_{2}\rangle=\hskip
1.0pt\mathrm{tr}\,(H_{1}H_{2})$, $\displaystyle
H_{1},H_{2}\in\mathbb{R}_{sym}^{d\times d}$. This inner product is also
related to the Frobenius norm through the identity $\displaystyle\langle
H,H\rangle=\hskip 1.0pt\mathrm{tr}\,(H^{2})=\|H\|^{2}$, $\displaystyle
H\in\mathbb{R}_{sym}^{d\times d}$.

Let $\displaystyle C^{1}(\mathbb{R}_{sym}^{d\times d})$ be the set of all
functions $\displaystyle f:\mathbb{R}_{sym}^{d\times d}\to\mathbb{R}$ such
that all first-order partial derivatives, $\displaystyle\partial
f(\Sigma)/\partial\sigma_{i,j}$, are continuous, where
$\displaystyle\Sigma=(\sigma_{i,j})$. Let $\displaystyle\delta_{i,j}$ denote
Kronecker’s delta, i.e., $\displaystyle\delta_{i,j}=1$ or $\displaystyle 0$
according as $\displaystyle i=j$ or $\displaystyle i\neq j$, respectively;
then we define $\displaystyle\nabla$, the gradient operator, by

$\displaystyle\nabla
f(\Sigma)=\bigg{(}\tfrac{1}{2}(1+\delta_{i,j})\frac{\partial}{\partial\sigma_{i,j}}f(\Sigma)\bigg{)}.$

The operator $\displaystyle\nabla$ has appeared in various articles in
multivariate analysis, e.g., Hadjicosta and Richards [16], Richards [24],
Sebastiani [25], and it arises in the theory of Taylor expansions on
$\displaystyle\mathbb{R}_{sym}^{d\times d}$, as follows: For
$\displaystyle\Sigma,H\in\mathbb{R}_{sym}^{d\times d}$ and $\displaystyle f\in
C^{1}(\mathbb{R}_{sym}^{d\times d})$, there holds the first-order Taylor
expansion,

$f(\Sigma+H)=f(\Sigma)+\langle H,\nabla
f(\Sigma)\rangle+O(\|H\|^{2}),\qquad\|H\|\to 0.$ (3.1)

Turning to the covariance matrix of the Bingham-distributed random vector
$\displaystyle X$, it follows from (1.1) that $\displaystyle X$ has the same
distribution as $\displaystyle-X$; therefore, by symmetry,
$\displaystyle\mathbb{E}(X)=0$. Hence $\displaystyle\mathrm{Cov}(X)$, the
covariance matrix of $\displaystyle X$, equals
$\displaystyle\mathbb{E}(XX^{\prime})$. By (3.2) and a straightforward
justification of the interchange of derivatives and integral, we obtain

$\displaystyle\displaystyle\mathrm{Cov}(X)$
$\displaystyle\displaystyle=[\Psi_{d}(\Sigma)]^{-1}\int_{\mathcal{S}^{d-1}}xx^{\prime}\exp(x^{\prime}\Sigma
x)\hskip 2.0pt\mathrm{d}x$
$\displaystyle\displaystyle=[\Psi_{d}(\Sigma)]^{-1}\nabla\int_{\mathcal{S}^{d-1}}\exp\big{(}\hskip
1.0pt\mathrm{tr}\,(\Sigma xx^{\prime})\big{)}\hskip 2.0pt\mathrm{d}x$
$\displaystyle\displaystyle=[\Psi_{d}(\Sigma)]^{-1}\nabla\Psi_{d}(\Sigma).$

Our asymptotic approximations to the matrix
$\displaystyle\nabla\Psi_{d}(\Sigma)$ through term-by-term differentiation of
the zonal polynomial series (2.3) will be justified by means of bounds for
$\displaystyle\|\nabla C_{(k)}(\Sigma)\|$, the Frobenius norm of the matrix
$\displaystyle\nabla C_{(k)}(\Sigma)$. In order to calculate bounds for each
$\displaystyle\nabla C_{(k)}(\Sigma)$, we will need the following properties
of the gradient operator.

######  Lemma 3.1.

Let $\displaystyle H$ be a $\displaystyle d\times d$ symmetric matrix. Then

$\nabla\exp\big{(}\hskip 1.0pt\mathrm{tr}\,(\Sigma
H)\big{)}=H\exp\big{(}\hskip 1.0pt\mathrm{tr}\,(\Sigma H)\big{)}.$ (3.2)

Also, for $\displaystyle k=0,1,2,\ldots$,

$\nabla[\hskip 1.0pt\mathrm{tr}\,(\Sigma)]^{k}=k[\hskip
1.0pt\mathrm{tr}\,(\Sigma)]^{k-1}I_{d}$ (3.3)

and

$\nabla\hskip 1.0pt\mathrm{tr}\,(\Sigma^{k})=k\Sigma^{k-1}.$ (3.4)

In the sequel, we will obtain asymptotic approximations of
$\displaystyle\mathrm{Cov}(X)$ involving terms of arbitrary degree in
$\displaystyle\Sigma$. In deriving those results we note that, as a
consequence of (2.1), (3.3), and (3.4), the matrix $\displaystyle\nabla
C_{(k)}(\Sigma)$ is a homogeneous polynomial of degree $\displaystyle k-1$ for
$\displaystyle k\geq 1$.

For $\displaystyle m\geq 2$, we use the asymptotic approximation

$\nabla\Psi_{d}(\Sigma)\approx\sum_{k=1}^{m-1}\frac{(1/2)_{k}}{(d/2)_{k}}\frac{\nabla
C_{(k)}(\Sigma)}{k!}.$ (3.5)

Note that we necessarily impose the condition $\displaystyle m\geq 2$ since
$\displaystyle\nabla C_{(0)}(\Sigma)=0$, the zero matrix.

Following on the results in Section 2, we naturally define

$\displaystyle\nabla
R_{m}(\Sigma)=\sum_{k=m}^{\infty}\frac{(1/2)_{k}}{(d/2)_{k}}\frac{\nabla
C_{(k)}(\Sigma)}{k!}.$

To establish the accuracy of the matrix approximation (3.5), we shall derive
the convergence properties as $\displaystyle d\to\infty$ (and also as
$\displaystyle m\to\infty$) of $\displaystyle\nabla R_{m}(\Sigma)$. As a
consequence of those results, we deduce that $\displaystyle\nabla
R_{m}(\Sigma)$ is well-defined.

######  Proposition 3.2.

Suppose that $\displaystyle\|\Sigma\|\leq\gamma_{0}d^{r/2}$, where
$\displaystyle\gamma_{0}>0$ and $\displaystyle 0\leq r<1$. Then, for all
$\displaystyle m\geq 2$ and all $\displaystyle
d\geq(2\gamma_{2}^{2})^{1/(1-r)}$,

$\displaystyle\displaystyle\|\nabla
R_{m}(\Sigma)\|\leq\frac{(2e)^{1/2}\,\gamma_{2}^{m-1}}{[(m-1)!]^{1/2}}d^{-[1+(m-1)(1-r)]/2}.$
(3.6)

In particular, $\displaystyle\|\nabla
R_{m}(\Sigma)\|=O\big{(}d^{-[1+(m-1)(1-r)]/2}\big{)}$ as $\displaystyle
d\to\infty$.

As a consequence of (3.6), the accuracy of the approximation (3.5) increases
as $\displaystyle d\to\infty$ (with $\displaystyle m$ held fixed) or as
$\displaystyle m$ increases (with $\displaystyle d$ held fixed). Also, as with
(2.7), the bound in (3.6) increases as $\displaystyle\gamma_{0}$ or
$\displaystyle r$ increases.

We now obtain an asymptotic expansion of $\displaystyle\mathrm{Cov}(X)$ up to
terms of arbitrary degree in $\displaystyle\Sigma$.

######  Theorem 3.3.

Suppose that the random vector $\displaystyle X$ has a Bingham distribution
with the probability density function (1.1). Suppose also that
$\displaystyle\|\Sigma\|\leq\gamma_{0}d^{r/2}$, where
$\displaystyle\gamma_{0}>0$ and $\displaystyle 0\leq r<1$. Then for all
$\displaystyle l,m\geq 2$ and all $\displaystyle
d>(6\gamma_{2}^{2})^{1/(1-r)}$,

$\mathrm{Cov}(X)=\bigg{(}1-\sum_{j=1}^{l-1}\frac{(1/2)_{j}}{(d/2)_{j}}\frac{C_{(j)}(\Sigma)}{j!}\bigg{)}\bigg{(}\sum_{k=1}^{m-1}\frac{(1/2)_{k}}{(d/2)_{k}}\frac{\nabla
C_{(k)}(\Sigma)}{k!}\bigg{)}+O(d^{-\alpha}),$ (3.7)

where

$\displaystyle\alpha=\tfrac{1}{2}\min\\{3-2r,1+(m-1)(1-r)\\}=\begin{cases}\tfrac{1}{2}(2-r),&m=2\\\
\tfrac{1}{2}(3-2r),&m\geq 3\end{cases}.$

As we showed before, $\displaystyle C_{(j)}(\Sigma)$ and $\displaystyle\nabla
C_{(k)}(\Sigma)$ are homogeneous polynomials of degrees $\displaystyle j-1$
and $\displaystyle k-1$, respectively. Therefore the highest degree of
$\displaystyle\Sigma$ appearing in (3.7) equals $\displaystyle l-1+m-2=l+m-3$.

Consider the asymptotic expansion of $\displaystyle\mathrm{Cov}(X)$ arising
from (3.7) with $\displaystyle l=2$ and $\displaystyle m=3$. Applying (3.3),
the rules for derivatives of polynomials in $\displaystyle\Sigma$, to (2.2),
the formulas for the zonal polynomials $\displaystyle C_{(1)}(\Sigma)$ and
$\displaystyle C_{(2)}(\Sigma)$, we have

$\displaystyle\nabla C_{(1)}(\Sigma)=\nabla\hskip
1.0pt\mathrm{tr}\,\Sigma=I_{d}$

and

$\displaystyle\nabla C_{(2)}(\Sigma)=\frac{1}{3}\nabla[(\hskip
1.0pt\mathrm{tr}\,\Sigma)^{2}+2\hskip
1.0pt\mathrm{tr}\,(\Sigma^{2})]=\frac{2}{3}[(\hskip
1.0pt\mathrm{tr}\,\Sigma)I_{d}+2\Sigma].$

Substituting these expressions into (3.7) and collecting terms we obtain

$\displaystyle\displaystyle\mathrm{Cov}(X)$
$\displaystyle\displaystyle=\Big{(}1-\frac{1}{d}\hskip
1.0pt\mathrm{tr}\,(\Sigma)\Big{)}\Big{(}\frac{1}{d}I_{d}+\frac{1}{d(d+2)}[(\hskip
1.0pt\mathrm{tr}\,\Sigma)I_{d}+2\Sigma]\Big{)}+O(d^{-(3-2r)/2})$
$\displaystyle\displaystyle=\frac{1}{d}I_{d}-\frac{2}{(d+2)}(\hskip
1.0pt\mathrm{tr}\,\Sigma)I_{d}+\frac{2}{d(d+2)}\Sigma$
$\displaystyle\displaystyle\quad-\frac{1}{d^{2}(d+2)}(\hskip
1.0pt\mathrm{tr}\,\Sigma)^{2}I_{d}-\frac{2}{d^{2}(d+2)}(\hskip
1.0pt\mathrm{tr}\,\Sigma)\Sigma+O(d^{-(3-2r)/2}).$

In some articles it was assumed, without loss of generality, that
$\displaystyle\hskip 1.0pt\mathrm{tr}\,(\Sigma)=0$; cf.,Kume and Walker [19].
In that case, the above approximation reduces to

$\displaystyle\mathrm{Cov}(X)=\frac{1}{d}I_{d}+\frac{2}{d(d+2)}\Sigma+O(d^{-(3-2r)/2}).$

##  4 Implications for applications

The results derived in Sections 2 and 3 provide explicit bounds on the
remainder terms in asymptotic expansions, for large values of $\displaystyle
d$, of the normalizing constant and the covariance matrix of the Bingham
distribution on $\displaystyle\mathcal{S}^{d-1}$. These results are, to the
best of our knowledge, the first such derivations in the literature, so we now
comment on the implications of those bounds for real-world applications.

To assess the implications of our results for applications, we surveyed the
literature on the Bingham distributions in a search for the range of reported
values of $\displaystyle d$. Dryden [11] applied the high-dimensional Bingham
distribution with $\displaystyle d=62,501$ to model cortical surfaces using
magnetic resonance images of the human brain. A related study by Brignell, et
al. [6] applied principal components analysis to reduce the data to dimensions
up to $\displaystyle 20$.

Sra [26] reviewed applications of high-dimensional spherical data analyses in
machine learning. Those areas include text clustering, gene expression data
analysis, feature extraction, and wireless communications, with dimensions up
to $\displaystyle d=1,000$. Brombin, Pesarin, and Salmaso [7] discussed
examples in the study of shape modeling, and simulations for which the
dimension of the data is as large as $\displaystyle d=50$.

Dai, Dorman, Dutta, and Maitra [10] applied spherical data analysis to the
characterization of changes in cerebral blood-flow when the human brain is not
subject to any stimuli or tasks, to the assessment of variability in
handwritten numerical digits, and to the identification of genetic pathways
that underlie cancer. Those topics involve dimensions ranging from
$\displaystyle d=20$ to $\displaystyle d=100$.

To summarize, the articles that we reviewed described applications with
dimensions ranging from small, $\displaystyle d=20$, to very large,
$\displaystyle d=62,501$.

Table 1: Bounds for $\displaystyle|R_{m}(\Sigma)|$ and $\displaystyle\|\nabla
R_{m}(\Sigma)\|$ with $\displaystyle(\gamma_{0},r)=(1,0.5)$.

(a) Bounds for $\displaystyle|R_{m}(\Sigma)|$ $\displaystyle d$ | $\displaystyle m=3$ | $\displaystyle m=6$ | $\displaystyle m=10$
---|---|---|---
20 | 0.18782 | $\displaystyle 0.00349$ | $\displaystyle 0.00001$
25 | 0.15887 | $\displaystyle 0.00250$ | $\displaystyle 0.00000$
50 | 0.09447 | $\displaystyle 0.00088$ | $\displaystyle 0.00000$
75 | 0.06970 | $\displaystyle 0.00048$ | $\displaystyle 0.00000$
100 | 0.05617 | $\displaystyle 0.00031$ | $\displaystyle 0.00000$
250 | 0.02825 | $\displaystyle 0.00008$ | $\displaystyle 0.00000$
500 | 0.01680 | $\displaystyle 0.00003$ | $\displaystyle 0.00000$
750 | 0.01239 | $\displaystyle 0.00002$ | $\displaystyle 0.00000$
1000 | 0.00999 | $\displaystyle 0.00001$ | $\displaystyle 0.00000$
5000 | 0.00299 | $\displaystyle 0.00000$ | $\displaystyle 0.00000$
10000 | 0.00178 | $\displaystyle 0.00000$ | $\displaystyle 0.00000$
25000 | 0.00089 | $\displaystyle 0.00000$ | $\displaystyle 0.00000$
50000 | 0.00053 | $\displaystyle 0.00000$ | $\displaystyle 0.00000$
62501 | 0.00045 | $\displaystyle 0.00000$ | $\displaystyle 0.00000$

(b) Bounds for $\displaystyle\|\nabla R_{m}(\Sigma)\|$ $\displaystyle d$ | $\displaystyle m=3$ | $\displaystyle m=6$ | $\displaystyle m=10$
---|---|---|---
20 | 0.15383 | $\displaystyle 0.00535$ | $\displaystyle 0.00002$
25 | 0.12306 | $\displaystyle 0.00362$ | $\displaystyle 0.00001$
50 | 0.06153 | $\displaystyle 0.00108$ | $\displaystyle 0.00000$
75 | 0.04102 | $\displaystyle 0.00053$ | $\displaystyle 0.00000$
100 | 0.03077 | $\displaystyle 0.00032$ | $\displaystyle 0.00000$
250 | 0.01231 | $\displaystyle 0.00006$ | $\displaystyle 0.00000$
500 | 0.00615 | $\displaystyle 0.00002$ | $\displaystyle 0.00000$
750 | 0.00410 | $\displaystyle 0.00001$ | $\displaystyle 0.00000$
1000 | 0.00308 | $\displaystyle 0.00001$ | $\displaystyle 0.00000$
5000 | 0.00062 | $\displaystyle 0.00000$ | $\displaystyle 0.00000$
10000 | 0.00031 | $\displaystyle 0.00000$ | $\displaystyle 0.00000$
25000 | 0.00012 | $\displaystyle 0.00000$ | $\displaystyle 0.00000$
50000 | 0.00006 | $\displaystyle 0.00000$ | $\displaystyle 0.00000$
62501 | 0.00005 | $\displaystyle 0.00000$ | $\displaystyle 0.00000$

Consider the bound in (2.7) for $\displaystyle R_{m}(\Sigma)$, the remainder
term in the expansion of $\displaystyle\Psi_{d}(\Sigma)$. With
$\displaystyle(\gamma_{0},r)=(1,0.5)$ and $\displaystyle m=3,6,10$ we
computed, for $\displaystyle d>(2\gamma_{2}^{2})^{1/(1-r)}\simeq 14$, that
bound. The computed values, given in Table 1(a), illustrate that for moderate
values of $\displaystyle r$, the bound (2.7) decreases quickly as
$\displaystyle m$ increases. In particular, for $\displaystyle m\geq 6$, all
entries in Table 1(a) are less than $\displaystyle 10^{-2}$.

We also computed values of the bound in (3.6) for $\displaystyle\|\nabla
R_{m}(\Sigma)\|$. Several such values are presented in Table 1(b), and they
also illustrate that, for moderate values of $\displaystyle r$, the values of
$\displaystyle\|\nabla R_{m}(\Sigma)\|$ can decrease quickly as $\displaystyle
d$ increases.

Table 2: Bounds for $\displaystyle|R_{m}(\Sigma)|$ and $\displaystyle\|\nabla
R_{m}(\Sigma)\|$ with $\displaystyle(\gamma_{0},r)=(1,0.75)$.

(a) Bounds for $\displaystyle|R_{m}(\Sigma)|$ $\displaystyle d$ | $\displaystyle m=3$ | $\displaystyle m=6$ | $\displaystyle m=10$
---|---|---|---
200 | 0.24357 | $\displaystyle 0.00587$ | $\displaystyle 0.00002$
225 | 0.23304 | $\displaystyle 0.00538$ | $\displaystyle 0.00001$
250 | 0.22401 | $\displaystyle 0.00497$ | $\displaystyle 0.00001$
275 | 0.21615 | $\displaystyle 0.00463$ | $\displaystyle 0.00001$
500 | 0.17274 | $\displaystyle 0.00295$ | $\displaystyle 0.00001$
750 | 0.14837 | $\displaystyle 0.00218$ | $\displaystyle 0.00000$
1000 | 0.13320 | $\displaystyle 0.00176$ | $\displaystyle 0.00000$
2000 | 0.10271 | $\displaystyle 0.00104$ | $\displaystyle 0.00000$
2500 | 0.09447 | $\displaystyle 0.00088$ | $\displaystyle 0.00000$
5000 | 0.07284 | $\displaystyle 0.00053$ | $\displaystyle 0.00000$
10000 | 0.05617 | $\displaystyle 0.00031$ | $\displaystyle 0.00000$
25000 | 0.03984 | $\displaystyle 0.00016$ | $\displaystyle 0.00000$
50000 | 0.03072 | $\displaystyle 0.00009$ | $\displaystyle 0.00000$
62501 | 0.02825 | $\displaystyle 0.00008$ | $\displaystyle 0.00000$

(b) Bounds for $\displaystyle\|\nabla R_{m}(\Sigma)\|$ $\displaystyle d$ | $\displaystyle m=3$ | $\displaystyle m=6$ | $\displaystyle m=10$
---|---|---|---
200 | 0.05785 | $\displaystyle 0.00261$ | $\displaystyle 0.00001$
225 | 0.05296 | $\displaystyle 0.00229$ | $\displaystyle 0.00001$
250 | 0.04893 | $\displaystyle 0.00203$ | $\displaystyle 0.00001$
275 | 0.04556 | $\displaystyle 0.00182$ | $\displaystyle 0.00001$
500 | 0.02910 | $\displaystyle 0.00093$ | $\displaystyle 0.00000$
750 | 0.02147 | $\displaystyle 0.00059$ | $\displaystyle 0.00000$
1000 | 0.01730 | $\displaystyle 0.00043$ | $\displaystyle 0.00000$
2000 | 0.01029 | $\displaystyle 0.00020$ | $\displaystyle 0.00000$
2500 | 0.00870 | $\displaystyle 0.00015$ | $\displaystyle 0.00000$
5000 | 0.00517 | $\displaystyle 0.00007$ | $\displaystyle 0.00000$
10000 | 0.00308 | $\displaystyle 0.00003$ | $\displaystyle 0.00000$
25000 | 0.00155 | $\displaystyle 0.00001$ | $\displaystyle 0.00000$
50000 | 0.00092 | $\displaystyle 0.00001$ | $\displaystyle 0.00000$
62501 | 0.00078 | $\displaystyle 0.00000$ | $\displaystyle 0.00000$

As we noted before, the upper bounds in (2.7) and (3.6) are increasing
functions of $\displaystyle r$. Therefore the upper bounds derived from any
$\displaystyle r<0.5$ will be less than the corresponding values in Table 1.

As regards values of $\displaystyle r$ greater than $\displaystyle 0.5$ we
note that if $\displaystyle r$ is very large, e.g., $\displaystyle r>0.9$, the
restriction $\displaystyle d>(2\gamma_{2}^{2})^{1/(1-r)}$ can result in
extremely large values of $\displaystyle d$, even in the tens of millions,
which have never arisen in real-world problems. Thus, we calculated the upper
bounds in (2.7) and (3.6) for $\displaystyle(\gamma_{0},r)=(1,0.75)$,
$\displaystyle m=3,6,10$, and $\displaystyle
d>(2\gamma_{2}^{2})^{1/(1-r)}\simeq 194$. These computed values are provided
in Table 2, and they also illustrate the rapid decrease in the values of
$\displaystyle|R_{m}(\Sigma)|$ and $\displaystyle\|\nabla R_{m}(\Sigma)\|$ as
$\displaystyle d$ increases.

It is noticeable that, in Table 1(a) with $\displaystyle m=3$, each entry is
greater than the corresponding entry in Table 1(b), and the order is reversed
for $\displaystyle m=6$ and $\displaystyle m=10$. Throughout Table 2(a),
however, every entry is greater than or equal to the corresponding entry in
Table 2(b). This phenomenon may be explained by comparing the upper bounds on
the right-hand sides of (2.7) and (3.6); on doing so, we determine that a
necessary and sufficient condition for the bound in (2.7) to be less than or
equal to the bound in (3.6) is that

$d^{r}\geq\frac{m(m+1)}{4\gamma_{0}}.$ (4.1)

Consider Table 1, where we have $\displaystyle(\gamma_{0},r)=(1,0.5)$. Then
the condition (4.1) is satisfied for $\displaystyle m=3$ and all
$\displaystyle d\geq 9$; this explains why each entry in the first column of
Table 1(a) is greater than the corresponding entry in Table 1(b). For
$\displaystyle m=6$ the condition (4.1) is satisfied in Table 1 for
$\displaystyle d\geq 111$, and this is illustrated by comparing the entries of
the second column of that table. Similar remarks can be made about Table 2
and, more generally, the condition (4.1) implies that, for all sufficiently
large $\displaystyle d$, the bound provided by (3.6) will be less than the
corresponding bound from (2.7).

We conclude this section by suggesting a method for choosing $\displaystyle
m$. A statistician undertaking an application of the high-dimensional Bingham
distribution first should conduct exploratory work to determine plausible
values for $\displaystyle\gamma_{0}$ and $\displaystyle r$. Next, the
statistician should choose a tolerance or cutoff level, denoted by
$\displaystyle\epsilon$, for the larger of the upper bounds in Propositions
2.1 and 3.2; the value of $\displaystyle\epsilon$ should reflect the
importance of the application, so that $\displaystyle\epsilon$ presumably
would be chosen to be small in the case of crucial biomedical research. Once
$\displaystyle\gamma_{0}$, $\displaystyle r$, and $\displaystyle\epsilon$ have
been chosen then $\displaystyle m$ should be chosen so that the larger of
those two upper bounds is no greater than $\displaystyle\epsilon$, i.e.,

$\displaystyle\max\left\\{\frac{\gamma_{3}\gamma_{2}^{m}}{[(m+1)!]^{1/2}}d^{-m(1-r)/2},\frac{(2e)^{1/2}\,\gamma_{2}^{m-1}}{[(m-1)!]^{1/2}}d^{-[1+(m-1)(1-r)]/2}\right\\}\leq\epsilon.$

Acknowledgments. The authors are grateful to the referees and the editors for
helpful comments on the initial version of this article.

Declarations. No funds, grants, or other support were received for conducting
this research. The authors have no relevant financial or non-financial
interests to disclose. The authors have no conflicts of interest to declare.

## References

  * [1]
  * [2] Bagyan, A. (2015). Central Limit Theorems for Randomly Modulated Sequences of Random Vectors with Resampling and Applications to Statistics. Doctoral dissertation, Penn State University.
  * [3] Bhattacharya, A., and Bhattacharya, R. (2012). Nonparametric Inference on Manifolds, With Applications to Shape Spaces. Cambridge University Press, New York.
  * [4] Bingham, C. (1974). An antipodally symmetric distribution on the sphere. Annals of Statistics, 2, 1201–1225.
  * [5] Bingham, C., Chang, T., and Richards, D. (1992). Approximating the matrix Fisher and Bingham distributions: Applications to spherical regression and Procrustes analysis. Journal of Multivariate Analysis, 41, 314–337.
  * [6] Brignell, C. J., Dryden, I. L., Gattone, S. A., Park, B., Leask, S., Browne, W. J., and Flynn, S. (2010). Surface shape analysis with an application to brain surface asymmetry in schizophrenia. Biostatistics, 11, 609–630.
  * [7] Brombin, C., Pesarin, F., and Salmaso, L. (2011). Dealing with more variables than the sample size: An application to shape analysis. In: Nonparametric Statistics and Mixture Models: A Festschrift in Honor of Thomas P. Hettmansperger (eds. D. Hunter, D. Richards, and J. Rosenberger), pp. 28–44. World Scientific Press, Singapore.
  * [8] Chikuse, Y. (2012). Statistics on Special Manifolds. Lecture Notes in Statistics, Volume 174, Springer, New York.
  * [9] Comtet, L. (1974). Advanced Combinatorics: The Art of Finite and Infinite Expansions. D. Reidel, Dordrecht.
  * [10] Dai, F., and Dorman, K. S., Dutta, S., and Maitra, R. (2021). Exploratory factor analysis of data on a sphere. Preprint, arXiv:2111.04940.
  * [11] Dryden, I. L. (2003). Statistical analysis on high-dimensional spheres and shape spaces. Annals of Statistics, 33, 1643–1665.
  * [12] Dryden, I. L., and Mardia, K. V. (2016). Statistical Shape Analysis, with Applications in R. Wiley, Chichester.
  * [13] Dwyer, P. S., and Macphail, M. S. (1948). Symbolic matrix derivatives. Annals of Mathematical Statistics, 19, 517–534.
  * [14] Fisher, N. I., Lewis, T., and Embleton, B. J. J. (1993). Statistical Analysis of Spherical Data. Cambridge University Press, New York.
  * [15] Gross, K. I., and Richards, D. St. P. (1987). Special functions of matrix argument. I. Algebraic induction, zonal polynomials, and hypergeometric functions. Transactions of the American Mathematical Society, 301, 781–811.
  * [16] Hadjicosta, E., and Richards, D. (2020). Integral transform methods in goodness-of-fit testing, II: The Wishart distributions. Annals of the Institute of Statistical Mathematics, 72, 1317–1370.
  * [17] Horn, R. A., and Johnson, C. R. (2013). Matrix Analysis, second edition. Cambridge University Press, New York.
  * [18] James, A. T. (1964). Distributions of matrix variates and latent roots derived from normal samples. Annals of Mathematical Statistics, 35, 475–501.
  * [19] Kume, A., and Walker, S. G. (2014). On the Bingham distribution with large dimension. Journal of Multivariate Analysis, 124, 345–352.
  * [20] Mardia, K. V., and Jupp, P. E. (2000). Directional Statistics. Wiley, Chichester.
  * [21] Muirhead, R. J. (1982). Aspects of Multivariate Statistical Theory. Wiley, New York.
  * [22] Patrangenaru, V., and Ellingson, L. (2015). Nonparametric Statistics on Manifolds and Their Applications to Object Data Analysis. CRC Press, Boca Raton, FL.
  * [23] Reznick, B. (1983). Some inequalities for products of power sums. Pacific Journal of Mathematics, 104, 443–463.
  * [24] Richards, D. St. P. (1982). Differential operators associated with zonal polynomials. II. Annals of the Institute of Statistical Mathematics, 34, 119–121.
  * [25] Sebastiani, P. (1996). On the derivatives of matrix powers. SIAM Journal on Matrix Analysis and Applications, 17, 640–648.
  * [26] Sra, S. (2018). Directional statistics in machine learning: A brief review. In: Applied Directional Statistics: Modern Methods and Case Studies (eds. C. Ley and T. Verdebout), pp. 275–292. Chapman and Hall/CRC Press, New York.

##  Appendix A Appendix: Proofs

###  A.1 The proofs of Proposition 2.1 and Corollary 2.2

Proof of Proposition 2.1: By (2.5),

$|R_{m}(\Sigma)|\leq\sum_{k=m}^{\infty}\frac{(1/2)_{k}}{(d/2)_{k}}\frac{|C_{(k)}(\Sigma)|}{k!},$
(A.1)

and by (2.1),

$|C_{(k)}(\Sigma)|\leq\frac{k!}{(1/2)_{k}}\
\sum_{i_{1}+2i_{2}+3i_{3}+\cdots+ki_{k}=k}\ \prod_{j=1}^{k}\frac{|\hskip
1.0pt\mathrm{tr}\,(\Sigma^{j})|^{i_{j}}}{i_{j}!\,(2j)^{i_{j}}}.$ (A.2)

For each vector of indices $\displaystyle(i_{1},i_{2},\ldots,i_{k})$ such that
$\displaystyle i_{1}+2i_{2}+\cdots+ki_{k}=k$, we have

$\displaystyle\displaystyle\prod_{j=1}^{k}|\hskip
1.0pt\mathrm{tr}\,(\Sigma^{j})|^{i_{j}}$
$\displaystyle\displaystyle=\|\Sigma\|^{k}\prod_{j=1}^{k}\frac{|\hskip
1.0pt\mathrm{tr}\,(\Sigma^{j})|^{i_{j}}}{\|\Sigma\|^{ji_{j}}}$
$\displaystyle\displaystyle=\|\Sigma\|^{k}\prod_{j=1}^{k}\frac{|\hskip
1.0pt\mathrm{tr}\,(\Sigma^{j})|^{i_{j}}}{(\hskip
1.0pt\mathrm{tr}\,(\Sigma^{2}))^{ji_{j}/2}}.$ (A.3)

Denote by $\displaystyle\lambda_{1},\ldots,\lambda_{d}$ the eigenvalues of
$\displaystyle\Sigma$. Then for each $\displaystyle j=1,\ldots,k$,

$\displaystyle|\hskip
1.0pt\mathrm{tr}\,(\Sigma^{j})|^{i_{j}}=\bigg{|}\sum_{i=1}^{d}\lambda_{i}^{j}\bigg{|}^{i_{j}}\leq\bigg{(}\sum_{i=1}^{d}|\lambda_{i}|^{j}\bigg{)}^{i_{j}},$

and by substituting this bound into (A.1) we obtain

$\prod_{j=1}^{k}|\hskip
1.0pt\mathrm{tr}\,(\Sigma^{j})|^{i_{j}}\leq\|\Sigma\|^{k}\prod_{j=1}^{k}\frac{\left(\sum_{i=1}^{d}|\lambda_{i}|^{j}\right)^{i_{j}}}{\left(\sum_{i=1}^{d}|\lambda_{i}|^{2}\right)^{ji_{j}/2}}.$
(A.4)

By a remarkable result of Reznick [23, p. 447, eq. (2.12)] we obtain, for each
$\displaystyle j=1,\ldots,k$ and all $\displaystyle d$,

$\frac{\left(\sum_{i=1}^{d}|\lambda_{i}|^{j}\right)^{i_{j}}}{\left(\sum_{i=1}^{d}|\lambda_{i}|^{2}\right)^{ji_{j}/2}}\leq
d^{\max\\{0,(2-j)i_{j}/2\\}}=\begin{cases}d^{i_{1}/2},&j=1\\\ 1,&j\geq
2\end{cases}.$ (A.5)

Substituting the bound in (A.5) into (A.4), we obtain

$\prod_{j=1}^{k}|\hskip 1.0pt\mathrm{tr}\,(\Sigma^{j})|^{i_{j}}\leq
d^{i_{1}/2}\|\Sigma\|^{k}\equiv
d^{i_{1}/2}\|\Sigma\|^{i_{1}+2i_{2}+3i_{3}+\cdots+ki_{k}}.$ (A.6)

On applying (A.6) to (A.2) we obtain

$|C_{(k)}(\Sigma)|\leq\frac{k!}{(1/2)_{k}}a_{k}\|\Sigma\|^{k},$ (A.7)

where

$a_{k}=\sum_{i_{1}+2i_{2}+3i_{3}+\cdots+ki_{k}=k}d^{i_{1}/2}\prod_{j=1}^{k}\frac{1}{i_{j}!\,(2j)^{i_{j}}}.$
(A.8)

We will derive in Lemma A.1 a single-sum expression for the multiple sum
$\displaystyle a_{k}$, and then we will deduce from that single-sum the
inequality

$a_{k}\leq\frac{(d^{1/2}/2)_{k}}{k!}.$ (A.9)

Substituting (A.9) into (A.7) and recalling that
$\displaystyle\|\Sigma\|\leq\gamma_{0}d^{r/2}$ for all $\displaystyle d$, we
obtain

$|C_{(k)}(\Sigma)|\leq\frac{(d^{1/2}/2)_{k}}{(1/2)_{k}}\|\Sigma\|^{k}\leq\frac{(d^{1/2}/2)_{k}}{(1/2)_{k}}(\gamma_{0}d^{r/2})^{k}.$
(A.10)

Inserting (A.10) into (A.1) we obtain

$|R_{m}(\Sigma)|\leq\sum_{k=m}^{\infty}\frac{(d^{1/2}/2)_{k}}{(d/2)_{k}}\frac{\gamma_{0}^{k}d^{rk/2}}{k!}.$
(A.11)

As we shall show in Lemma A.2, with
$\displaystyle\gamma_{1}=\tfrac{1}{2}(\sqrt{3}+1)\simeq 1.366025$,

$\frac{(d^{1/2}/2)_{k}}{(d/2)_{k}}\leq\gamma_{1}^{k-1}((k-1)!)^{1/2}d^{-k/2}$
(A.12)

for all $\displaystyle d\geq 1$, $\displaystyle k\geq 1$. Letting
$\displaystyle\gamma_{2}=\gamma_{1}\gamma_{0}$ and substituting (A.12) into
(A.11), we find that for all $\displaystyle d$,

$|R_{m}(\Sigma)|\leq\gamma_{1}^{-1}\sum_{k=m}^{\infty}\frac{((k-1)!)^{1/2}}{k!}\,\gamma_{2}^{k}d^{-k(1-r)/2}.$
(A.13)

By applying the ratio test, we find that the series (A.13) converges
absolutely for all fixed $\displaystyle d$.

Note that although the $\displaystyle k$th term in (A.13) is $\displaystyle
O(d^{-k(1-r)/2})$, it is not evident that the sum of the resulting series is
$\displaystyle O(d^{-m(1-r)/2})$ as $\displaystyle d\to\infty$. To establish
that stated convergence rate we apply to (A.13) the Cauchy-Schwarz inequality,
obtaining

$\displaystyle\displaystyle|R_{m}(\Sigma)|$
$\displaystyle\displaystyle\leq\gamma_{1}^{-1}\Big{(}\sum_{k=m}^{\infty}\frac{(k-1)!}{(k!)^{2}}\Big{)}^{1/2}\Big{(}\sum_{k=m}^{\infty}\gamma_{2}^{2k}d^{-k(1-r)}\Big{)}^{1/2}$
$\displaystyle\displaystyle=\gamma_{1}^{-1}\Big{(}\sum_{k=m}^{\infty}\frac{1}{k\cdot
k!}\Big{)}^{1/2}\cdot\gamma_{2}^{m}d^{-m(1-r)/2}\big{(}1-\gamma_{2}^{2}d^{-(1-r)}\big{)}^{-1/2}.$
(A.14)

For $\displaystyle k\geq m$, it is straightforward that

$\displaystyle\frac{m}{k}\leq 1\leq\binom{k}{m}\equiv\frac{k!}{m!(k-m)!},$

equivalently,

$\displaystyle\frac{1}{k\cdot k!}\leq\frac{1}{m\cdot
m!}\cdot\frac{1}{(k-m)!},$

also that

$\displaystyle\frac{1}{m}\leq\frac{2}{m+1}.$

Therefore

$\displaystyle\frac{1}{k\cdot k!}\leq\frac{1}{m\cdot
m!}\cdot\frac{1}{(k-m)!}\leq\frac{2}{(m+1)!}\cdot\frac{1}{(k-m)!}.$

Summing this inequality over all $\displaystyle k\geq m$ we obtain

$\sum_{k=m}^{\infty}\frac{1}{k\cdot
k!}\leq\frac{2}{(m+1)!}\sum_{k=m}^{\infty}\frac{1}{(k-m)!}=\frac{2e}{(m+1)!}.$
(A.15)

Noting that the inequality $\displaystyle d\geq(2\gamma_{2}^{2})^{1/(1-r)}$ is
equivalent to

$\big{(}1-\gamma_{2}^{2}d^{-(1-r)}\big{)}^{-1}\leq 2,$ (A.16)

and applying (A.15) and (A.16) to (A.1) we obtain, for all $\displaystyle
d\geq(2\gamma_{2}^{2})^{1/(1-r)}$, the inequality

$|R_{m}(\Sigma)|\leq
2\gamma_{1}^{-1}\Big{(}\frac{2e}{(m+1)!}\Big{)}^{1/2}\gamma_{2}^{m}d^{-m(1-r)/2}.$
(A.17)

This establishes (2.7) and proves that $\displaystyle
R_{m}(\Sigma)=O(d^{-m(1-r)/2})$ as $\displaystyle d\to\infty$, and the proof
of the proposition now is complete. $\displaystyle\Box$

Proof of Corollary 2.2: By (A.17) with $\displaystyle m=1$ we have, for all
$\displaystyle d\geq(2\gamma_{2}^{2})^{1/(1-r)}$,

$|R_{1}(\Sigma)|\leq 2\gamma_{1}^{-1}e^{1/2}\gamma_{2}d^{-(1-r)/2}.$ (A.18)

For $\displaystyle d>(6\gamma_{2}^{2})^{1/(1-r)}$, we use the values of the
constants given in (2.6) to determine that the right-hand side of (A.18) is
strictly less than $\displaystyle 1$. Then we apply the geometric series to
obtain, for all $\displaystyle l\geq 2$,

$\displaystyle\displaystyle[\Psi_{d}(\Sigma)]^{-1}$
$\displaystyle\displaystyle=[1+R_{1}(\Sigma)]^{-1}$
$\displaystyle\displaystyle=1-R_{1}(\Sigma)+\sum_{j=2}^{\infty}(-1)^{j}\big{(}R_{1}(\Sigma)\big{)}^{j}$
$\displaystyle\displaystyle=1-\bigg{(}\sum_{j=1}^{l-1}\frac{(1/2)_{j}}{(d/2)_{j}}\frac{C_{(j)}(\Sigma)}{j!}+R_{l}(\Sigma)\bigg{)}+\sum_{j=2}^{\infty}(-1)^{j}\big{(}R_{1}(\Sigma)\big{)}^{j}.$
(A.19)

By Proposition 2.1, $\displaystyle R_{l}(\Sigma)=O(d^{-l(1-r)/2})$. On
applying (A.17) with $\displaystyle m=1$ we obtain, for all $\displaystyle
d>(6\gamma_{2}^{2})^{1/(1-r)}$,

$\displaystyle\displaystyle\left|\sum_{j=2}^{\infty}(-1)^{j}\big{(}R_{1}(\Sigma)\big{)}^{j}\right|$
$\displaystyle\displaystyle\leq\sum_{j=2}^{\infty}|R_{1}(\Sigma)|^{j}$
$\displaystyle\displaystyle\leq\sum_{j=2}^{\infty}\big{[}2\gamma_{1}^{-1}\gamma_{2}e^{1/2}d^{-(1-r)/2}\big{]}^{j}.$
(A.20)

By summing this geometric series, we deduce that (A.1) is $\displaystyle
O(d^{-(1-r)})$ as $\displaystyle d\to\infty$. Therefore by (A.1), for all
$\displaystyle l\geq 2$ and all $\displaystyle d>(6\gamma_{2}^{2})^{1/(1-r)}$,

$\displaystyle[\Psi_{d}(\Sigma)]^{-1}=1-\sum_{j=1}^{l-1}\frac{(1/2)_{j}}{(d/2)_{j}}\frac{C_{(j)}(\Sigma)}{j!}+O(d^{-l(1-r)/2})+O(d^{-(1-r)}).$

On applying the well-known property,

$O(d^{-a})+O(d^{-b})=O(d^{-\min(a,b)})$ (A.21)

for $\displaystyle a,b>0$, we obtain (2.8). $\displaystyle\Box$

###  A.2 The proof of Proposition 3.2 and Theorem 3.3

Proof of Lemma 3.1. The result (3.2) follows by straightforwardly applying
each element of the matrix $\displaystyle\nabla$ to the function
$\displaystyle\exp\big{(}\hskip 1.0pt\mathrm{tr}\,(\Sigma H)\big{)}$.

The formula (3.3) follows from the chain rule:

$\displaystyle\nabla[\hskip 1.0pt\mathrm{tr}\,(\Sigma)]^{k}=k[\hskip
1.0pt\mathrm{tr}\,(\Sigma)]^{k-1}\nabla\hskip
1.0pt\mathrm{tr}\,(\Sigma)=k[\hskip 1.0pt\mathrm{tr}\,(\Sigma)]^{k-1}I_{d}.$

As for (3.4), that result can be deduced by arguments similar to those given
by Sebastiani [25, Lemmas 3.1 or 5.1]. Also see Dwyer and Macphail [13, p.
528, Section 14] for the derivatives of the trace of powers of square, non-
symmetric matrices; calculations similar to theirs can also lead to (3.4).

A succinct derivation of (3.4) is obtained from the Taylor expansion (3.1), as
follows. By using the commutativity property of the trace, i.e.,
$\displaystyle\hskip 1.0pt\mathrm{tr}\,(\Sigma H)=\hskip
1.0pt\mathrm{tr}\,(H\Sigma)$, we find that

$\displaystyle\displaystyle\hskip 1.0pt\mathrm{tr}\,[(\Sigma+H)^{k}]$
$\displaystyle\displaystyle=\hskip
1.0pt\mathrm{tr}\,[(\Sigma+H)(\Sigma+H)\cdots(\Sigma+H)]$
$\displaystyle\displaystyle=\hskip 1.0pt\mathrm{tr}\,(\Sigma^{k})+k\hskip
1.0pt\mathrm{tr}\,(H\Sigma^{k-1})+O(\|H\|^{2})$
$\displaystyle\displaystyle\equiv\hskip 1.0pt\mathrm{tr}\,(\Sigma^{k})+\langle
H,k\Sigma^{k-1}\rangle+O(\|H\|^{2})$ (A.22)

as $\displaystyle H\to 0$. Setting $\displaystyle f(\Sigma)=\hskip
1.0pt\mathrm{tr}\,(\Sigma^{k})$ in (3.1) and comparing the result with (A.2),
we obtain (3.4). $\displaystyle\Box$

Proof of Proposition 3.2: By the subadditivity property of the Frobenius norm,
we have

$\displaystyle\displaystyle\|\nabla R_{m}(\Sigma)\|$
$\displaystyle\displaystyle=\bigg{\|}\sum_{k=m}^{\infty}\frac{(1/2)_{k}}{(d/2)_{k}}\frac{\nabla
C_{(k)}(\Sigma)}{k!}\bigg{\|}$
$\displaystyle\displaystyle\leq\sum_{k=m}^{\infty}\bigg{\|}\frac{(1/2)_{k}}{(d/2)_{k}}\frac{\nabla
C_{(k)}(\Sigma)}{k!}\bigg{\|}.$

For $\displaystyle k\geq 2$, it follows from (2.1) that

$\displaystyle\displaystyle C_{(k)}(\Sigma)$
$\displaystyle\displaystyle=\frac{k!}{(1/2)_{k}}\Bigg{[}\frac{[\hskip
1.0pt\mathrm{tr}\,(\Sigma)]^{k}}{k!\,2^{k}}+\sum_{\begin{subarray}{c}i_{1}+2i_{2}+3i_{3}+\cdots+ki_{k}=k\\\
i_{1}\leq k-2\end{subarray}}\ \prod_{j=1}^{k}\frac{[\hskip
1.0pt\mathrm{tr}\,(\Sigma^{j})]^{i_{j}}}{i_{j}!\,(2j)^{i_{j}}}\Bigg{]}$
$\displaystyle\displaystyle\equiv\frac{[\hskip
1.0pt\mathrm{tr}\,(\Sigma)]^{k}}{(1/2)_{k}\,2^{k}}+\frac{k!}{(1/2)_{k}}p_{k}(\Sigma),$
(A.23)

where

$p_{k}(\Sigma)=\sum_{\begin{subarray}{c}i_{1}+2i_{2}+3i_{3}+\cdots+ki_{k}=k\\\
i_{1}\leq k-2\end{subarray}}\ \prod_{j=1}^{k}\frac{[\hskip
1.0pt\mathrm{tr}\,(\Sigma^{j})]^{i_{j}}}{i_{j}!\,(2j)^{i_{j}}}$ (A.24)

is a polynomial in $\displaystyle\\{\hskip 1.0pt\mathrm{tr}\,(\Sigma),\hskip
1.0pt\mathrm{tr}\,(\Sigma^{2}),\ldots,\hskip
1.0pt\mathrm{tr}\,(\Sigma^{k})\\}$. It follows from (A.24) that $\displaystyle
p_{k}(\Sigma)$ is homogeneous of degree $\displaystyle k$ in
$\displaystyle\Sigma$ and its coefficients do not depend on $\displaystyle d$;
moreover, because of the restriction $\displaystyle i_{1}\leq k-2$,d the
highest power of $\displaystyle\hskip 1.0pt\mathrm{tr}\,(\Sigma)$ which can
appear in $\displaystyle p_{k}(\Sigma)$ is $\displaystyle k-2$.

By (3.3) and (A.2),

$\displaystyle\displaystyle\frac{(1/2)_{k}}{(d/2)_{k}}\frac{\nabla
C_{(k)}(\Sigma)}{k!}$
$\displaystyle\displaystyle=\frac{(1/2)_{k}}{(d/2)_{k}\,k!}\nabla\bigg{[}\frac{[\hskip
1.0pt\mathrm{tr}\,(\Sigma)]^{k}}{(1/2)_{k}\,2^{k}}+\frac{k!}{(1/2)_{k}}p_{k}(\Sigma)\bigg{]}$
$\displaystyle\displaystyle=\frac{1}{(d/2)_{k}}\bigg{[}\frac{1}{(k-1)!\,2^{k}}[\hskip
1.0pt\mathrm{tr}\,(\Sigma)]^{k-1}I_{d}+\nabla p_{k}(\Sigma)\bigg{]}.$ (A.25)

On applying to (A.2) the subadditivity property of the Frobenius norm we
obtain

$\displaystyle\displaystyle\bigg{\|}\frac{(1/2)_{k}}{(d/2)_{k}}$
$\displaystyle\displaystyle\frac{\nabla C_{(k)}(\Sigma)}{k!}\bigg{\|}$
$\displaystyle\displaystyle\leq\frac{1}{(d/2)_{k}}\bigg{[}\frac{1}{(k-1)!\,2^{k}}|\hskip
1.0pt\mathrm{tr}\,(\Sigma)|^{k-1}\|I_{d}\|+\|\nabla p_{k}(\Sigma)\|\bigg{]}$
$\displaystyle\displaystyle=\frac{1}{(d/2)_{k}}\bigg{[}\frac{d^{1/2}}{(k-1)!\,2^{k}}|\hskip
1.0pt\mathrm{tr}\,(\Sigma)|^{k-1}+\|\nabla p_{k}(\Sigma)\|\bigg{]}.$ (A.26)

By (A.24) and the product rule for derivatives,

$\displaystyle\displaystyle\nabla p_{k}(\Sigma)$
$\displaystyle\displaystyle=\sum_{\begin{subarray}{c}i_{1}+2i_{2}+3i_{3}+\cdots+ki_{k}=k\\\
i_{1}\leq k-2\end{subarray}}\,\nabla\prod_{j=1}^{k}\frac{[\hskip
1.0pt\mathrm{tr}\,(\Sigma^{j})]^{i_{j}}}{i_{j}!\,(2j)^{i_{j}}}$
$\displaystyle\displaystyle=\sum_{\begin{subarray}{c}i_{1}+2i_{2}+3i_{3}+\cdots+ki_{k}=k\\\
i_{1}\leq
k-2\end{subarray}}\,\sum_{l=1}^{k}\Big{(}\prod_{\begin{subarray}{c}j=1\\\
j\neq l\end{subarray}}^{k}\frac{[\hskip
1.0pt\mathrm{tr}\,(\Sigma^{j})]^{i_{j}}}{i_{j}!\,(2j)^{i_{j}}}\Big{)}\cdot\frac{\nabla[\hskip
1.0pt\mathrm{tr}\,(\Sigma^{l})]^{i_{l}}}{i_{l}!\,(2l)^{i_{l}}}.$

By (3.4) and the chain rule,

$\displaystyle\nabla[\hskip
1.0pt\mathrm{tr}\,(\Sigma^{l})]^{i_{l}}=i_{l}[\hskip
1.0pt\mathrm{tr}\,(\Sigma^{l})]^{i_{l}-1}\nabla[\hskip
1.0pt\mathrm{tr}\,(\Sigma^{l})]=l\,i_{l}[\hskip
1.0pt\mathrm{tr}\,(\Sigma^{l})]^{i_{l}-1}\Sigma^{l-1};$

hence,

$\displaystyle\nabla
p_{k}(\Sigma)=\sum_{\begin{subarray}{c}i_{1}+2i_{2}+\cdots+ki_{k}=k\\\
i_{1}\leq
k-2\end{subarray}}\sum_{l=1}^{k}\frac{l\,i_{l}}{i_{l}!(2l)^{i_{l}}}\Big{(}\prod_{\begin{subarray}{c}j=1\\\
j\neq l\end{subarray}}^{k}\frac{[\hskip
1.0pt\mathrm{tr}\,(\Sigma^{j})]^{i_{j}}}{i_{j}!(2j)^{i_{j}}}\Big{)}[\hskip
1.0pt\mathrm{tr}\,(\Sigma^{l})]^{i_{l}-1}\Sigma^{l-1},$

which also reveals that $\displaystyle\nabla p_{k}(\Sigma)$ is homogeneous of
degree $\displaystyle k-1$. On applying the subadditivity property of the
Frobenius norm, we obtain

$\displaystyle\displaystyle\|\nabla p_{k}(\Sigma)\|$
$\displaystyle\displaystyle\
\leq\sum_{\begin{subarray}{c}i_{1}+2i_{2}+\cdots+ki_{k}=k\\\ i_{1}\leq
k-2\end{subarray}}\bigg{(}\prod_{j=1}^{k}\frac{1}{i_{j}!(2j)^{i_{j}}}\bigg{)}\sum_{l=1}^{k}l\,i_{l}\Big{(}\prod_{\begin{subarray}{c}j=1\\\
j\neq l\end{subarray}}^{k}|\hskip
1.0pt\mathrm{tr}\,(\Sigma^{j})|^{i_{j}}\Big{)}|\hskip
1.0pt\mathrm{tr}\,(\Sigma^{l})|^{i_{l}-1}\|\Sigma^{l-1}\|.$

On applying (A.6) we have, for all $\displaystyle d$,

$\displaystyle\Big{(}\prod_{\begin{subarray}{c}j=1\\\ j\neq
l\end{subarray}}^{k}|\hskip
1.0pt\mathrm{tr}\,(\Sigma^{j})|^{i_{j}}\Big{)}|\hskip
1.0pt\mathrm{tr}\,(\Sigma^{l})|^{i_{l}-1}\leq
d^{i_{1}/2}\|\Sigma\|^{i_{1}+2i_{2}+3i_{3}+\cdots+ki_{k}-l}=d^{i_{1}/2}\|\Sigma\|^{k-l}.$

By also applying the submultiplicative inequality,
$\displaystyle\|\Sigma^{l-1}\|\leq\|\Sigma\|^{l-1}$, we obtain

$\displaystyle\displaystyle\|\nabla p_{k}(\Sigma)\|$
$\displaystyle\displaystyle\leq\sum_{\begin{subarray}{c}i_{1}+2i_{2}+\cdots+ki_{k}=k\\\
i_{1}\leq
k-2\end{subarray}}\bigg{(}\prod_{j=1}^{k}\frac{1}{i_{j}!\,(2j)^{i_{j}}}\bigg{)}\sum_{l=1}^{k}li_{l}\,d^{i_{1}/2}\|\Sigma\|^{k-l}\|\Sigma\|^{l-1}$
$\displaystyle\displaystyle=\|\Sigma\|^{k-1}\sum_{\begin{subarray}{c}i_{1}+2i_{2}+\cdots+ki_{k}=k\\\
i_{1}\leq
k-2\end{subarray}}d^{i_{1}/2}\bigg{(}\prod_{j=1}^{k}\frac{1}{i_{j}!\,(2j)^{i_{j}}}\bigg{)}\sum_{l=1}^{k}l\,i_{l}$
$\displaystyle\displaystyle=k\|\Sigma\|^{k-1}\sum_{\begin{subarray}{c}i_{1}+2i_{2}+\cdots+ki_{k}=k\\\
i_{1}\leq
k-2\end{subarray}}d^{i_{1}/2}\prod_{j=1}^{k}\frac{1}{i_{j}!\,(2j)^{i_{j}}}.$
(A.27)

By (A.8) and (A.9),

$\displaystyle\displaystyle\sum_{\begin{subarray}{c}i_{1}+2i_{2}+\cdots+ki_{k}=k\\\
i_{1}\leq
k-2\end{subarray}}d^{i_{1}/2}\prod_{j=1}^{k}\frac{1}{i_{j}!\,(2j)^{i_{j}}}$
$\displaystyle\displaystyle\equiv a_{k}-\frac{d^{k/2}}{k!\,2^{k}}$
$\displaystyle\displaystyle\leq\frac{1}{k!}\big{[}(d^{1/2}/2)_{k}-(d^{1/2}/2)^{k}\big{]},$

and by substituting this bound into (A.2), we obtain

$\displaystyle\|\nabla
p_{k}(\Sigma)\|\leq\frac{1}{(k-1)!}\big{[}(d^{1/2}/2)_{k}-(d^{1/2}/2)^{k}\big{]}\|\Sigma\|^{k-1}.$

On applying the latter inequality at (A.2), we find that

$\bigg{\|}\frac{(1/2)_{k}}{(d/2)_{k}}\frac{\nabla
C_{(k)}(\Sigma)}{k!}\bigg{\|}\\\
\leq\frac{1}{(d/2)_{k}\,(k-1)!}\bigg{[}\frac{d^{1/2}}{2^{k}}|\hskip
1.0pt\mathrm{tr}\,(\Sigma)|^{k-1}+\big{[}(d^{1/2}/2)_{k}-(d^{1/2}/2)^{k}\big{]}\|\Sigma\|^{k-1}\bigg{]}.$

Setting $\displaystyle(i_{1},i_{2},\ldots,i_{k})=(1,0,\ldots,0)$ in (A.6), we
obtain

$\displaystyle|\hskip 1.0pt\mathrm{tr}\,(\Sigma)|\leq d^{1/2}[\hskip
1.0pt\mathrm{tr}\,(\Sigma^{2})]^{1/2}=d^{1/2}\|\Sigma\|;$

hence

$\displaystyle\displaystyle\bigg{\|}$
$\displaystyle\displaystyle\frac{(1/2)_{k}}{(d/2)_{k}}\frac{\nabla
C_{(k)}(\Sigma)}{k!}\bigg{\|}$
$\displaystyle\displaystyle\leq\frac{1}{(d/2)_{k}\,(k-1)!}\bigg{[}\frac{d^{1/2}}{2^{k}}(d^{1/2}\|\Sigma\|)^{k-1}+\big{[}(d^{1/2}/2)_{k}-(d^{1/2}/2)^{k}\big{]}\|\Sigma\|^{k-1}\bigg{]}$
$\displaystyle\displaystyle=\frac{(d^{1/2}/2)_{k}}{(d/2)_{k}}\frac{\|\Sigma\|^{k-1}}{(k-1)!}.$

On applying the inequalities $\displaystyle\|\Sigma\|\leq\gamma_{0}d^{r/2}$
and (A.32) we obtain

$\displaystyle\displaystyle\bigg{\|}\frac{(1/2)_{k}}{(d/2)_{k}}\frac{\nabla
C_{(k)}(\Sigma)}{k!}\bigg{\|}$
$\displaystyle\displaystyle\leq\frac{(d^{1/2}/2)_{k}}{(d/2)_{k}}\frac{(\gamma_{0}d^{r/2})^{k-1}}{(k-1)!}$
$\displaystyle\displaystyle\leq\gamma_{1}^{k-1}((k-1)!)^{1/2}d^{-k/2}\frac{(\gamma_{0}d^{r/2})^{k-1}}{(k-1)!}$
$\displaystyle\displaystyle=d^{-1/2}\frac{(\gamma_{2}d^{-(1-r)/2})^{k-1}}{((k-1)!)^{1/2}}.$

Therefore

$\displaystyle\displaystyle\|\nabla R_{m}(\Sigma)\|$
$\displaystyle\displaystyle\leq\sum_{k=m}^{\infty}\bigg{\|}\frac{(1/2)_{k}}{(d/2)_{k}}\frac{\nabla
C_{(k)}(\Sigma)}{k!}\bigg{\|}$ $\displaystyle\displaystyle\leq
d^{-1/2}\sum_{k=m}^{\infty}\frac{(\gamma_{2}d^{-(1-r)/2})^{k-1}}{((k-1)!)^{1/2}},$
(A.28)

and, by applying the ratio test, we find that the latter series converges for
all fixed $\displaystyle d$.

On applying the Cauchy-Schwarz inequality to (A.2), we obtain

$\displaystyle\|\nabla R_{m}(\Sigma)\|\leq
d^{-1/2}\Big{(}\sum_{k=m}^{\infty}\frac{1}{(k-1)!}\Big{)}^{1/2}\Big{(}\sum_{k=m}^{\infty}(\gamma_{2}^{2}d^{-(1-r)})^{k-1}\Big{)}^{1/2}.$

For any nonnegative integers $\displaystyle a$ and $\displaystyle b$, it is
elementary that

$\displaystyle\frac{1}{(a+b)!}\leq\frac{1}{a!b!};$

therefore

$\displaystyle\sum_{k=m}^{\infty}\frac{1}{(k-1)!}=\sum_{k=0}^{\infty}\frac{1}{(k+m-1)!}\leq\sum_{k=0}^{\infty}\frac{1}{k!\,(m-1)!}=\frac{e}{(m-1)!}.$

Also, for all $\displaystyle d$ such that
$\displaystyle\gamma_{2}^{2}d^{-(1-r)}\leq 1/2$, equivalently, $\displaystyle
d\geq(2\gamma_{2}^{2})^{1/(1-r)}$, we have

$\displaystyle\displaystyle\sum_{k=m}^{\infty}(\gamma_{2}^{2}d^{-(1-r)})^{k-1}$
$\displaystyle\displaystyle=(\gamma_{2}^{2}d^{-(1-r)})^{m-1}(1-\gamma_{2}^{2}d^{-(1-r)})^{-1}$
$\displaystyle\displaystyle\leq 2(\gamma_{2}^{2}d^{-(1-r)})^{m-1},$

and then we obtain

$\displaystyle\displaystyle\|\nabla R_{m}(\Sigma)\|$
$\displaystyle\displaystyle\leq
d^{-1/2}\Big{(}\frac{e}{(m-1)!}\Big{)}^{1/2}\big{(}2(\gamma_{2}^{2}d^{-(1-r)})^{m-1}\big{)}^{1/2}$
$\displaystyle\displaystyle=(2e)^{1/2}[(m-1)!]^{-1/2}\gamma_{2}^{m-1}d^{-[1+(m-1)(1-r)]/2}.$

Therefore $\displaystyle\|\nabla R_{m}(\Sigma)\|=O(d^{-[1+(m-1)(1-r)]/2})$ as
$\displaystyle d\to\infty$. $\displaystyle\Box$

Proof of Theorem 3.3: Since $\displaystyle\nabla C_{(0)}(\Sigma)=0$ then by
the zonal polynomial expansion (2.3) of $\displaystyle\Psi_{d}(\Sigma)$,

$\nabla\Psi_{d}(\Sigma)=\sum_{k=1}^{m-1}\frac{(1/2)_{k}}{(d/2)_{k}}\frac{\nabla
C_{(k)}(\Sigma)}{k!}+\nabla R_{m}(\Sigma)$ (A.29)

and, by Proposition 3.2, $\displaystyle\|\nabla
R_{m}(\Sigma)\|=O(d^{-[1+(m-1)(1-r)]/2})$ as $\displaystyle d\to\infty$.

On applying the asymptotic expansion of $\displaystyle[\Psi_{d}(\Sigma)]^{-1}$
given in (2.8), and using (A.29), we obtain

$\displaystyle\displaystyle\mathrm{Cov}(X)$
$\displaystyle\displaystyle=[\Psi_{d}(\Sigma)]^{-1}\nabla\Psi_{d}(\Sigma)$
$\displaystyle\displaystyle=\bigg{(}1-\sum_{j=1}^{l-1}\frac{(1/2)_{j}}{(d/2)_{j}}\frac{C_{(j)}(\Sigma)}{j!}+O(d^{-(1-r)})\bigg{)}$
$\displaystyle\displaystyle\qquad\cdot\bigg{(}\sum_{k=1}^{m-1}\frac{(1/2)_{k}}{(d/2)_{k}}\frac{\nabla
C_{(k)}(\Sigma)}{k!}+O(d^{-[1+(m-1)(1-r)]/2})\bigg{)}.$

Expanding this product, we obtain

$\displaystyle\displaystyle\mathrm{Cov}(X)$
$\displaystyle\displaystyle=\bigg{(}1-\sum_{j=1}^{l-1}\frac{(1/2)_{j}}{(d/2)_{j}}\frac{C_{(j)}(\Sigma)}{j!}\bigg{)}\bigg{(}\sum_{k=1}^{m-1}\frac{(1/2)_{k}}{(d/2)_{k}}\frac{\nabla
C_{(k)}(\Sigma)}{k!}\bigg{)}$
$\displaystyle\displaystyle\qquad+O(d^{-(1-r)})\cdot\bigg{(}\sum_{k=1}^{m-1}\frac{(1/2)_{k}}{(d/2)_{k}}\frac{\nabla
C_{(k)}(\Sigma)}{k!}\bigg{)}$ $\displaystyle\displaystyle\qquad+1\cdot
O(d^{-[1+(m-1)(1-r)]/2})$
$\displaystyle\displaystyle\qquad+\bigg{(}\sum_{j=1}^{l-1}\frac{(1/2)_{j}}{(d/2)_{j}}\frac{C_{(j)}(\Sigma)}{j!}\bigg{)}\cdot
O(d^{-[1+(m-1)(1-r)]/2})$ $\displaystyle\displaystyle\qquad+O(d^{-(1-r)})\cdot
O(d^{-[1+(m-1)(1-r)]/2}).$ (A.30)

On applying Proposition 3.2, we obtain

$\displaystyle\displaystyle\bigg{\|}\sum_{k=1}^{m-1}\frac{(1/2)_{k}}{(d/2)_{k}}\frac{\nabla
C_{(k)}(\Sigma)}{k!}\bigg{\|}$ $\displaystyle\displaystyle=\|\nabla
R_{1}(\Sigma)-\nabla R_{m}(\Sigma)\|$ $\displaystyle\displaystyle\leq\|\nabla
R_{1}(\Sigma)\|+\|\nabla R_{m}(\Sigma)\|$
$\displaystyle\displaystyle=O(d^{-1/2})+O(d^{-[1+(m-1)(1-r)]/2})$
$\displaystyle\displaystyle=O(d^{-1/2});$

therefore

$\displaystyle\displaystyle
O(d^{-(1-r)})\cdot\bigg{\|}\sum_{k=1}^{m-1}\frac{(1/2)_{k}}{(d/2)_{k}}\frac{\nabla
C_{(k)}(\Sigma)}{k!}\bigg{\|}$ $\displaystyle\displaystyle=O(d^{-(1-r)})\cdot
O(d^{-1/2})$ $\displaystyle\displaystyle=O(d^{-(3-2r)/2}).$

Similarly, by Proposition 2.1,

$\displaystyle\displaystyle\sum_{j=1}^{l-1}\frac{(1/2)_{j}}{(d/2)_{j}}\frac{C_{(j)}(\Sigma)}{j!}$
$\displaystyle\displaystyle=R_{1}(\Sigma)-R_{l}(\Sigma)$
$\displaystyle\displaystyle=O(d^{-(1-r)/2})+O(d^{-l(1-r)/2})$
$\displaystyle\displaystyle=O(d^{-(1-r)/2}),$

and therefore

$\displaystyle\displaystyle\bigg{(}\sum_{j=1}^{l-1}\frac{(1/2)_{j}}{(d/2)_{j}}\frac{C_{(j)}(\Sigma)}{j!}$
$\displaystyle\displaystyle\bigg{)}\cdot O(d^{-[1+(m-1)(1-r)]/2})$
$\displaystyle\displaystyle=O(d^{-(1-r)/2})\cdot O(d^{-[1+(m-1)(1-r)]/2})$
$\displaystyle\displaystyle=O(d^{-[1+m(1-r)]/2}).$

The last $\displaystyle O(\cdot)$ term in (A.2) is

$\displaystyle O(d^{-(1-r)})\cdot
O(d^{-[1+(m-1)(1-r)]/2})=O(d^{-[1+(m+1)(1-r)]/2}).$

Collecting all $\displaystyle O(\cdot)$ terms in (A.2), we obtain

$\displaystyle\displaystyle\mathrm{Cov}(X)$
$\displaystyle\displaystyle=\bigg{(}1-\sum_{j=1}^{m-1}\frac{(1/2)_{j}}{(d/2)_{j}}\frac{C_{(j)}(\Sigma)}{j!}\bigg{)}\bigg{(}\sum_{k=1}^{m-1}\frac{(1/2)_{k}}{(d/2)_{k}}\frac{\nabla
C_{(k)}(\Sigma)}{k!}\bigg{)}$
$\displaystyle\displaystyle\qquad+O(d^{-(3-2r)/2})+O(d^{-[1+(m-1)(1-r)]/2})$
$\displaystyle\displaystyle\qquad+O(d^{-[1+m(1-r)]/2})+O(d^{-[1+(m+1)(1-r)]/2}).$

On applying (A.21), we obtain (3.7). $\displaystyle\Box$

###  A.3 The proofs of (A.9) and (A.12)

First, we establish the inequality (A.9) for the coefficients $\displaystyle
a_{k}$ defined in (A.8).

######  Lemma A.1.

For $\displaystyle k=0,1,2,\ldots$,

$a_{k}=\sum_{l=0}^{k}\frac{\big{(}(d^{1/2}-1)/2\big{)}^{l}}{l!}\frac{(1/2)_{k-l}}{(k-l)!}\leq\frac{(d^{1/2}/2)_{k}}{k!}.$
(A.31)

Proof. First, we follow the approach of Comtet [9, p. 97, Eq. 2d] to derive an
explicit formula for $\displaystyle a_{k}$.

Let $\displaystyle t=(t_{1},t_{2},t_{3},\ldots)$ be a vector of
indeterminates, and define

$\displaystyle
a_{k}(t)=\sum_{i_{1}+2i_{2}+3i_{3}+\cdots+ki_{k}=k}\prod_{j=1}^{k}\frac{t_{j}^{i_{j}}}{i_{j}!\,(2j)^{i_{j}}},$

$\displaystyle k\geq 0$. For an indeterminate $\displaystyle u$, a formal
generating-function for the sequence
$\displaystyle\\{a_{k}(t),k=0,1,2,\ldots\\}$ is

$\displaystyle\displaystyle G_{t}(u)=\sum_{k=0}^{\infty}a_{k}(t)u^{k}$
$\displaystyle\displaystyle=\sum_{k=0}^{\infty}u^{k}\sum_{i_{1}+2i_{2}+3i_{3}+\cdots+ki_{k}=k}\prod_{j=1}^{k}\frac{t_{j}^{i_{j}}}{i_{j}!\,(2j)^{i_{j}}}$
$\displaystyle\displaystyle=\sum_{k=0}^{\infty}\sum_{i_{1}+2i_{2}+3i_{3}+\cdots+ki_{k}=k}\prod_{j=1}^{k}\frac{t_{j}^{i_{j}}u^{ji_{j}}}{i_{j}!\,(2j)^{i_{j}}}$
$\displaystyle\displaystyle=\prod_{j=1}^{\infty}\bigg{(}\sum_{i_{j}=0}^{\infty}\frac{t_{j}^{i_{j}}u^{ji_{j}}}{i_{j}!\,(2j)^{i_{j}}}\bigg{)}$
$\displaystyle\displaystyle=\prod_{j=1}^{\infty}\exp\Big{(}\frac{t_{j}u^{j}}{2j}\Big{)}=\exp\bigg{(}\sum_{j=1}^{\infty}\frac{t_{j}u^{j}}{2j}\bigg{)}.$

Now set $\displaystyle t_{1}=d^{1/2}$ and $\displaystyle t_{j}=1$ for all
$\displaystyle j\geq 2$. Then $\displaystyle a_{k}(t)$ reduces to
$\displaystyle a_{k}$; $\displaystyle G_{t}(u)$ reduces to $\displaystyle
G(u)$, the formal generating-function for the sequence
$\displaystyle\\{a_{k}\\}$; and we obtain

$\displaystyle\displaystyle G(u)=\sum_{k=0}^{\infty}a_{k}u^{k}$
$\displaystyle\displaystyle=\exp\bigg{(}\tfrac{1}{2}d^{1/2}u+\tfrac{1}{2}\sum_{j=2}^{\infty}\frac{u^{j}}{j}\bigg{)}$
$\displaystyle\displaystyle=\exp\bigg{(}\tfrac{1}{2}(d^{1/2}-1)u+\tfrac{1}{2}\sum_{j=1}^{\infty}\frac{u^{j}}{j}\bigg{)}$
$\displaystyle\displaystyle=\exp\big{(}\tfrac{1}{2}(d^{1/2}-1)u-\tfrac{1}{2}\log(1-u)\big{)}$
$\displaystyle\displaystyle=\exp\big{(}(d^{1/2}-1)u/2\big{)}\cdot(1-u)^{-1/2}.$

Expanding both of the latter functions in infinite series, we obtain

$\displaystyle\displaystyle G(u)$
$\displaystyle\displaystyle=\bigg{(}\sum_{l=0}^{\infty}\frac{\big{(}(d^{1/2}-1)/2\big{)}^{l}}{l!}u^{l}\bigg{)}\bigg{(}\sum_{m=0}^{\infty}\frac{(1/2)_{m}}{m!}u^{m}\bigg{)}$
$\displaystyle\displaystyle=\sum_{k=0}^{\infty}u^{k}\sum_{l=0}^{k}\frac{\big{(}(d^{1/2}-1)/2\big{)}^{l}}{l!}\frac{(1/2)_{k-l}}{(k-l)!}.$

By comparing the coefficients of $\displaystyle u^{k}$ in $\displaystyle
G_{t}(u)$ and $\displaystyle G(u)$, we obtain the equality in (A.31).

Since
$\displaystyle\big{(}(d^{1/2}-1)/2\big{)}^{l}\leq\big{(}(d^{1/2}-1)/2\big{)}_{l}$
then, by applying the well-known convolution identity,

$\displaystyle\sum_{l=0}^{k}\frac{(v_{1})_{l}}{l!}\frac{(v_{2})_{k-l}}{(k-l)!}=\frac{(v_{1}+v_{2})_{k}}{k!},$

$\displaystyle v_{1},v_{2}\in\mathbb{R}$, we obtain

$\displaystyle
a_{k}\leq\sum_{l=0}^{k}\frac{\big{(}(d^{1/2}-1)/2\big{)}_{l}}{l!}\frac{(1/2)_{k-l}}{(k-l)!}=\frac{(d^{1/2}/2)_{k}}{k!}.$

The proof of (A.31) now is complete. $\displaystyle\Box$

Next, we establish (A.12).

######  Lemma A.2.

Let $\displaystyle\gamma_{1}=\tfrac{1}{2}(\sqrt{3}+1)$. Then, for all
$\displaystyle k\geq 1$ and $\displaystyle d\geq 1$,

$\frac{(d^{1/2}/2)_{k}}{(d/2)_{k}}\leq\gamma_{1}^{k-1}((k-1)!)^{1/2}d^{-k/2}.$
(A.32)

Proof. Denote by $\displaystyle r_{k}$ the left-hand side of (A.32). Then
$\displaystyle r_{1}=d^{-1/2}$ and, for $\displaystyle j\geq 1$,

$\displaystyle\frac{r_{j+1}}{r_{j}}=\frac{(d^{1/2}/2)_{j+1}}{(d/2)_{j+1}}\cdot\frac{(d/2)_{j}}{(d^{1/2}/2)_{j}}=\frac{d^{1/2}+2j}{d+2j}.$

We claim that

$\frac{d^{1/2}+2j}{d+2j}\leq\gamma_{1}j^{1/2}d^{-1/2},$ (A.33)

$\displaystyle j\geq 1$, and we prove this as follows. Define $\displaystyle
y\equiv y(t)=(t^{2}+2jt)/(t^{2}+2j)$, $\displaystyle t\geq 1$. It is simple to
verify that $\displaystyle y(t)$ attains its maximum at $\displaystyle
t_{j}:=1+(2j+1)^{1/2}$ and that

$\displaystyle y(t_{j})=\frac{1}{2}\big{(}1+(2j+1)^{1/2}\big{)}.$

For $\displaystyle j\geq 1$, it is also straightforward to show that

$\displaystyle\frac{1}{2}\big{(}1+(2j+1)^{1/2}\big{)}\leq\gamma_{1}j^{1/2}.$

Consequently,

$\displaystyle\frac{d+2jd^{1/2}}{d+2j}=y(d^{1/2})\leq
y(t_{j})\leq\gamma_{1}j^{1/2}.$

Multiplying the left- and right-hand sides of the above inequality by
$\displaystyle d^{-1/2}$, we obtain (A.33). $\displaystyle\Box$