# High-dimensional copula-based Wasserstein dependence.

Steven De Keyser  Irène Gijbels  Department of Mathematics, KU Leuven,
Celestijnenlaan 200B, B-3001 Leuven (Heverlee), Belgium

###### Abstract

We generalize $2$-Wasserstein dependence coefficients to measure dependence
between a finite number of random vectors. This generalization includes
theoretical properties, and in particular focuses on an interpretation of
maximal dependence and an asymptotic normality result for a proposed semi-
parametric estimator under a Gaussian copula assumption. In addition, we
discuss general axioms for dependence measures between multiple random
vectors, other plausible normalizations, and various examples. Afterwards, we
look into plug-in estimators based on penalized empirical covariance matrices
in order to deal with high dimensionality issues and take possible marginal
independencies into account by inducing (block) sparsity. The latter ideas are
investigated via a simulation study, considering other dependence coefficients
as well. We illustrate the use of the developed methods in two real data
applications.

###### keywords:

Copula , Normal scores rank correlation , Regularization , Sparsity ,
Wasserstein dependence
2020 MSC: Primary 62Axx, 62Hxx; Secondary 62Exx, 62Gxx.

##  1 Introduction

A prominent line of research in statistics considers measuring dependence
between groups of variables. In case of two groups, greatly celebrated is the
canonical correlation analysis of [27] relying on the Pearson correlation
coefficient. To step away from the assumption of Gaussianity, concordance
measures as studied in, e.g., [22, 37, 42] among many others, are also used
for quantifying general monotonic associations between two random vectors in
[23]. In [26], measures of association computed from collapsed random
variables are used to measure the dependence between random vectors.
Fundamental is copula theory (e.g., [38, 43]), allowing to split multivariate
distributions into marginal distributions on the one hand, and a dependence
structure described by the copula on the other hand. Especially when the
marginals are continuous, the preference often goes to copula-based dependence
measures since then, by Sklar’s theorem [43], the copula is unique, and hence
margin-free dependence can unequivocally be quantified.

Statistical independence between random vectors holds if and only if the true
underlying copula is the product of the marginal copulas (where a one
dimensional copula is basically a uniform distribution on $[0,1]$), yielding
zero dependence. However, the dependence measures mentioned above do not
detect all types of departure from independence, meaning that they might
vanish while the independence product copula is misspecified. Since the work
of [44], there has been a growing interest for dependence measures that
completely characterize independence. Some recent developments are, e.g., the
Hoeffding’s Phi-Square measure of [35], the $\Phi$-dependence measures of [11]
(of which the Hellinger correlation [20] and essential dependence [51] are
particular cases), the coefficient of Chatterjee [8, 2, 18, 1], and the
$2$-Wasserstein coefficients of [36].

The aim of this article is to elaborate further on the optimal transport
measures of [36]. First, the focus will be on extending their dependence
coefficients from two to finitely many random vectors. We do this from a
copula point of view. This includes generalizing the results of [36], and also
verifying the axioms stated in [11] (see also A). Some additional examples,
focusing on, e.g., the impact of the normalization, are provided as well.
Afterwards, we dive into the Bures-Wasserstein dependence making a Gaussian
copula assumption. This yields dependence measures that are attractive, and
more amenable for estimation. The results are a far from straightforward
extension of results of [36] to the case of a finite number of random vectors,
and require significant mathematical care.

The proposed semi-parametric estimation approach of the Bures-Wasserstein
coefficients relies on the sample matrix of normal scores rank correlations
(see, e.g., [24]). Since we extend the theory to a general finite amount of
groups of variables, high dimensional cases with a large number of variables
compared to the sample size are of study interest as well. Acclaimed
penalization techniques are known to significantly improve (inverse)
covariance matrix estimation (see, e.g., [29] and references therein). We
utilize these ideas in our Gaussian copula context in order to correct for
high dimensionality bias and possibly enforce sparsity at the individual level
or group level to take plausible marginal independencies into account.

The outline of this paper is the following. Section 2 explains how optimal
transport theory is combined with copula theory in order to arrive at a
dependence measure between multiple groups of random variables that completely
characterizes independence and is invariant with respect to the univariate
marginal distributions. The verification of the properties postulated in [11]
for such dependence measures is also part of this section. The Gaussian copula
approach is discussed in Section 3, with special attention to the meaning of
maximal dependence, and asymptotic normality of the proposed semi-parametric
estimator. Afterwards, different regularization techniques for Gaussian copula
covariance matrices are discussed in Section 4. Next to an empirical
illustration of the asymptotic normality result, simulations are performed in
Section 5 to investigate the performance of these regularization techniques on
various dependence coefficients for random vectors. Two real data applications
are discussed in Section 6. All proofs are deferred to the Appendix. Any
experiments reported can be reproduced via the source code available at
https://github.com/StevenDeKeyser98/VecDep. Additional figures are included in
the Supplementary Material.

##  2 General $2$-Wasserstein dependence

We consider a $q$-dimensional random vector
$\mathbf{X}=(\mathbf{X}_{1},\dots,\mathbf{X}_{k})$ defined on
$(\mathbb{R}^{q},\mathcal{B}(\mathbb{R}^{q}),\lambda^{q})$ consisting of $k$
marginal random vectors $\mathbf{X}_{i}=(X_{i1},\dots,X_{id_{i}})$ for
$i=1,\dots,k$ having $d_{i}$ continuous univariate marginal random variables
$X_{ij}$ for $j=1,\dots,d_{i}$. The numbers
$d_{1},\dots,d_{k}\in\mathbb{Z}_{>0}$ are such that $q=d_{1}+\cdots+d_{k}$,
and $\lambda^{q}$ denotes the $q$-dimensional Lebesgue measure defined on
$\mathcal{B}(\mathbb{R}^{q})$, the Borel sigma-algebra on $\mathbb{R}^{q}$.
Let $\mathbb{P}:\mathcal{B}(\mathbb{R}^{q})\rightarrow\mathbb{R}$ be a
probability measure. Our aim is to measure the dependence between
$\mathbf{X}_{1},\dots,\mathbf{X}_{k}$. For $\mathbb{I}=[0,1]$ and
$\mathcal{P}(\mathbb{I}^{q})$ the set of Borel probability measures on
$\mathbb{I}^{q}$, the random vector $\mathbf{X}$ is assigned a copula
probability law $\mu_{C}\in\mathcal{P}(\mathbb{I}^{q})$ having corresponding
marginal copula laws $\mu_{C_{i}}\in\mathcal{P}(\mathbb{I}^{d_{i}})$ of
$\mathbf{X}_{i}$ for $i=1,\dots,k$. Note that in case $d_{i}=1$, the law
$\mu_{C_{i}}$ is that of a uniform distribution on $\mathbb{I}$. We denote
$\Gamma(\mu_{C_{1}},\dots,\mu_{C_{k}})$ for the set of measures
$\gamma\in\mathcal{P}(\mathbb{I}^{q})$ such that

$\mu_{C_{i}}(B_{i})=\gamma\left(\mathbb{I}^{d_{1}}\times\cdots\times\mathbb{I}^{d_{i-1}}\times
B_{i}\times\mathbb{I}^{d_{i+1}}\times\cdots\times\mathbb{I}^{d_{k}}\right)$

for all $B_{i}\in\mathcal{B}(\mathbb{I}^{d_{i}})$ and $i=1,\dots,k$. This is a
natural generalization of the set of coupling measures. Obviously,
$\mu_{C}\in\Gamma(\mu_{C_{1}},\dots,\mu_{C_{k}})$. Quantifying the intensity
of relation between $\mathbf{X}_{1},\dots,\mathbf{X}_{k}$ can be done by
measuring the difference between $\mu_{C}$ and the product measure
$\mu_{C_{1}}\times\cdots\times\mu_{C_{k}}$. We use the $2$-Wasserstein
distance, whose square is, for certain measures
$\pi,\widetilde{\pi}\in\mathcal{P}(\mathbb{I}^{q})$ given by

$W_{2}^{2}(\pi,\widetilde{\pi})=\underset{\gamma\in\Gamma(\pi,\widetilde{\pi})}{\inf}\int_{\mathbb{I}^{2q}}||\mathbf{u}-\mathbf{v}||^{2}d\gamma(\mathbf{u},\mathbf{v})=\underset{\begin{subarray}{c}\mathbf{U}\sim\pi\\\
\mathbf{V}\sim\widetilde{\pi}\end{subarray}}{\inf}\mathbb{E}\left(||\mathbf{U}-\mathbf{V}||^{2}\right).$
(1)

The interpretation of the metric (1) is optimal transport, see, e.g., [39]. It
is the minimal effort (cost) required to transform the mass of $\pi$ into the
mass of $\widetilde{\pi}$, i.e., for every $(\mathbf{u},\mathbf{v})$ transport
an infinitesimal amount of mass $d\gamma(\mathbf{u},\mathbf{v})$ from
$\mathbf{u}$ to $\mathbf{v}$ at a distance cost of
$||\mathbf{u}-\mathbf{v}||^{2}$. Aggregating the mass
$\gamma(\\{\mathbf{u}\\}\times\mathbb{I}^{q})$ that leaves $\mathbf{u}$ gives
$\pi(\mathbf{u})$ and the total mass
$\gamma(\mathbb{I}^{q}\times\\{\mathbf{v}\\})$ that reaches $\mathbf{v}$
equals $\widetilde{\pi}(\mathbf{v})$. For certain non-degenerate (i.e. not
Dirac delta distributions) reference laws
$\nu_{i}\in\mathcal{P}(\mathbb{I}^{d_{i}})$ for $i=1,\dots,k$, we now define

$\begin{split}T_{d_{1},\dots,d_{k}}(\mu_{C};\nu_{1},\dots,\nu_{k})&=W_{2}^{2}(\mu_{C},\nu_{1}\times\cdots\times\nu_{k})-W_{2}^{2}(\mu_{C_{1}}\times\cdots\times\mu_{C_{k}},\nu_{1}\times\cdots\times\nu_{k})\\\
&=W_{2}^{2}(\mu_{C},\nu_{1}\times\cdots\times\nu_{k})-\sum_{i=1}^{k}W_{2}^{2}(\mu_{C_{i}},\nu_{i}),\end{split}$

where the second identity is known to be true (see, e.g., [39]). We call a
subset $A\subset\mathcal{P}(\mathbb{I}^{q})$ $W_{2}$-compact if every sequence
in the metric space $(A,W_{2})$ has a convergent subsequence with limit in
$A$. Lemma 1 gives the main properties of $T_{d_{1},\dots,d_{k}}$. Proofs of
Lemma 1, and all other theoretical results of Section 2, can be found in B.

######  Lemma 1.

It holds that

  * (a)

$T_{d_{1},\dots,d_{k}}(\mu_{C};\nu_{1},\dots,\nu_{k})\geq 0$

  * (b)

$T_{d_{1},\dots,d_{k}}(\mu_{C_{1}}\times\cdots\times\mu_{C_{k}};\nu_{1},\dots,\nu_{k})=0$

  * (c)

If either $\mu_{C_{i}}=\nu_{i}$ for all $i=1,\dots,k$ or $\nu_{i}$ is
absolutely continuous (with respect to $\lambda^{d_{i}}$) for all
$i=1,\dots,k$, then $T_{d_{1},\dots,d_{k}}(\mu_{C};\nu_{1},\dots,\nu_{k})=0$
implies $\mu_{C}=\mu_{C_{1}}\times\cdots\times\mu_{C_{k}}$

  * (d)

The set $\Gamma(\mu_{C_{1}},\dots,\mu_{C_{k}})$ is $W_{2}$-compact in
$\mathcal{P}(\mathbb{I}^{q})$ and the mapping $T_{d_{1},\dots,d_{k}}(\hskip
2.84544pt\cdot\hskip
2.84544pt;\nu_{1},\dots,\nu_{k}):(\Gamma(\mu_{C_{1}},\dots,\mu_{C_{k}}),W_{2})\rightarrow(\mathbb{R},|\cdot|)$
is continuous.

Its interpretation, together with its mathematical properties, make
$T_{d_{1},\dots,d_{k}}(\mu_{C};\nu_{1},\dots,\nu_{k})$ an appealing measure of
dependence between $k$ random vectors $\mathbf{X}_{1},\dots,\mathbf{X}_{k}$.
In what follows, we assume that $\nu_{i}$ is absolutely continuous for
$i=1,\dots,k$, and let $G_{C}\subset\Gamma(\mu_{C_{1}},\dots,\mu_{C_{k}})$ be
a compact set such that $\mu_{C}\in G_{C}$. General axioms for dependence
measures between multiple random vectors are formulated in [11], see also A.
Lemma 1 offers aid in proving them here, see Proposition 1.

######  Proposition 1.

Consider Axioms (A1)-(A8) given in A, and a normalized version of
$T_{d_{1},\dots,d_{k}}(\mu_{C};\nu_{1},\dots,\nu_{k})$ given by

$\mathcal{D}(\mu_{C};\nu_{1},\dots,\nu_{k})=\frac{T_{d_{1},\dots,d_{k}}(\mu_{C};\nu_{1},\dots,\nu_{k})}{\sup_{\pi\in
G_{C}}T_{d_{1},\dots,d_{k}}(\pi;\nu_{1},\dots,\nu_{k})}.$ (2)

Then, $\mathcal{D}$ satisfies (A1)-(A3) and (A5)-(A7). Axioms (A4) and (A8)
are satisfied by the non-normalized version $T_{d_{1},\dots,d_{k}}$.

The supremum in (2) is attained when $G_{C}$ is $W_{2}$-compact (e.g.,
$G_{C}=\Gamma(\mu_{C_{1}},\dots,\mu_{C_{k}})$) because of (d) in Lemma 1. It
represents the case of maximal dependence, which characterization (and hence
the overall behaviour of (2)), largely depends on $G_{C}$. There is still
freedom in choosing the normalization by picking the set $G_{C}$. It might
impose additional constraints (in addition to having marginals $\mu_{C_{i}}$)
on the $\pi\in\Gamma(\mu_{C_{1}},\dots,\mu_{C_{k}})$ that characterizes
maximal dependence (for example $\pi$ should be in the same copula family as
$\mu_{C}$). Strictly speaking, we can have a different $G_{C}$ for every
copula $C$, even if the marginals are the same.

Regarding axioms (A4) and (A8), we make the following remark.

######  Remark 1.

While (2) does not satisfy axiom (A4) in general (that is for every possible
choice of $G_{C}$), there might still be some specific choices for $G_{C}$
such that (A4) is satisfied. We illustrate this further in Example 1.

Also, when considering $C_{n}\to C$ uniformly as $n\to\infty$, axiom (A8) can
be satisfied by (2) under some extra constraints. The numerator converges if
$G_{C}$ is chosen such that $C\mapsto\sup_{\pi\in
G_{C}}T_{d_{1},\dots,d_{k}}(\pi;\nu_{1},\dots,\nu_{k})$ is continuous
(considering the uniform metric on the space of copulas). In Section 3, we see
that this holds in the class of Gaussian copulas when taking
$G_{C}=\Gamma(\mu_{C_{1}},\dots,\mu_{C_{k}})$.

We now arrive at a natural generalization of the Wasserstein dependence
coefficients of [36], which come from (2) with two particular choices of
reference measures.

######  Definition 1.

For $m\in\mathbb{Z}_{>0}$, let $\gamma_{m}$ denote the measure of an
$m$-variate Gaussian copula with identity correlation matrix (i.e., the
$m$-variate independence copula). For $q=d_{1}+\cdots+d_{k}$ with
$d_{1},\dots,d_{k}\in\mathbb{Z}_{>0}$ and for $\mu_{C}\in
G_{C}\subset\Gamma(\mu_{C_{1}},\dots,\mu_{C_{k}})$, where $G_{C}$ is
$W_{2}$-compact and non-degenerate
$\mu_{C_{i}}\in\mathcal{P}(\mathbb{I}^{d_{i}})$ for $i=1,\dots,k$, define

$\begin{split}\mathcal{D}_{1}^{d_{1},\dots,d_{k}}(\mu_{C})&=\mathcal{D}(\mu_{C};\gamma_{d_{1}},\dots,\gamma_{d_{k}})=\frac{W_{2}^{2}(\mu_{C},\gamma_{q})-\sum_{i=1}^{k}W_{2}^{2}(\mu_{C_{i}},\gamma_{d_{i}})}{\sup_{\pi\in
G_{C}}W_{2}^{2}(\pi,\gamma_{q})-\sum_{i=1}^{k}W_{2}^{2}(\mu_{C_{i}},\gamma_{d_{i}})}\\\
\mathcal{D}_{2}^{d_{1},\dots,d_{k}}(\mu_{C})&=\mathcal{D}(\mu_{C};\mu_{C_{1}},\dots,\mu_{C_{k}})=\frac{W_{2}^{2}(\mu_{C},\mu_{C_{1}}\times\cdots\times\mu_{C_{k}})}{\sup_{\pi\in
G_{C}}W_{2}^{2}(\pi,\mu_{C_{1}}\times\cdots\times\mu_{C_{k}})}.\end{split}$

If the context is clear, we just write $\mathcal{D}_{r}(\mu_{C})$ for
$r\in\\{1,2\\}$, or also
$\mathcal{D}_{r}(\mathbf{X}_{1};\cdots;\mathbf{X}_{k})$ to emphasize that we
are measuring the dependence between $k$ random vectors (having joint copula
$C$).

Let us consider an example illustrating that $\mathcal{D}$ does not
necessarily satisfy Axiom (A4) in general.

######  Example 1.

Consider a random vector $(X_{1},X_{2},X_{3})$ having a trivariate Gaussian
copula $\mu_{C}$ with correlation matrix

$\mathbf{R}=\begin{pmatrix}1&\rho&0\\\ \rho&1&0\\\ 0&0&1\end{pmatrix},\hskip
11.38092pt\text{where}\hskip 5.69046pt-1\leq\rho\leq 1.$

Let $\mu_{C_{i}}$ be the marginal copula measure of $X_{i}$ for $i=1,2,3$,
which in this case is in fact the Lebesgue measure corresponding to a
$\mathcal{U}[0,1]$ distribution. The product measure
$\mu_{C_{1}}\times\mu_{C_{2}}\times\mu_{C_{3}}$ is the three dimensional
independence copula, being equal to the trivariate Gaussian copula with
identity correlation matrix $\mathbf{I}_{3}$. Also note that $X_{3}$ is
independent of $(X_{1},X_{2})$. One can quickly check that (using (3), see
further)

$\mathcal{D}_{2}(X_{1};X_{2};X_{3})=\frac{W_{2}^{2}(\mu_{C};\mu_{C_{1}}\times\mu_{C_{2}}\times\mu_{C_{3}})}{\sup_{\pi\in
G_{C}}W_{2}^{2}(\pi,\mu_{C_{1}}\times\mu_{C_{2}}\times\mu_{C_{3}})}=\frac{4-2\sqrt{1-\rho}-2\sqrt{1+\rho}}{\text{sup}_{\pi\in
G_{C}}W_{2}^{2}(\pi,\mu_{C_{1}}\times\mu_{C_{2}}\times\mu_{C_{3}})}.$

The remaining question is what to pick for $G_{C}$, i.e., which quantity do we
put in the denominator and defines the maximal amount of dependence. Well, if
$G_{C}=\Gamma(\mu_{C_{1}},\mu_{C_{2}},\mu_{C_{3}})$, we should find the
squared $2$-Wasserstein distance between the independence copula and any other
trivariate distribution having $\mu_{C_{1}},\mu_{C_{2}}$ and $\mu_{C_{3}}$ as
marginals, that is every possible trivariate copula. This is, as far as we are
concerned, an open problem. However, in this context, it is reasonable to
restrict $G_{C}$ to the Gaussian copula family. Doing so, one has

$\text{sup}_{\pi\in
G_{C}}W_{2}^{2}(\pi,\mu_{C_{1}}\times\mu_{C_{2}}\times\mu_{C_{3}})=W_{2}^{2}(\mu_{C^{\text{co}}},\mu_{C_{1}}\times\mu_{C_{2}}\times\mu_{C_{3}})=6-2\sqrt{3},$

where $\mu_{C^{\text{co}}}$ stands for the comonotonicity copula measure,
i.e., the limit of an equicorrelated
($\mathbf{R}=(1-\rho)\mathbf{I}_{3}+\rho\mathbf{1}_{3}\mathbf{1}_{3}^{\text{T}}$,
where $\mathbf{1}_{3}^{\text{T}}=(1,1,1)$ and $\rho\in(-1/2,1)$) Gaussian
copula with correlation $\rho$ tending to $1$. With this choice of $G_{C}$,
$(X_{1},X_{2},X_{3})$ can never reach the maximum dependence, since

$4-2\sqrt{1-\rho}-2\sqrt{1+\rho}\leq 4-2\sqrt{2}<6-2\sqrt{3}.$

Another way to put it, is that

$\mathcal{D}_{2}(X_{1};X_{2};X_{3})=\frac{2-\sqrt{1-\rho}-\sqrt{1+\rho}}{3-\sqrt{3}}<\frac{2-\sqrt{1-\rho}-\sqrt{1+\rho}}{2-\sqrt{2}}=\mathcal{D}_{2}(X_{1};X_{2}),$

where $\mathcal{D}_{2}(X_{1};X_{2})$ is computed in a similar way, also
restricting the couplings $\Gamma(\mu_{C_{1}},\mu_{C_{2}})$ to Gaussian ones.
We thus see that, when adding an independent component $X_{3}$ into
consideration, the dependence has decreased and hence axiom (A4) is definitely
not fulfilled. Taking a look at $\mathbf{R}$, it is maybe more tempting to
have maximal dependence when $|\rho|=1$ and restrict $G_{C}$ further to only
those $\pi\in\Gamma(\mu_{C_{1}},\mu_{C_{2}},\mu_{C_{3}})$ that are Gaussian
and furthermore satisfy $\pi(B_{1}\times B_{2}\times B_{3})=\pi(B_{1}\times
B_{2}\times\mathbb{I})\cdot\pi(\mathbb{I}\times\mathbb{I}\times B_{3})$ for
all $B_{1},B_{2},B_{3}\in\mathcal{B}(\mathbb{I})$. Then, it is quickly seen
that

$\text{sup}_{\pi\in
G_{C}}W_{2}^{2}(\pi,\mu_{C_{1}}\times\mu_{C_{2}}\times\mu_{C_{3}})=2-\sqrt{2},$

and hence $\mathcal{D}_{2}(X_{1};X_{2};X_{3})=\mathcal{D}_{2}(X_{1};X_{2})$,
in harmony with axiom (A4). So, for actual computation, it is better to
restrict $G_{C}$ to the Gaussian copula family, and if some additional
information is given (like zeroes in the correlation matrix), incorporating
this in $G_{C}$ can lead to a more interpretative dependence quantification.

Except for some families like normal distributions, computing the Wasserstein
distance is very involved and tools and theory for statistical inference are
still scarce. The authors of [36] give an overview of the literature so far,
concluding that additional theory is still needed, and propose a quasi-
Gaussian (based on covariance matrices) approach instead. We assume that the
copula of $\mathbf{X}$ is Gaussian.

##  3 A Gaussian copula approach

In this section, we assume a Gaussian copula model for $\mathbf{X}$, elaborate
more on maximal dependence, and discuss statistical inference within this
framework.

###  3.1 The Bures-Wasserstein distance

The main incentive is the well-known formula for the squared $2$-Wasserstein
distance between Gaussian distributions, say with covariance matrices
$\mathbf{R}$ and $\mathbf{S}$, becoming the so-called squared Bures-
Wasserstein distance (see, e.g., [45]) between $\mathbf{R}$ and $\mathbf{S}$:

$d_{W}^{2}(\mathbf{R},\mathbf{S})=\text{tr}(\mathbf{R})+\text{tr}(\mathbf{S})-2\text{tr}\left\\{\left(\mathbf{R}^{1/2}\mathbf{S}\mathbf{R}^{1/2}\right)^{1/2}\right\\},$
(3)

were tr stands for the trace of a matrix. We denote by
$\mathbb{S}^{q}=\\{\mathbf{R}\in\mathbb{R}^{q\times
q}:\mathbf{R}^{\text{T}}=\mathbf{R}\\}$ the set of symmetric $q\times q$
matrices, $\mathbb{S}_{\geq}^{q}\subset\mathbb{S}^{q}$ the set of positive
semi-definite ones and $\mathbb{S}^{q}_{>}\subset\mathbb{S}_{\geq}^{q}$ the
set of positive definite ones. Let again $q=d_{1}+\cdots+d_{k}$ and consider
$\mathbf{R}_{ii}\in\mathbb{S}^{d_{i}}_{\geq}$ for $i=1,\dots,k$. We also
define the set

$\Gamma(\mathbf{R}_{11},\dots,\mathbf{R}_{kk})=\left\\{\mathbf{A}\in\mathbb{S}^{q}_{\geq}:\mathbf{A}=\begin{pmatrix}\mathbf{R}_{11}&\boldsymbol{\Psi}_{12}&\cdots&\boldsymbol{\Psi}_{1k}\\\
\boldsymbol{\Psi}_{12}^{\text{T}}&\mathbf{R}_{22}&\cdots&\boldsymbol{\Psi}_{2k}\\\
\vdots&\vdots&\ddots&\vdots\\\
\boldsymbol{\Psi}_{1k}^{\text{T}}&\boldsymbol{\Psi}_{2k}^{\text{T}}&\cdots&\mathbf{R}_{kk}\end{pmatrix}\hskip
5.69046pt\text{for some}\hskip
5.69046pt\boldsymbol{\Psi}_{ij}\in\mathbb{R}^{d_{i}\times d_{j}}\right\\},$

as the set of all covariance matrices of random vectors
$\mathbf{Z}=(\mathbf{Z}_{1},\dots,\mathbf{Z}_{k})$ such that the covariance
matrix of $\mathbf{Z}_{i}$, being $\mathbf{R}_{ii}$, remains fixed for all
$i=1,\dots,k$ and

$\mathbf{R}_{0}=\begin{pmatrix}\mathbf{R}_{11}&\mathbf{0}_{12}&\dots&\mathbf{0}_{1k}\\\
\mathbf{0}_{12}^{\text{T}}&\mathbf{R}_{22}&\dots&\mathbf{0}_{2k}\\\
\vdots&\vdots&\ddots&\vdots\\\
\mathbf{0}_{1k}^{\text{T}}&\mathbf{0}_{2k}^{\text{T}}&\dots&\mathbf{R}_{kk}\end{pmatrix},$
(4)

with $\mathbf{0}_{ij}\in\mathbb{R}^{d_{i}\times d_{j}}$ a matrix of zeroes, as
the covariance matrix when the $\mathbf{Z}_{i}$ are all independent.

Consider now a random vector
$\mathbf{X}=(\mathbf{X}_{1},\dots,\mathbf{X}_{k})$ having a Gaussian copula
with covariance matrix

$\mathbf{R}=\begin{pmatrix}\mathbf{R}_{11}&\mathbf{R}_{12}&\dots&\mathbf{R}_{1k}\\\
\mathbf{R}_{12}^{\text{T}}&\mathbf{R}_{22}&\dots&\mathbf{R}_{2k}\\\
\vdots&\vdots&\ddots&\vdots\\\
\mathbf{R}_{1k}^{\text{T}}&\mathbf{R}_{2k}^{\text{T}}&\dots&\mathbf{R}_{kk}\end{pmatrix}.$
(5)

This means that (5) is the usual covariance matrix of the random vector
$\mathbf{Z}=(\mathbf{Z}_{1},\dots,\mathbf{Z}_{k})$, with
$\mathbf{Z}_{i}=(Z_{i1},\dots,Z_{id_{i}})$ and $Z_{ij}=(\Phi^{-1}\circ
F_{ij})(X_{ij})$ for $i=1,\dots,k$ and $j=1,\dots,d_{i}$, where $F_{ij}$ is
the marginal distribution of $X_{ij}$ and $\Phi^{-1}$ the univariate standard
normal quantile function. Measuring the dependence between
$\mathbf{X}_{1},\dots,\mathbf{X}_{k}$ can be done by utilizing the
$2$-Wasserstein dependence coefficients of Definition 1, now taking $d_{W}$
for $W_{2}$.

######  Definition 2.

For $q=d_{1}+\cdots+d_{k}$ with $d_{1},\dots,d_{k}\in\mathbb{Z}_{>0}$,
$\mathbf{R}\in\Gamma(\mathbf{R}_{11},\dots,\mathbf{R}_{kk})$ with
$\mathbf{R}_{ii}\in\mathbb{S}^{d_{i}}_{\geq}\setminus\\{\mathbf{0}\\}$ for
$i=1,\dots,k$ and a $d_{W}$-compact set
$G_{\mathbf{R}}\subset\Gamma(\mathbf{R}_{11},\dots,\mathbf{R}_{kk})$ with
$\mathbf{R}\in G_{\mathbf{R}}$, define

$\begin{split}\mathcal{D}_{1}(\mathbf{R};d_{1},\dots,d_{k})&=\frac{d_{W}^{2}(\mathbf{R},\mathbf{I}_{q})-\sum_{i=1}^{k}d_{W}^{2}(\mathbf{R}_{ii},\mathbf{I}_{d_{i}})}{\text{sup}_{\mathbf{A}\in
G_{\mathbf{R}}}d_{W}^{2}(\mathbf{A},\mathbf{I}_{q})-\sum_{i=1}^{k}d_{W}^{2}(\mathbf{R}_{ii},\mathbf{I}_{d_{i}})}\\\
\mathcal{D}_{2}(\mathbf{R};d_{1},\dots,d_{k})&=\frac{d_{W}^{2}(\mathbf{R},\mathbf{R}_{0})}{\sup_{\mathbf{A}\in
G_{\mathbf{R}}}d_{W}^{2}(\mathbf{A},\mathbf{R}_{0})},\end{split}$

where $\mathbf{R}_{0}$ is the matrix given in (4). If the context is clear, we
also write $\mathcal{D}_{r}(\mathbf{R})$ or
$\mathcal{D}_{r}(\mathbf{X}_{1};\cdots;\mathbf{X}_{k})$ for $r\in\\{1,2\\}$.

If the true copula is indeed Gaussian, the adequacy of the Bures-Wasserstein
dependence remains, and we obtain something way more easy to handle. In order
to make them fully practically usable, that is to say suitably attractive for
estimation, we ought to find explicit expressions for the suprema in the
denominator of the dependence measures. When
$G_{\mathbf{R}}=\Gamma(\mathbf{R}_{11},\dots,\mathbf{R}_{kk})$ and $k=2$, the
authors of [36] found elegant solutions to this problem, which we generalize
to general $k$.

###  3.2 Maximal Bures-Wasserstein dependence

We need the definition of majorization of two vectors and its behaviour under
convex functions, as studied in [33].

######  Definition 3.

For two vectors $\mathbf{x},\mathbf{y}\in\mathbb{R}^{q}$, we say that
$\mathbf{y}$ majorizes $\mathbf{x}$, denoted as $\mathbf{x}\prec\mathbf{y}$,
if

$\begin{cases}\sum_{i=1}^{\ell}x_{[i]}\leq\sum_{i=1}^{\ell}y_{[i]}&\mbox{for
}\ell=1,\dots,q-1\\\ \sum_{i=1}^{q}x_{[i]}=\sum_{i=1}^{q}y_{[i]}\end{cases},$

where $x_{[1]}\geq\dots\geq x_{[q]}$ are the components of $\mathbf{x}$ in
decreasing order, and similarly for $\mathbf{y}$.

When $\boldsymbol{\lambda}$ and $\boldsymbol{\mu}$ are the vectors of
eigenvalues of two correlation matrices, $\boldsymbol{\lambda}$ being
majorized by $\boldsymbol{\mu}$ means that the proportion of the total
variance explained by the $\ell$ first principal components is larger for the
correlation matrix with eigenvalues $\boldsymbol{\mu}$, for any
$\ell\in\\{1,\dots,q-1\\}$, than for the correlation matrix with eigenvalues
$\boldsymbol{\lambda}$. Fixing $k$ covariance matrices
$\mathbf{R}_{11},\dots,\mathbf{R}_{kk}$ with
$\mathbf{R}_{ii}\in\mathbb{S}^{d_{i}}_{\geq}$, the goal is now to find
$\mathbf{R}_{m}\in\Gamma(\mathbf{R}_{11},\dots,\mathbf{R}_{kk})$ whose ordered
eigenvalues majorize those of any matrix
$\mathbf{A}\in\Gamma(\mathbf{R}_{11},\dots,\mathbf{R}_{kk})$. Together with
Lemma 2 (see Proposition 3.C.1 in [33]), this will enable us to characterize
maximal dependence between $k$ random vectors in terms of covariance matrices.

######  Lemma 2.

If $g:I\rightarrow\mathbb{R}$ is convex, with $I\subset\mathbb{R}$ an
interval, then

$\mathbf{x}\prec\mathbf{y}\hskip 5.69046pt\implies\hskip
5.69046pt\sum_{i=1}^{q}g(x_{i})\leq\sum_{i=1}^{q}g(y_{i})$

for all $\mathbf{x},\mathbf{y}\in I^{q}$.

Since the $2$-Wasserstein dependence coefficients satisfy axiom (A1), we can
assume that $d_{1}\leq d_{2}\leq\dots\leq d_{k}$ without loss of generality.
Suppose that

$\mathbf{R}_{ii}=\mathbf{U}_{ii}\boldsymbol{\Lambda}_{ii}\mathbf{U}_{ii}^{\text{T}}$
(6)

is the eigendecomposition of $\mathbf{R}_{ii}$, i.e., with
$\boldsymbol{\Lambda}_{ii}=\text{diag}(\lambda_{1,ii},\lambda_{2,ii},\dots,\lambda_{d_{i},ii})$
the $d_{i}\times d_{i}$ diagonal matrix with $d_{i}$ ordered eigenvalues
$\lambda_{1,ii}\geq\lambda_{2,ii}\geq\dots\geq\lambda_{d_{i},ii}$ on the
diagonal (counting multiplicities), and $\mathbf{U}_{ii}$ an orthogonal matrix
containing the corresponding eigenvectors for $i=1,\dots,k$. The proof of
Proposition 2, and all other theoretical results of Section 3, are provided in
C.

######  Proposition 2.

Let $\mathbf{R}_{ii}\in\mathbb{S}^{d_{i}}_{\geq}$ have eigendecomposition (6)
for $i=1,\dots,k$. Define the matrix $\mathbf{R}_{m}$ as

$\mathbf{R}_{m}=\begin{pmatrix}\mathbf{R}_{11}&\boldsymbol{\Psi}_{12}&\cdots&\boldsymbol{\Psi}_{1k}\\\
\boldsymbol{\Psi}_{12}^{\text{T}}&\mathbf{R}_{22}&\cdots&\boldsymbol{\Psi}_{2k}\\\
\vdots&\vdots&\ddots&\vdots\\\
\boldsymbol{\Psi}_{1k}^{\text{T}}&\boldsymbol{\Psi}_{2k}^{\text{T}}&\cdots&\mathbf{R}_{kk}\end{pmatrix}\in\mathbb{R}^{q\times
q},$ (7)

with $q=d_{1}+\cdots+d_{k}$ and $d_{1}\leq d_{2}\leq\dots\leq d_{k}$, and off-
diagonal blocks

$\boldsymbol{\Psi}_{ij}=\mathbf{U}_{ii}\boldsymbol{\Lambda}_{ii}^{1/2}\boldsymbol{\Pi}_{ij}\boldsymbol{\Lambda}_{jj}^{1/2}\mathbf{U}_{jj}^{\text{T}}\in\mathbb{R}^{d_{i}\times
d_{j}},$

where

$\boldsymbol{\Pi}_{ij}=\left(\hskip 2.84544pt\mathbf{I}_{d_{i}}\hskip
5.69046pt\mathbf{0}_{d_{i}\times(d_{j}-d_{i})}\hskip
2.84544pt\right)\in\mathbb{R}^{d_{i}\times d_{j}},$

the $d_{i}\times d_{j}$ upper left block of $\mathbf{I}_{d_{i}+d_{j}}$ for
$i=1,\dots,k,j=i+1,\dots,k$ (denoting $\mathbf{0}_{d_{i}\times(d_{j}-d_{i})}$
for the $d_{i}\times(d_{j}-d_{i})$ matrix of zeroes). If we define
$\lambda_{j,ii}=0$ for $j=d_{i}+1,\dots,q$, the eigenvalues of
$\mathbf{R}_{m}$ are

$\boldsymbol{\lambda}(\mathbf{R}_{m})=(\lambda_{j,11}+\lambda_{j,22}+\dots+\lambda_{j,kk})_{j=1}^{q}.$
(8)

Furthermore, for any
$\mathbf{A}\in\Gamma(\mathbf{R}_{11},\dots,\mathbf{R}_{kk})$ with eigenvalues
$\boldsymbol{\lambda}(\mathbf{A})=(\lambda_{j})_{j=1}^{q}$, it holds that

$\boldsymbol{\lambda}(\mathbf{A})\prec\boldsymbol{\lambda}(\mathbf{R}_{m}).$

Example 2 gives the matrix $\mathbf{R}_{m}$ for some specific cases.

######  Example 2.

Some expressions for $\mathbf{R}_{m}$ in case $k=2$ can be found in [36]. If
$d_{i}=1$ with $\mathbf{R}_{ii}=1$ for all $i=1,\dots,k$, the matrix
$\mathbf{R}_{m}$ is obviously given by $\mathbf{1}_{q\times q}$, a matrix full
of ones. Consider next $\mathbf{Z}_{i}=(Z_{i1},Z_{i2})$ for $i=1,\dots,k$,
i.e., $k$ bivariate random vectors, with covariance matrix of
$(Z_{i1},Z_{i2})$ given by

$\mathbf{R}_{ii}=\begin{pmatrix}1&\rho_{i}\\\ \rho_{i}&1\end{pmatrix}.$

Assuming $\rho_{i},\rho_{j}\geq 0$, one can check that
$\boldsymbol{\Psi}_{ij}$ of $\mathbf{R}_{m}$ in (7) for $i=1,\dots,k-1$ and
$j=i+1,\dots,k$ is given by

$\boldsymbol{\Psi}_{ij}=\begin{pmatrix}\frac{\sqrt{1+\rho_{i}}\sqrt{1+\rho_{j}}+\sqrt{1-\rho_{i}}\sqrt{1-\rho_{j}}}{2}&\frac{\sqrt{1+\rho_{i}}\sqrt{1+\rho_{j}}-\sqrt{1-\rho_{i}}\sqrt{1-\rho_{j}}}{2}\vspace{0.1cm}\\\
\frac{\sqrt{1+\rho_{i}}\sqrt{1+\rho_{j}}-\sqrt{1-\rho_{i}}\sqrt{1-\rho_{j}}}{2}&\frac{\sqrt{1+\rho_{i}}\sqrt{1+\rho_{j}}+\sqrt{1-\rho_{i}}\sqrt{1-\rho_{j}}}{2}\end{pmatrix}.$
(9)

The result is similar in case $\rho_{i}\leq 0$ or $\rho_{j}\leq 0$ up to some
signs (orthogonal transformations, to which $d_{W}$ is invariant). The
principal components of $(Z_{i1},Z_{i2})$ are

$Y_{i1}=\frac{1}{\sqrt{2}}(Z_{i1}+Z_{i2})\hskip 5.69046pt\text{and}\hskip
5.69046ptY_{i2}=\frac{1}{\sqrt{2}}(Z_{i2}-Z_{i1}),$

corresponding to the eigenvalues $1+\rho_{i}\geq 1-\rho_{i}$ respectively,
with

$\text{Corr}(Y_{i1},Y_{j1})=\frac{\text{Corr}(Z_{i2},Z_{j2})+\text{Corr}(Z_{i2},Z_{j1})+\text{Corr}(Z_{i1},Z_{j2})+\text{Corr}(Z_{i1},Z_{j1})}{2\sqrt{1+\text{Corr}(Z_{i1},Z_{i2})}\sqrt{1+\text{Corr}(Z_{j1},Z_{j2})}}$

and

$\text{Corr}(Y_{i2},Y_{j2})=\frac{\text{Corr}(Z_{i2},Z_{j2})-\text{Corr}(Z_{i2},Z_{j1})-\text{Corr}(Z_{i1},Z_{j2})+\text{Corr}(Z_{i1},Z_{j1})}{2\sqrt{1-\text{Corr}(Z_{i1},Z_{i2})}\sqrt{1-\text{Corr}(Z_{j1},Z_{j2})}}.$

A quick check then verifies that if $((Z_{11},Z_{12}),\dots,(Z_{k1},Z_{k2}))$
has correlation matrix $\mathbf{R}_{m}$ with blocks (9), it holds that
$\text{Corr}(Y_{i1},Y_{j1})=\text{Corr}(Y_{i2},Y_{j2})=1$ for all
$i=1,\dots,k-1$ and $j=i+1,\dots,k$, i.e., all first principal components are
perfectly correlated and all second principal components as well.

Proposition 3 states that the matrix $\mathbf{R}_{m}$ in (7) maximizes the
intensity of dependence, i.e.,
$\mathcal{D}_{1}(\mathbf{R}_{m})=\mathcal{D}_{2}(\mathbf{R}_{m})=1$ when
taking $G_{\mathbf{R}}=\Gamma(\mathbf{R}_{11},\dots,\mathbf{R}_{kk})$ for
fixed marginal covariance matrices $\mathbf{R}_{11},\dots,\mathbf{R}_{kk}$.

######  Proposition 3.

Let $q=d_{1}+\cdots+d_{k}$ and $\mathbf{R}_{ii}\in\mathbb{S}_{\geq}^{d_{i}}$
for $i=1,\dots,k$. The matrix $\mathbf{R}_{m}$ in (7) maximizes
$d_{W}(\mathbf{R},\mathbf{I}_{q})$ and $d_{W}(\mathbf{R},\mathbf{R}_{0})$ with
$\mathbf{R}_{0}$ given in (4) among all
$\mathbf{R}\in\Gamma(\mathbf{R}_{11},\dots,\mathbf{R}_{kk})$.

A general interpretation of $\mathbf{R}_{m}$ is given in Remark 2. .

######  Remark 2.

Assuming again that $d_{1}\leq d_{2}\leq\cdots\leq d_{k}$, the matrix
$\mathbf{R}_{m}$ in (7) is the covariance matrix of

$\begin{pmatrix}\vspace{0.1cm}\mathbf{U}_{11}\boldsymbol{\Lambda}_{11}^{1/2}\mathbf{Z}_{1}\\\
\mathbf{U}_{22}\boldsymbol{\Lambda}_{22}^{1/2}\mathbf{Z}_{2}\\\ \vdots\\\
\mathbf{U}_{kk}\boldsymbol{\Lambda}_{kk}^{1/2}\mathbf{Z}_{k}\end{pmatrix}\in\mathbb{R}^{q\times
1},$

where
$\mathbf{Z}_{1}=(Z_{11},\dots,Z_{1d_{1}})^{\text{T}},\mathbf{Z}_{2}=(Z_{21},\dots,Z_{2d_{2}})^{\text{T}},\dots,\mathbf{Z}_{k}=(Z_{k1},\dots,Z_{kd_{k}})^{\text{T}}$
such that for all $i=1,\dots,k$ we have
$\mathbf{Z}_{i}\sim\mathcal{N}_{d_{i}}(\mathbf{0}_{d_{i}},\mathbf{I}_{d_{i}})$,
and in addition $Z_{ij}=Z_{(i+1)j}$ for all $i=1,\dots,k-1$,
$j=1,\dots,d_{i}$, i.e., $\mathbf{Z}_{i}$ and $\mathbf{Z}_{j}$ have the first
$\min\\{d_{i},d_{j}\\}$ components in common for all $i,j\in\\{1,\dots,k\\}$.
The correlation matrix of the principal components of
$\mathbf{U}_{11}\boldsymbol{\Lambda}_{11}^{1/2}\mathbf{Z}_{1},\dots,\mathbf{U}_{kk}\boldsymbol{\Lambda}_{kk}^{1/2}\mathbf{Z}_{k}$
is

$\begin{pmatrix}\mathbf{I}_{d_{1}}&\boldsymbol{\Pi}_{12}&\dots&\boldsymbol{\Pi}_{1k}\vspace{0.1cm}\\\
\boldsymbol{\Pi}_{12}^{\text{T}}&\mathbf{I}_{d_{2}}&\cdots&\boldsymbol{\Pi}_{2k}\\\
\vdots&\vdots&\ddots&\vdots\\\
\boldsymbol{\Pi}_{1k}^{\text{T}}&\boldsymbol{\Pi}_{2k}^{\text{T}}&\dots&\mathbf{I}_{d_{k}}\end{pmatrix},$

with $\boldsymbol{\Pi}_{ij}$ as in Proposition 2. Hence, if
$(\mathbf{X}_{1},\dots,\mathbf{X}_{k})$ has a Gaussian copula with covariance
matrix $\mathbf{R}_{m}$, then for $\ell=1,\dots,\min\\{d_{i},d_{j}\\}$ the
$\ell$-th principal components of $((\Phi^{-1}\circ
F_{i1})(X_{i1}),\dots,(\Phi^{-1}\circ F_{id_{i}})(X_{id_{i}}))$ and
$((\Phi^{-1}\circ F_{j1})(X_{j1}),\dots,$ $(\Phi^{-1}\circ
F_{jd_{j}})(X_{jd_{j}}))$ are perfectly correlated for all
$i,j\in\\{1,\dots,k\\}$. This is the interpretation of maximal dependence for
the Bures-Wasserstein dependence measures.

In the upcoming section, we also assume that
$G_{\mathbf{R}}=\Gamma(\mathbf{R}_{11},\dots,\mathbf{R}_{kk})$.

###  3.3 Statistical inference

In practice, we have an i.i.d. sample
$\mathbf{X}^{(\ell)}=(\mathbf{X}_{1}^{(\ell)},\dots,\mathbf{X}_{k}^{(\ell)})$
for $\ell=1,\dots,n$ from $\mathbf{X}$, where
$\mathbf{X}_{i}^{(\ell)}=(X_{i1}^{(\ell)},\dots,X_{id_{i}}^{(\ell)})$ for
$\ell=1,\dots,n$ is a sample from $\mathbf{X}_{i}$ for $i=1,\dots,k$. A
natural estimator for the Gaussian copula covariance matrix is known as the
matrix of sample normal scores rank correlation coefficients (see, e.g.,
[24]),

$\widehat{\mathbf{R}}_{n}=\begin{pmatrix}\widehat{\mathbf{R}}_{11}&\widehat{\mathbf{R}}_{12}&\cdots&\widehat{\mathbf{R}}_{1k}\\\
\widehat{\mathbf{R}}_{12}^{\text{T}}&\widehat{\mathbf{R}}_{22}&\cdots&\widehat{\mathbf{R}}_{2k}\\\
\vdots&\vdots&\ddots&\vdots\\\
\widehat{\mathbf{R}}_{1k}^{\text{T}}&\widehat{\mathbf{R}}_{2k}^{\text{T}}&\cdots&\widehat{\mathbf{R}}_{kk}\end{pmatrix}\hskip
5.69046pt\text{with}\hskip
5.69046pt(\widehat{\mathbf{R}}_{im})_{jt}=\widehat{\rho}_{ij,mt}=\frac{\frac{1}{n}\sum_{\ell=1}^{n}\widehat{Z}_{ij}^{(\ell)}\widehat{Z}_{mt}^{(\ell)}}{\frac{1}{n}\sum_{\ell=1}^{n}\left\\{\Phi^{-1}\left(\frac{\ell}{n+1}\right)\right\\}^{2}},$
(10)

defined by computing normal scores

$\widehat{Z}_{ij}^{(\ell)}=\Phi^{-1}\left(\frac{n}{n+1}\widehat{F}_{ij}\left(X_{ij}^{(\ell)}\right)\right)$

obtained through the empirical cdf
$\widehat{F}_{ij}(x_{ij})=\frac{1}{n}\sum_{\ell=1}^{n}\mathbbm{1}\\{X_{ij}^{(\ell)}\leq
x_{ij}\\}$ for $i=1,\dots,k$ and $j=1,\dots,d_{i}$. The quantity
$\widehat{\rho}_{ij,mt}$ is calculated as the conventional Pearson correlation
of the bivariate scores
$((\widehat{Z}_{ij}^{(1)},\widehat{Z}_{mt}^{(1)}),\dots,(\widehat{Z}_{ij}^{(n)},\widehat{Z}_{mt}^{(n)}))$
and by observing that

$\frac{1}{n}\sum_{\ell=1}^{n}\widehat{Z}_{ij}^{(\ell)}=\frac{1}{n}\sum_{\ell=1}^{n}\Phi^{-1}\left(\frac{\ell}{n+1}\right)=0\hskip
5.69046pt\text{and}\hskip
5.69046pt\frac{1}{n}\sum_{\ell=1}^{n}\left(\widehat{Z}_{ij}^{(\ell)}\right)^{2}=\frac{1}{n}\sum_{\ell=1}^{n}\left\\{\Phi^{-1}\left(\frac{\ell}{n+1}\right)\right\\}^{2},$

which holds because $\Phi^{-1}(\alpha)=-\Phi^{-1}(1-\alpha)$ for
$\alpha\in[0,1]$ and $n\widehat{F}_{ij}(X_{ij}^{(\ell)})$ is the rank of
$X_{ij}^{(\ell)}$ in the sample $X_{ij}^{(1)},\dots,X_{ij}^{(n)}$.

A next natural step in estimating $\mathcal{D}_{r}(\mathbf{R})$ is to plug in
$\widehat{\mathbf{R}}_{n}$ for the unknown $\mathbf{R}$. Define the map
$\varphi$ by
$\varphi(\mathbf{\Sigma})=\mathbf{D}_{\mathbf{\Sigma}}^{-1/2}\mathbf{\Sigma}\mathbf{D}_{\mathbf{\Sigma}}^{-1/2}$
for $\mathbf{\Sigma}\in\mathbb{S}^{q}$, where
$\mathbf{D}_{\boldsymbol{\Sigma}}$ is the diagonal matrix containing the
diagonal of $\boldsymbol{\Sigma}$, and let
$||\mathbf{\Sigma}||_{\text{F}}=\text{tr}^{1/2}(\mathbf{\Sigma}^{\text{T}}\mathbf{\Sigma})$
be the Frobenius norm. Fréchet differentiability of the mapping

$(\mathbb{S}^{q},||\cdot||_{\text{F}})\rightarrow(\mathbb{R},|\cdot|):\mathbf{\Sigma}\mapsto(\mathcal{D}_{r}\circ\varphi)(\mathbf{\Sigma})$

on $\mathbb{S}^{q}_{>}$ suffices in order for the delta method to transform an
asymptotic normality result for $\widehat{\mathbf{R}}_{n}$ into an asymptotic
normality result for $\mathcal{D}_{r}(\widehat{\mathbf{R}}_{n})$. We first
highlight some notation. We assume $d_{1}\leq d_{2}\leq\cdots\leq d_{k}$, the
matrix $\mathbf{R}_{m}$ is again defined by (7) based on the
eigendecompositions
$\mathbf{R}_{ii}=\mathbf{U}_{ii}\boldsymbol{\Lambda}_{ii}\mathbf{U}_{ii}^{\text{T}}$
as in (6), and $\mathbf{R}_{0}$ is the matrix in (4). Further, let

$\mathbf{P}_{i}=\left(\mathbf{0}_{d_{i}\times d_{1}}\cdots\hskip
2.84544pt\mathbf{0}_{d_{i}\times d_{i-1}}\hskip
2.84544pt\mathbf{I}_{d_{i}}\hskip 2.84544pt\mathbf{0}_{d_{i}\times
d_{i+1}}\cdots\hskip 2.84544pt\mathbf{0}_{d_{i}\times
d_{k}}\right)\in\mathbb{R}^{d_{i}\times q}$ (11)

be the projection matrix onto the $d_{i}$ coordinates, satisfying
$\mathbf{R}_{ii}=\mathbf{P}_{i}\mathbf{R}\mathbf{P}_{i}^{\text{T}}$. Partition
the matrix $\boldsymbol{\Lambda}_{ii}$ as

$\boldsymbol{\Lambda}_{ii}=\begin{pmatrix}\boldsymbol{\Lambda}_{ii,1}&\mathbf{0}_{d_{1}\times(d_{2}-d_{1})}&\cdots&\mathbf{0}_{d_{1}\times(d_{i}-d_{i-1})}\\\
\mathbf{0}_{d_{1}\times(d_{2}-d_{1})}^{\text{T}}&\boldsymbol{\Lambda}_{ii,2}&\cdots&\mathbf{0}_{(d_{2}-d_{1})\times(d_{i}-d_{i-1})}\\\
\vdots&\vdots&\ddots&\vdots\\\
\mathbf{0}_{d_{1}\times(d_{i}-d_{i-1})}^{\text{T}}&\mathbf{0}_{(d_{2}-d_{1})\times(d_{i}-d_{i-1})}^{\text{T}}&\cdots&\boldsymbol{\Lambda}_{ii,i}\end{pmatrix}\in\mathbb{R}^{d_{i}\times
d_{i}},$ (12)

with
$\boldsymbol{\Lambda}_{ii,j}\in\mathbb{R}^{(d_{j}-d_{j-1})\times(d_{j}-d_{j-1})}$
for $j=1,\dots,i$ and defining $d_{0}=0$. Based on these partitions, further
define

$\begin{split}\boldsymbol{\Delta}_{1}&=\left(\boldsymbol{\Lambda}_{11,1}+\boldsymbol{\Lambda}_{22,1}+\cdots+\boldsymbol{\Lambda}_{kk,1}\right)^{-1/2}\in\mathbb{R}^{d_{1}\times
d_{1}}\\\
\widetilde{\boldsymbol{\Delta}}_{1}&=\left(\boldsymbol{\Lambda}_{11,1}^{2}+\boldsymbol{\Lambda}_{22,1}^{2}+\cdots+\boldsymbol{\Lambda}_{kk,1}^{2}\right)^{-1/2}\boldsymbol{\Lambda}_{11,1}\in\mathbb{R}^{d_{1}\times
d_{1}}\\\
\boldsymbol{\Delta}_{i}&=\begin{pmatrix}\boldsymbol{\Delta}_{i-1}&\mathbf{0}_{d_{i-1}\times(d_{i}-d_{i-1})}\\\
\mathbf{0}_{d_{i-1}\times(d_{i}-d_{i-1})}^{\text{T}}&\left(\boldsymbol{\Lambda}_{ii,i}+\boldsymbol{\Lambda}_{(i+1)(i+1),i}+\cdots+\boldsymbol{\Lambda}_{kk,i}\right)^{-1/2}\end{pmatrix}\in\mathbb{R}^{d_{i}\times
d_{i}}\\\
\widetilde{\boldsymbol{\Delta}}_{i}&=\begin{pmatrix}\widetilde{\boldsymbol{\Delta}}_{i-1}\mathbf{D}_{i}&\mathbf{0}_{d_{i-1}\times(d_{i}-d_{i-1})}\\\
\mathbf{0}_{d_{i-1}\times(d_{i}-d_{i-1})}^{\text{T}}&\left(\boldsymbol{\Lambda}_{ii,i}^{2}+\boldsymbol{\Lambda}_{(i+1)(i+1),i}^{2}+\cdots+\boldsymbol{\Lambda}_{kk,i}^{2}\right)^{-1/2}\boldsymbol{\Lambda}_{ii,i}\end{pmatrix}\in\mathbb{R}^{d_{i}\times
d_{i}}\end{split}$ (13)

for $i=2,\dots,k$ with
$\mathbf{D}_{i}=\text{diag}\left(\boldsymbol{\Lambda}_{(i-1)(i-1),1}^{-1}\boldsymbol{\Lambda}_{ii,1},\dots,\boldsymbol{\Lambda}_{(i-1)(i-1),i-1}^{-1}\boldsymbol{\Lambda}_{ii,i-1}\right)\in\mathbb{R}^{d_{i-1}\times
d_{i-1}}$. Finally, we will need the matrices

$\begin{split}\mathbf{J}=\mathbf{R}_{0}^{-1/2}\left(\mathbf{R}_{0}^{1/2}\mathbf{R}\mathbf{R}_{0}^{1/2}\right)^{1/2}\mathbf{R}_{0}^{-1/2}=\begin{pmatrix}\mathbf{J}_{11}&\mathbf{J}_{12}&\cdots&\mathbf{J}_{1k}\\\
\mathbf{J}_{21}&\mathbf{J}_{22}&\cdots&\mathbf{J}_{2k}\\\
\vdots&\vdots&\ddots&\vdots\\\
\mathbf{J}_{k1}&\mathbf{J}_{k2}&\cdots&\mathbf{J}_{kk}\end{pmatrix}\in\mathbb{R}^{q\times
q}\end{split}$ (14)

with $\mathbf{J}_{ij}\in\mathbb{R}^{d_{i}\times d_{j}}$ and

$\mathbf{J}_{0}=\text{diag}(\mathbf{J}_{11},\dots,\mathbf{J}_{kk})\in\mathbb{R}^{q\times
q}.$ (15)

Theorem 1 formally states an asymptotic normality result for the estimator
$\mathcal{D}_{r}(\widehat{\mathbf{R}}_{n})$.

######  Theorem 1.

Let $\mathbf{X}$ have a Gaussian copula with correlation matrix
$\mathbf{R}\in\mathbb{S}^{q}_{>}$ with $\mathbf{R}\neq\mathbf{R}_{0}$ such
that $\mathbf{R}_{ii}\in\mathbb{S}_{>}^{d_{i}}$ has $d_{i}$ distinct
eigenvalues, and let $\widehat{\mathbf{R}}_{n}$ be given by (10) based on
which the plug-in estimator
$\widehat{\mathcal{D}}_{r,n}=\mathcal{D}_{r}(\widehat{\mathbf{R}}_{n})$ is
constructed. Then for $r\in\\{1,2\\}$, it holds that

$\sqrt{n}\left(\widehat{\mathcal{D}}_{r,n}-\mathcal{D}_{r}(\mathbf{R})\right)\xrightarrow{d}\mathcal{N}\left(\mathbf{0}_{q},\zeta_{r}^{2}\right)$

as $n\to\infty$, with asymptotic variance

$\zeta_{r}^{2}=2\text{tr}\left[\left\\{\mathbf{R}\left(\mathbf{M}_{r}-\mathbf{D}_{\mathbf{M}_{r}\mathbf{R}}\right)\right\\}^{2}\right],$

where $\mathbf{D}_{\mathbf{M}_{r}\mathbf{R}}$ is the diagonal matrix
consisting of the diagonal of $\mathbf{M}_{r}\mathbf{R}$, and

$\mathbf{M}_{1}=\frac{1}{2C_{1}}\left(-\mathbf{R}^{-1/2}+(1-\mathcal{D}_{1}(\mathbf{R}))\mathbf{R}_{0}^{-1/2}+\mathcal{D}_{1}(\mathbf{R})\boldsymbol{\Upsilon}_{1}\right),\hskip
11.38092pt\mathbf{M}_{2}=\frac{1}{C_{2}}\left(-\frac{1}{2}\left(\mathbf{J}_{0}+\mathbf{J}^{-1}\right)+\left(1-\mathcal{D}_{2}(\mathbf{R})\right)\mathbf{I}_{q}+\mathcal{D}_{2}(\mathbf{R})\boldsymbol{\Upsilon}_{2}\right),$

with

$\begin{split}C_{1}=\sum_{i=1}^{k}\text{tr}\left(\mathbf{R}_{ii}^{1/2}\right)-\text{tr}\left(\mathbf{R}_{m}^{1/2}\right),\hskip
28.45274pt\boldsymbol{\Upsilon}_{1}=\begin{pmatrix}\mathbf{U}_{11}\boldsymbol{\Delta}_{1}\mathbf{U}_{11}^{\text{T}}&\mathbf{0}_{d_{1}\times
d_{2}}&\cdots&\mathbf{0}_{d_{1}\times d_{k}}\\\ \mathbf{0}_{d_{1}\times
d_{2}}^{\text{T}}&\mathbf{U}_{22}\boldsymbol{\Delta}_{2}\mathbf{U}_{22}^{\text{T}}&\cdots&\mathbf{0}_{d_{2}\times
d_{k}}\\\ \vdots&\vdots&\ddots&\vdots\\\ \mathbf{0}_{d_{1}\times
d_{k}}^{\text{T}}&\mathbf{0}_{d_{2}\times
d_{k}}^{\text{T}}&\cdots&\mathbf{U}_{kk}\boldsymbol{\Delta}_{k}\mathbf{U}_{kk}^{\text{T}}\end{pmatrix}\end{split}$
(16)

and

$\begin{split}C_{2}=\text{tr}(\mathbf{R})-\text{tr}\left\\{\left(\mathbf{R}_{0}^{1/2}\mathbf{R}_{m}\mathbf{R}_{0}^{1/2}\right)^{1/2}\right\\},\hskip
28.45274pt\boldsymbol{\Upsilon}_{2}=\begin{pmatrix}\mathbf{U}_{11}\widetilde{\boldsymbol{\Delta}}_{1}\mathbf{U}_{11}^{\text{T}}&\mathbf{0}_{d_{1}\times
d_{2}}&\cdots&\mathbf{0}_{d_{1}\times d_{k}}\\\ \mathbf{0}_{d_{1}\times
d_{2}}^{\text{T}}&\mathbf{U}_{22}\widetilde{\boldsymbol{\Delta}}_{2}\mathbf{U}_{22}^{\text{T}}&\cdots&\mathbf{0}_{d_{2}\times
d_{k}}\\\ \vdots&\vdots&\ddots&\vdots\\\ \mathbf{0}_{d_{1}\times
d_{k}}^{\text{T}}&\mathbf{0}_{d_{2}\times
d_{k}}^{\text{T}}&\cdots&\mathbf{U}_{kk}\widetilde{\boldsymbol{\Delta}}_{k}\mathbf{U}_{kk}^{\text{T}}\end{pmatrix}.\end{split}$
(17)

######  Remark 3.

When $\mathbf{R}=\mathbf{R}_{0}$, it holds that $\zeta_{r}=0$ and
$\sqrt{n}\widehat{\mathcal{D}}_{r,n}\xrightarrow{p}0$ for $n\to\infty$. The
higher-order delta method can however still provide a weak convergence result
for $n\widehat{\mathcal{D}}_{r,n}$. A detailed study of this is research in
progress.

We look at another example, which is also studied in [12], but in the context
of $\Phi$-dependence measures.

######  Example 3.

Consider a four dimensional random vector $(X_{1},X_{2},X_{3},X_{4})$ having a
Gaussian copula with covariance matrix

$\begin{pmatrix}1&\rho_{1}&\rho_{2}&\rho_{2}\\\
\rho_{1}&1&\rho_{2}&\rho_{2}\\\ \rho_{2}&\rho_{2}&1&\rho_{1}\\\
\rho_{2}&\rho_{2}&\rho_{1}&1\end{pmatrix},\hskip 8.5359pt\text{where}\hskip
8.5359pt\rho_{1}\geq 2|\rho_{2}|-1.$ (18)

Then one can check that (recall also Example 2 for finding the matrix
$\mathbf{R}_{m}$)

Fig. 1: Dependence coefficients (19) (left) and (20) (right) as a function of
$\rho_{2}$ for different values of $\rho_{1}$.

$\begin{split}\mathcal{D}_{1}\left((X_{1},X_{2});(X_{3},X_{4})\right)=\frac{2(1+\rho_{1})^{1/2}-(\rho_{1}-2\rho_{2}+1)^{1/2}-(\rho_{1}+2\rho_{2}+1)^{1/2}}{\left(2-\sqrt{2}\right)\left\\{(1+\rho_{1})^{1/2}+(1-\rho_{1})^{1/2}\right\\}}\end{split}$
(19)

and

$\mathcal{D}_{2}\left((X_{1},X_{2});(X_{3},X_{4})\right)=\frac{4-2\left|1-\rho_{1}\right|-\left(\rho_{1}^{2}+2\rho_{1}-2\rho_{1}\rho_{2}-2\rho_{2}+1\right)^{1/2}-\left(\rho_{1}^{2}+2\rho_{1}+2\rho_{1}\rho_{2}+2\rho_{2}+1\right)^{1/2}}{4-\sqrt{2}\left\\{\left|1-\rho_{1}\right|+\rho_{1}+1\right\\}}.$
(20)

Fig. 1 shows how (19) and (20) depend on $\rho_{2}$ for different values of
$\rho_{1}$. Clearly, independence holds if and only if $\rho_{2}=0$, as it
should. When $\rho_{1}=2|\rho_{2}|-1$, the second principal component of
$(Z_{1},Z_{2})=((\Phi^{-1}\circ F_{1})(X_{1}),(\Phi^{-1}\circ F_{2})(X_{2}))$
is perfectly correlated with the second principal component of
$(Z_{3},Z_{4})=((\Phi^{-1}\circ F_{3})(X_{3}),(\Phi^{-1}\circ F_{4})(X_{4}))$,
where $F_{i}$ is the marginal distribution of $X_{i}$ for $i=1,2,3,4$, see
also [12]. This causes the $\Phi$-dependence measures studied in [12] to reach
their maximum value, because a singularity is attained. In particular, taking
for instance $\rho_{1}=-0.4$ and $\rho_{2}=0.3$, all (normalized)
$\Phi$-dependence measures equal $1$, while $\mathcal{D}_{1}=0.396$ and
$\mathcal{D}_{2}=0.3$.

Fig. 2: Asymptotic standard deviation of dependence coefficients (19) (left)
and (20) (right) as a function of $\rho_{2}$ for different values of
$\rho_{1}$.

The reason for $\mathcal{D}_{1}$ and $\mathcal{D}_{2}$ still being small, is
that $\rho_{2}$ is pretty small and, recalling Remark 2, not both first and
second principal components of $(Z_{1},Z_{2})$ and $(Z_{3},Z_{4})$ are
perfectly correlated, only one of them is. Only when $|\rho_{2}|=1$ and thus
also $\rho_{1}=1$, we have $\mathcal{D}_{1}=\mathcal{D}_{2}=1$. Maximal Bures-
Wasserstein dependence is not attainable for the family (18) if $\rho_{1}\neq
1$, because it imposes additional restrictions on the correlations (i.e., not
every element in $\Gamma(\mathbf{R}_{11},\dots,\mathbf{R}_{kk})$ is a member
of (18)). Picking a set
$G_{\mathbf{R}}\subset\Gamma(\mathbf{R}_{11},\dots,\mathbf{R}_{kk})$ for a
normalization that adjusts according to these restrictions, would lead to more
cases of maximal dependence.

Fig. 2 depicts the asymptotic standard deviation $\zeta_{r}$ of
$\mathcal{D}_{r}$ for $r\in\\{1,2\\}$ as given in Theorem 1 for this specific
example. We mainly observe increasing behaviour when the strength of
dependence increases. However, in some cases where $\rho_{1}$ and $\rho_{2}$
get close to satisfying $|\rho_{2}|=(\rho_{1}+1)/2$, we see the asymptotic
standard deviation going down. For example, if $\rho_{1}=0$, the asymptotic
standard deviation $\zeta_{1}$ is maximal ($\approx 0.275$) at
$|\rho_{2}|\approx 0.426$, after which it converges to $\approx 0.214$ for
$|\rho_{2}|\to 0.5$. When keeping $\rho_{1}$ fixed and letting
$|\rho_{2}|\to(\rho_{1}+1)/2$, the dependence $\mathcal{D}_{r}$ attains a
local maximum, which it cannot transcend, resulting in (slightly) lower
asymptotic variance. This behaviour was also noticed in the same example in
[12] for their (normalized) $\Phi$-dependence measures, whose asymptotic
variance tends to zero when $|\rho_{2}|\to(\rho_{1}+1)/2$, because the
singularity forces them to reach their global maximum. For the optimal
transport dependence coefficients, the global maximum is only attained when
$|\rho_{2}|,\rho_{1}\to 1$, in which case the asymptotic variance also tends
to zero.

##  4 A Gaussian copula approach: regularized estimation

Empirical covariance matrices are known to approach singularity when the
dimension is close to the sample size. The estimator
$\widehat{\mathcal{D}}_{r,n}$ does not require an inverted covariance matrix,
but it inquires about eigenvalue dispersion, and this tends to be biased when
using the empirical covariance matrix, see, e.g., [31]. Typically, estimates
of large eigenvalues tend to be biased upwards, and estimates of small
eigenvalues tend to be biased downwards. Increasing dimensionality aggravates
this, but penalization techniques can be used to restrain. We briefly discuss
ridge regularization, as in [50], but now in a Gaussian copula context.

Of course, the ridge estimator will not completely shrink elements of the
empirical covariance matrix to zero. However, the task might be to find a
likely estimate for which multiple variables are marginally independent, i.e.,
a sparse estimate of the covariance matrix having zero entries. For this
purpose, we look at a Gaussian copula formulation of the penalization ideas
discussed in [29]. Finally, penalties can also be applied to groups of
elements instead of just individual elements. A group lasso penalty, for
instance, enables to shrink entire blocks of the covariance matrix to zero.

In Section 4.1, we discuss penalization methods for the Gaussian copula
covariance matrix in case $q$ remains fixed with the sample size. Afterwards,
in Section 4.2, we briefly touch upon the case where $q$ depends on $n$.

###  4.1 Fixed dimension

Denote $\widehat{\boldsymbol{\Sigma}}_{\text{ML}}$ for the maximum likelihood
estimator of a $\mathcal{N}_{q}(\mathbf{0}_{q},\boldsymbol{\Sigma})$ model,
and
$\mathbf{R}=\varphi(\boldsymbol{\Sigma})=\mathbf{D}_{\boldsymbol{\Sigma}}^{-1/2}\boldsymbol{\Sigma}\mathbf{D}_{\boldsymbol{\Sigma}}^{-1/2}$
the correlation matrix corresponding to $\boldsymbol{\Sigma}$ (recall that
$\mathbf{D}_{\boldsymbol{\Sigma}}$ is the diagonal matrix containing the
diagonal of $\boldsymbol{\Sigma}$). By adding a certain penalty, say
$P_{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\scriptscriptstyle\bullet$}}}}}}(\boldsymbol{\Sigma},\omega_{n})$
to the Gaussian log-likelihood, where $\omega_{n}$ is a certain penalty
parameter depending on $n$, a general penalized optimization problem for
estimating $\boldsymbol{\Sigma}$ (under the constraint of positive
definiteness $\boldsymbol{\Sigma}\succ\mathbf{0}$), is given by

$\widehat{\boldsymbol{\Sigma}}_{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\scriptscriptstyle\bullet$}}}}}}\in\text{arg
min}_{\boldsymbol{\Sigma}\succ\mathbf{0}}\left\\{\ln\left|\boldsymbol{\Sigma}\right|+\text{tr}\left(\boldsymbol{\Sigma}^{-1}\widehat{\boldsymbol{\Sigma}}_{\text{ML}}\right)+P_{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\scriptscriptstyle\bullet$}}}}}}(\boldsymbol{\Sigma},\omega_{n})\right\\},$
(21)

and corresponding correlation matrix
$\widehat{\mathbf{R}}_{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\scriptscriptstyle\bullet$}}}}}}=\varphi(\widehat{\boldsymbol{\Sigma}}_{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\scriptscriptstyle\bullet$}}}}}})$.
However, the core of this paper is that we merely assume a Gaussian copula.
Hence, instead of making use of $\widehat{\boldsymbol{\Sigma}}_{\text{ML}}$,
we compute $\widehat{\boldsymbol{\Sigma}}_{n}$, being the (block) matrix of
sample normal scores rank covariances with entries
$n^{-1}\sum_{\ell=1}^{n}\widehat{Z}_{ij}^{(\ell)}\widehat{Z}_{mt}^{(\ell)}$
for $i,m\in\\{1,\dots,k\\},j\in\\{1,\dots,d_{i}\\},t\in\\{1,\dots,d_{m}\\}$
(similar block notation as in (10)). The main difference is that we do not
have true Gaussian scores, but only non-parametrically estimated Gaussian
scores
$\widehat{Z}_{ij}^{(\ell)}=\Phi^{-1}(n/(n+1)\widehat{F}_{ij}(X_{ij}^{(\ell)}))$.
The copula formulation of (21) becomes

$\widehat{\boldsymbol{\Sigma}}_{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\scriptscriptstyle\bullet$}}}}},n}\in\text{arg
min}_{\boldsymbol{\Sigma}\succ\mathbf{0}}\left\\{\ln\left|\boldsymbol{\Sigma}\right|+\text{tr}\left(\boldsymbol{\Sigma}^{-1}\widehat{\boldsymbol{\Sigma}}_{n}\right)+P_{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\scriptscriptstyle\bullet$}}}}}}(\boldsymbol{\Sigma},\omega_{n})\right\\},$
(22)

where we use the additional subscript $n$ in
$\widehat{\boldsymbol{\Sigma}}_{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\scriptscriptstyle\bullet$}}}}},n}$
to indicate that we are in the copula context. A genuine Gaussian copula
correlation matrix is then
$\widehat{\mathbf{R}}_{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\scriptscriptstyle\bullet$}}}}},n}=\varphi(\widehat{\boldsymbol{\Sigma}}_{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\scriptscriptstyle\bullet$}}}}},n})$.

We go deeper into three choices for the penalty
$P_{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\scriptscriptstyle\bullet$}}}}}}$:

  * –

the ridge penalty
$P_{\text{R}}(\boldsymbol{\Sigma},\omega_{n})=(1-\omega_{n})\text{tr}(\mathbf{R}^{-1})$,
where $\omega_{n}\in(0,1]$

  * –

(adaptive) lasso-type penalties
$P_{\text{LT}}(\boldsymbol{\Sigma},\omega_{n})=\sum_{i,j,m,t}p_{\omega_{n}}(\Delta_{ij,mt}|\sigma_{ij,mt}|)$,
where $\sigma_{ij,mt}$ is the $(j,t)$’th element of the $(i,m)$’th block of
$\boldsymbol{\Sigma}$ (similar block notation as in (10)), $\Delta_{ij,mt}\geq
0$ is a weight (e.g., zero when $(i,j)=(m,t)$ in order to not shrink diagonal
elements, or larger for smaller preliminary estimated entries in order to
shrink these more) and $p_{\omega_{n}}$ is a certain penalty function
depending on $\omega_{n}\geq 0$

  * –

group lasso-type penalties
$P_{\text{GLT}}(\boldsymbol{\Sigma},\omega_{n})=2\sum_{i,m=1,m>i}^{k}p_{\omega_{n}}(\sqrt{d_{i}d_{m}}\left|\left|\boldsymbol{\Sigma}_{im}\right|\right|_{\text{F}})+\sum_{i=1}^{k}p_{\omega_{n}}(\sqrt{d_{i}(d_{i}-1)}\left|\left|\boldsymbol{\Delta}_{i}*\boldsymbol{\Sigma}_{ii}\right|\right|_{\text{F}})$,
where $||\cdot||_{\text{F}}$ is the Frobenius matrix norm,
$\boldsymbol{\Delta}_{i}\in\mathbb{R}^{d_{i}\times d_{i}}$ a matrix with ones
as off-diagonal elements and zeroes on the diagonal (in order to avoid
shrinking the variances), and $p_{\omega_{n}}$ a certain penalty function
depending on $\omega_{n}\geq 0$.

The ridge penalty is different from the other ones (and will also be
considered separately in the simulations in Section 5) in the sense that it
only focuses on improving the estimation of the Gaussian copula covariance
matrix (and corresponding dependence coefficients) when $q$ is large compared
to $n$, and not on a sparsity assumption. Asymptotic properties are centred
around consistency. Having $q>n$ is definitely manageable for ridge
penalization. For the latter two penalties, the set $\mathcal{A}$ defined as
$\mathcal{A}=\\{\alpha:\boldsymbol{\sigma}_{\alpha}\neq
0,\alpha=1,\dots,\widetilde{q}\\}$, where
$\boldsymbol{\sigma}=\text{vech}(\boldsymbol{\Sigma})\in\mathbb{R}^{\widetilde{q}}$
is the vector of upper diagonal elements of $\boldsymbol{\Sigma}$ and
$\widetilde{q}=q(q-1)/2$, is of crucial importance. Indeed, next to
consistency, we hope that $\mathbb{P}(\mathcal{A}_{n}=\mathcal{A})\to 1$ for
$n\to\infty$, where
$\mathcal{A}_{n}=\\{\alpha:\widehat{\boldsymbol{\sigma}}_{n,\alpha}\neq
0,\alpha=1,\dots,\widetilde{q}\\}$, with
$\widehat{\boldsymbol{\sigma}}_{n}=\text{vech}(\widehat{\boldsymbol{\Sigma}}_{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\scriptscriptstyle\bullet$}}}}},n})$,
a property called sparsistency. Having $q>n$ leads to degeneracy of the lasso-
type estimators, which is why we restrict ourselves to $q\leq n$ in the
simulations.

Example 4 shows that the sparsity of $\boldsymbol{\Sigma}$ (which is obviously
preserved by $\mathbf{R}$), i.e., the entries belonging to $\mathcal{A}^{c}$,
can manifest itself in different forms, calling for different shrinkage
penalties.

######  Example 4.

A covariance graph is a graphical model that represents variables as nodes and
marginal dependencies as edges (similar to a Markov random field representing
conditional dependencies), see, e.g., Section 1 of [4] for several references.
Marginal independencies correspond to individual zeroes in the correlation
matrix and many different sparsity patterns can occur. The first plot of Fig.
3 shows the sparsity structure of a correlation matrix of a random covariance
graph for a total of $20$ variables, where $32.5\%$ of the elements equal
zero, corresponding to the TRUE entries. A penalty of the form $P_{\text{LT}}$
allows reflecting these marginal independencies in the estimated correlation
matrix, and could result in a better plug-in estimator for our dependence
measures.

Fig. 3: Different sparsity patterns of a $20$ dimensional correlation matrix.

Next, imagine a person answering a total of twenty questions in the form of
seven short questionnaires $\mathbf{X}_{1},\dots,\mathbf{X}_{7}$, all
consisting of three questions, except for $\mathbf{X}_{7}$, which only
contains two. The interest is in the relationship between the answers of the
different questionnaires. Furthermore, assume that
$\mathbf{X}_{1},\dots,\mathbf{X}_{6}$ are on completely different topics, and
all three questions are each time self-contained. Only the very last two
questions contained in $\mathbf{X}_{7}$ are related to each other, and to the
other questions in the first six questionnaires as well. Then, it can be
expected that only the blocks $\mathbf{R}_{i7}$ for $i=1,\dots,6$ and
$\mathbf{R}_{77}$ are different from zero (and the diagonal of $\mathbf{R}$,
of course), i.e., $\mathbf{R}$ has a sparsity pattern as visualized in the
second plot of Fig. 3, where $76.5\%$ of the elements equal zero. Such a
pattern calls for a $P_{\text{GLT}}$ penalty, enabling shrinkage of entire
blocks, and possibly more accurate estimation of the Gaussian copula
dependence coefficients compared to using the matrix of normal scores rank
correlations $\widehat{\mathbf{R}}_{n}$, or a $P_{\text{LT}}$ penalty.

We shall now elaborate more on the theoretical properties and computational
aspects of the estimators
$\widehat{\boldsymbol{\Sigma}}_{\text{R},n},\widehat{\boldsymbol{\Sigma}}_{\text{LT},n}$
and $\widehat{\boldsymbol{\Sigma}}_{\text{GLT},n}$ given in (22) corresponding
to the three types of penalties $P_{\text{R}},P_{\text{LT}}$ and
$P_{\text{GLT}}$.

Ridge regularization

Warton [50] extensively studies the estimator
$\widehat{\boldsymbol{\Sigma}}_{\text{R}}$, i.e., the estimator (21) under the
assumption of a normal distribution, with penalty term
$P_{\text{R}}(\boldsymbol{\Sigma},\omega_{n})=(1-\omega_{n})\text{tr}(\mathbf{R}^{-1})$
for
$\mathbf{R}=\varphi(\boldsymbol{\Sigma})=\mathbf{D}_{\mathbf{\Sigma}}^{-1/2}\mathbf{\Sigma}\mathbf{D}_{\mathbf{\Sigma}}^{-1/2}$
and $\omega_{n}\in(0,1]$. His Theorem 1 tells us that

$\widehat{\mathbf{R}}_{\text{R}}=\omega_{n}\widehat{\mathbf{R}}_{\text{ML}}+(1-\omega_{n})\mathbf{I}_{q},$

where
$\widehat{\mathbf{R}}_{\text{ML}}=\varphi(\widehat{\boldsymbol{\Sigma}}_{\text{ML}})$
is the corresponding maximum penalized likelihood estimator. The value of
$\omega_{n}$ is picked through $K$-fold cross-validation with the normal
likelihood as objective function. If $\widehat{\lambda}_{\text{ML}}$ is an
eigenvalue of $\widehat{\mathbf{R}}_{\text{ML}}$, then
$\omega_{n}\widehat{\lambda}_{\text{ML}}+1-\omega_{n}$ is an eigenvalue of
$\widehat{\mathbf{R}}_{\text{R}}$, i.e., all eigenvalues smaller than one are
expanded, and all eigenvalues larger than one are shrunk towards one (recall
the discrepancy of biased eigenvalue dispersion of the non-regularized
empirical covariance matrix), at a pace that increases as $\omega_{n}$
decreases. Furthermore, the value of $\omega_{n}$ tends to one in probability
if $n\to\infty$, and there is zero probability that $\omega_{n}=1$ (no
regularization) if $\widehat{\mathbf{R}}_{\text{ML}}$ does not have full rank,
i.e., $\widehat{\mathbf{R}}_{\text{R}}$ is guaranteed to be non-singular, even
if $q>n$ (see Theorem 2 and Theorem 3 in [50]). In our Gaussian copula
context, we suggest the estimator

$\widehat{\mathbf{R}}_{\text{R},n}=\omega_{n}\widehat{\mathbf{R}}_{n}+(1-\omega_{n})\mathbf{I}_{q},$
(23)

and the same $K$-fold cross-validation procedure for selecting $\omega_{n}$,
but now based on an estimated Gaussian sample. Heuristically (we do not go
into detail here), the asymptotic properties of
$\widehat{\mathbf{R}}_{\text{R}}$ carry over to
$\widehat{\mathbf{R}}_{\text{R},n}$, since we know that (see, e.g., the proof
of Theorem 1)
$||\widehat{\mathbf{R}}_{\text{ML}}-\widehat{\mathbf{R}}_{n}||_{\infty}=\mathchoice{{\scriptstyle\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{{\scriptscriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptscriptstyle\mathcal{O}$}}_{p}(n^{-1/2})$
for $n\to\infty$.

Regularized estimation and sparsity

Concentrating on the fully Gaussian setting, the general penalization
criterion
$P_{\text{LT}}(\boldsymbol{\Sigma},\omega_{n})=\sum_{i,j,m,t}p_{\omega_{n}}(\Delta_{ij,mt}|\sigma_{ij,mt}|)$
is considered in [29]. Typically, one takes $\Delta_{ij,mt}=0$ if
$(i,j)=(m,t)$ and $1$ otherwise in order to avoid penalization of the diagonal
elements, which do not vanish. Another possibility is
$\Delta_{ij,mt}=|\widehat{\sigma}_{ij,mt}|^{-1}\mathds{1}(|\widehat{\rho}_{ij,mt}|<\rho_{n})$,
where $\widehat{\sigma}_{ij,mt}$ is a preliminary estimate for
$\sigma_{ij,mt}$ with corresponding correlation $\widehat{\rho}_{ij,mt}$ and
$0<\rho_{n}<1$ a certain threshold value, i.e., we only shrink those elements
that have a sufficiently small preliminary estimated correlation, and the
amount of shrinkage is proportional to the size of the preliminary estimated
covariance. This is an idea similar to the adaptive lasso of [13] and used in,
e.g., [16]. A similar problem often arises when estimation of the precision
matrix (inverse covariance matrix) is of interest, think of graphical models
for example (see, e.g., [29] for some specific references). It is usually
solved by first performing a local linear approximation to the penalty
function (see, e.g., [52] in case of a precision matrix):

$p_{\omega_{n}}\left(\Delta_{ij,mt}|\sigma_{ij,mt}|\right)\approx
p_{\omega_{n}}\left(\Delta_{ij,mt}\left|\widehat{\sigma}_{ij,mt}^{(w)}\right|\right)+\Delta_{ij,mt}\hskip
1.70709ptp_{\omega_{n}}^{\prime}\left(\Delta_{ij,mt}\left|\widehat{\sigma}_{ij,mt}^{(w)}\right|\right)\left(|\sigma_{ij,mt}|-\left|\widehat{\sigma}_{ij,mt}^{(w)}\right|\right),$

with $\widehat{\sigma}_{ij,mt}^{(w)}$ a current estimated entry of
$\boldsymbol{\Sigma}$ in step $w$. Hence,
$\widehat{\boldsymbol{\Sigma}}_{\text{LT}}^{(w+1)}$ should be taken as
(typically, one iteration already suffices for satisfactory results)

$\widehat{\boldsymbol{\Sigma}}_{\text{LT}}^{(w+1)}\in\text{arg
min}_{\boldsymbol{\Sigma}}\left\\{\ln\left|\boldsymbol{\Sigma}\right|+\text{tr}\left(\boldsymbol{\Sigma}^{-1}\widehat{\boldsymbol{\Sigma}}_{\text{ML}}\right)+\sum_{i,m,j,t}\Delta_{ij,mt}\hskip
1.70709ptp_{\omega_{n}}^{\prime}\left(\Delta_{ij,mt}\left|\widehat{\sigma}_{ij,mt}^{(w)}\right|\right)|\sigma_{ij,mt}|\right\\},$
(24)

which is a weighted covariance graphical lasso problem with weights

$\widetilde{\Delta}_{ij,mt}=\Delta_{ij,mt}\hskip
1.70709ptp_{\omega_{n}}^{\prime}\left(\Delta_{ij,mt}\left|\widehat{\sigma}_{ij,mt}^{(w)}\right|\right).$

We can summarize the $\widetilde{\Delta}_{ij,mt}$ values in a block matrix
$\widetilde{\boldsymbol{\Delta}}$ (again similar as in (10)), allowing us to
rewrite (24) as

$\widehat{\boldsymbol{\Sigma}}_{\text{LT}}^{(w+1)}\in\text{arg
min}_{\boldsymbol{\Sigma}}\left\\{\ln\left|\boldsymbol{\Sigma}\right|+\text{tr}\left(\boldsymbol{\Sigma}^{-1}\widehat{\boldsymbol{\Sigma}}_{\text{ML}}\right)+\left|\left|\widetilde{\boldsymbol{\Delta}}*\boldsymbol{\Sigma}\right|\right|_{1}\right\\},$
(25)

where
$||\mathbf{A}||_{1}=||\text{vec}(\mathbf{A})||_{1}=\sum_{ij}|\mathbf{A}_{ij}|$
is the $L_{1}$-norm of the vector of all elements contained in $\mathbf{A}$,
and $*$ stands for elementwise multiplication. As illustrated in [4], the
optimization (25) consists of a convex part
$\text{tr}(\boldsymbol{\Sigma}^{-1}\widehat{\boldsymbol{\Sigma}}_{\text{ML}})+||\widetilde{\boldsymbol{\Delta}}*\boldsymbol{\Sigma}||_{1}$,
and a concave part $\ln|\boldsymbol{\Sigma}|$, making the entire problem non-
convex (big difference with the precision matrix case, where graphical lasso
is convex), i.e., convergence to a global minimum is not guaranteed. Also,
when $q>n$, the solution to (25) will be degenerate because
$\widehat{\boldsymbol{\Sigma}}_{\text{ML}}$ is not full rank. The authors of
[4] suggest to use
$\widehat{\boldsymbol{\Sigma}}_{\text{ML}}+\delta\mathbf{I}_{q}$ for some
$\delta>0$ in such cases, where $\delta$ is chosen such that, e.g., the
resulting matrix has condition number equal to $q$. Still, they encounter
difficulties of estimation when $q>n$. For this reason, we restrict ourselves
to $q\leq n$ in the simulations.

Using a majorize approach, they propose in [4] to solve convex approximations
to the original problem in an iterative way. Next to sparsity, their algorithm
achieves positive definiteness. Another optimization technique for solving
(25) is developed in [49], who uses coordinate descent, resulting in a faster
and more stable algorithm. We use this algorithm as it is implemented in [16].

Under a set of typical assumptions, consistency and sparsistency results for
the estimator $\widehat{\boldsymbol{\Sigma}}_{\text{LT}}$ are given in [29].
One of their main conclusions is the preference for non-convex penalty
functions such as the scad penalty of [14]:

$p_{\omega_{n}}^{\text{scad}}(t)=\begin{cases}\omega_{n}t&\mbox{for
}t\leq\omega_{n}\\\
\frac{1}{2(a-1)}\left(2a\omega_{n}t-t^{2}-\omega_{n}^{2}\right)&\mbox{for
}\omega_{n}<t\leq a\omega_{n}\\\ (a+1)\omega_{n}^{2}/2&\mbox{for
}t>a\omega_{n}\end{cases},$

where $a>2$. We take $a=3.7$ because of the arguments given in [14]. Such
penalties shrink less entries that are large in magnitude, and as such reduce
the bias. Moreover, a strong theoretical upper bound on the tuning parameter
$\omega_{n}$, as is needed for, e.g., the lasso penalty
$p^{\text{lasso}}_{\omega_{n}}(t)=\omega_{n}t$ (being the limit of
$p_{\omega_{n}}^{\text{scad}}(t)$ for $a\to\infty$) in order to guarantee
consistency, is not needed, yielding better sparsity properties. So, also in
case of sparse covariance matrix estimation, the scad has the oracle property
in the sense of [14] (zeroes are asymptotically estimated as zero, and
estimated non-zeroes are asymptotically normal).

So far, we have been assuming a multivariate normal model, while the core of
this paper is that we merely assume a semi-parametric normal copula model,
i.e., within reach is the estimator
$\widehat{\boldsymbol{\Sigma}}_{\text{LT},n}$, and not
$\widehat{\boldsymbol{\Sigma}}_{\text{LT}}$. As mentioned in [15], who use the
term sparse M-estimator, “such estimation problems have benefited from a very
limited attention so far”, and “the large sample analysis amply differs from
the fully parametric viewpoint”. Recall the sets $\mathcal{A}_{n}$ and
$\mathcal{A}$ introduced for the property of sparsistency. Define also

$a_{n}=\max_{1\leq\alpha\leq\widetilde{q}}\left\\{\frac{dp_{\omega_{n}}}{dt}\left(\left|\widehat{\boldsymbol{\sigma}}_{n,\alpha}\right|\right):\alpha\in\mathcal{A}_{n}\right\\}\hskip
8.5359pt\text{and}\hskip
8.5359ptb_{n}=\max_{1\leq\alpha\leq\widetilde{q}}\left\\{\frac{d^{2}p_{\omega_{n}}}{dt^{2}}\left(\left|\widehat{\boldsymbol{\sigma}}_{n,\alpha}\right|\right):\alpha\in\mathcal{A}_{n}\right\\},$
(26)

where $dp_{\omega_{n}}(|\widehat{\boldsymbol{\sigma}}_{n,\alpha}|)/dt$ denotes
the derivative of $p_{\omega_{n}}(t)$ with respect to $t$, evaluated in
$|\widehat{\boldsymbol{\sigma}}_{n,\alpha}|$, and similarly for the second
derivative. The sequence $a_{n}$ is related to the asymptotic bias of the
penalized estimator, and equals $a_{n}=\omega_{n}$ for the lasso penalty.
Theorem 2 states the consistency and sparsistency of the estimator
$\widehat{\boldsymbol{\Sigma}}_{\text{LT},n}$.

######  Theorem 2.

Assume that there exist $M,\widetilde{M}>0$ such that

$\left|\frac{d^{2}p_{\omega_{n}}}{dt^{2}}\left(t_{1}\right)-\frac{d^{2}p_{\omega_{n}}}{dt^{2}}\left(t_{2}\right)\right|\leq
M\left|t_{1}-t_{2}\right|,$ (27)

for any $t_{1},t_{2}\in\mathbb{R}$ such that
$t_{1},t_{2}>\widetilde{M}\omega_{n}$. Suppose that for $a_{n}$ and $b_{n}$
defined in (26), it holds that $a_{n},b_{n}\to 0$ when $n\to\infty$. Then,
there exists a solution $\widehat{\boldsymbol{\Sigma}}_{\text{LT},n}$ to (22)
with penalty
$P_{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\scriptscriptstyle\bullet$}}}}}}=P_{\text{LT}}$,
for each $n$, satisfying

$\left|\left|\widehat{\boldsymbol{\sigma}}_{n}-\boldsymbol{\sigma}\right|\right|_{2}=\mathcal{O}_{p}\left(\ln(\ln(n))n^{-1/2}+a_{n}\right)$

for $n\to\infty$, where
$\widehat{\boldsymbol{\sigma}}_{n}=\text{vech}(\widehat{\boldsymbol{\Sigma}}_{\text{LT},n})$
and $\boldsymbol{\sigma}=\text{vech}(\boldsymbol{\Sigma})$. If, in addition,
$\omega_{n}\to 0$, $a_{n}\omega_{n}^{-1}\to 0$,
$\sqrt{n}\omega_{n}(\ln(\ln(n)))^{-1}\to\infty$ when $n\to\infty$, and

$\lim\inf_{n\to\infty}\lim\inf_{t\to
0}\omega_{n}^{-1}\frac{dp_{\omega_{n}}}{dt}\left(\left|t\right|\right)>0,$
(28)

it also holds that

$\lim_{n\to\infty}\mathbb{P}\left(\mathcal{A}_{n}=\mathcal{A}\right)=1.$

###### Proof.

This follows immediately from Theorem 3.1. and Theorem 3.2. in [15], noting
that the necessary assumptions are satisfied, since they are verified for the
Gaussian copula case and Gaussian likelihood loss function in Section E of the
Appendix in [15]. ∎

The conditions (27) and (28) are satisfied by, e.g., the lasso and scad
penalty. The condition in (27) is a smoothing condition on the penalty
function, and (28) guarantees sparsity in the estimates. However, note that
the condition $a_{n}\omega_{n}^{-1}\to 0$ for $n\to\infty$ cannot be fulfilled
by the lasso since then $a_{n}=\omega_{n}$. Hence, Theorem 2 does not
guarantee sparsistency of the lasso estimator.

Regularized estimation and group sparsity

Recall the second plot in Fig. 3 of Example 4, where entire groups of
variables are independent of each other, and/or have no dependence within,
resulting in a correlation sparsity pattern with entire zero blocks. A penalty
of the form
$P_{\text{GLT}}(\boldsymbol{\Sigma},\omega_{n})=2\sum_{i,m=1,m>i}^{k}p_{\omega_{n}}(\sqrt{d_{i}d_{m}}\left|\left|\boldsymbol{\Sigma}_{im}\right|\right|_{\text{F}})+\sum_{i=1}^{k}p_{\omega_{n}}(\sqrt{d_{i}(d_{i}-1)}\left|\left|\boldsymbol{\Delta}_{i}*\boldsymbol{\Sigma}_{ii}\right|\right|_{\text{F}})$,
where $||\cdot||_{\text{F}}$ is the Frobenius norm and
$\boldsymbol{\Delta}_{i}\in\mathbb{R}^{d_{i}\times d_{i}}$ a matrix with ones
as off-diagonal elements and zeroes on the diagonal in order to avoid
shrinkage of the variances, might lead to an estimator
$\widehat{\boldsymbol{\Sigma}}_{\text{GLT},n}$ that incorporates this sparsity
structure. The non-differentiability of $||\cdot||_{\text{F}}$ at $\mathbf{0}$
allows to simultaneously shrink all entries of a certain block, and hence
performs a group selection. For actually computing
$\widehat{\boldsymbol{\Sigma}}_{\text{GLT},n}$, a local linear approximation
of $p_{\omega_{n}}$ can again be performed, which will boil down to a similar
problem as finding $\widehat{\boldsymbol{\Sigma}}_{\text{GLT},n}$ in case of
the lasso penalty $p^{\text{lasso}}_{\omega_{n}}(t)=\omega_{n}t$. So, let us
assume from now on that $p_{\omega_{n}}=p_{\omega_{n}}^{\text{lasso}}$. The
ideas explained in [49], which we use for numerically finding
$\widehat{\boldsymbol{\Sigma}}_{\text{LT},n}$, can in general not be used for
finding $\widehat{\boldsymbol{\Sigma}}_{\text{GLT},n}$, because some of the
arguments do not apply to the Frobenius norm. Nevertheless, the problem

$\widehat{\boldsymbol{\Sigma}}_{\text{GLT},n}^{*}=\text{arg
min}_{\boldsymbol{\Sigma}}\left\\{\frac{1}{2}\left|\left|\boldsymbol{\Sigma}-\widehat{\boldsymbol{\Sigma}}_{n}\right|\right|_{\text{F}}^{2}+2\omega_{n}\sum_{\begin{subarray}{c}i,m=1\\\
m>i\end{subarray}}^{k}\sqrt{d_{i}d_{m}}\left|\left|\boldsymbol{\Sigma}_{im}\right|\right|_{\text{F}}+\omega_{n}\sum_{i=1}^{k}\sqrt{d_{i}(d_{i}-1)}\left|\left|\boldsymbol{\Delta}_{i}*\boldsymbol{\Sigma}_{ii}\right|\right|_{\text{F}}\right\\}$
(29)

has a solution in the form of an elementwise soft thresholding operation.
Denoting $(\widehat{\boldsymbol{\Sigma}}^{*}_{\text{GLT},n})_{ij,mt}$ for the
$(j,t)$’th element of the $(i,m)$’th block of
$\widehat{\boldsymbol{\Sigma}}^{*}_{\text{GLT},n}$, similarly for
$\widehat{\boldsymbol{\Sigma}}_{n}$, and further

$\mathcal{S}_{i,m}=\sqrt{\mathop{\sum_{j=1}^{d_{i}}\sum_{t=1}^{d_{m}}}_{(i,j)\neq(m,t)}\left\\{\left(\widehat{\boldsymbol{\Sigma}}_{n}\right)_{ij,mt}\right\\}^{2}},\hskip
8.5359pt\text{and}\hskip
8.5359pt\gamma_{i,m}=\begin{cases}\sqrt{d_{i}d_{m}}&\mbox{if }i\neq m\\\
\sqrt{d_{i}(d_{i}-1)}&\mbox{if }i=m,\end{cases}$

it is known that (similar to, e.g., Proposition 1 in [6])

$\left(\widehat{\boldsymbol{\Sigma}}^{*}_{\text{GLT},n}\right)_{ij,mt}=\begin{cases}0&\mbox{if
}\mathcal{S}_{i,m}\leq\omega_{n}\gamma_{i,m}\hskip 2.84544pt\text{and}\hskip
2.84544pt(i,j)\neq(m,t)\\\
\left(\widehat{\boldsymbol{\Sigma}}_{n}\right)_{ij,mt}&\mbox{if
}(i,j)=(m,t)\\\
\left(\widehat{\boldsymbol{\Sigma}}_{n}\right)_{ij,mt}\left(1-\frac{\omega_{n}\gamma_{i,m}}{\mathcal{S}_{i,m}}\right)&\mbox{if
}\mathcal{S}_{i,m}>\omega_{n}\gamma_{i,m}\hskip 2.84544pt\text{and}\hskip
2.84544pt(i,j)\neq(m,t).\end{cases}$

Hence, we can numerically compute
$\widehat{\boldsymbol{\Sigma}}_{\text{GLT},n}$ by using ideas similar to the
optimization approach of [4], i.e., by majorizing $\ln|\boldsymbol{\Sigma}|$
by its tangent plane and using generalized gradient descent steps, which comes
down to iteratively solving (convex) problems of the form (29). For the actual
implementation, we used to code available in [5], and fine-tuned it such that
it can also cope with a group lasso penalty. Regarding the asymptotic
properties, it is intuitively clear that a result similar to Theorem 2 will
hold, with sparsistency formulated at the level of blocks instead of
individual elements. For asymptotic properties in case of truly independent
copies, we refer to Section 3 of [6].

###  4.2 Dimension depending on the sample size

So far, we have been assuming that the total dimension $q$ of the random
vector $\mathbf{X}$ remains fixed with $n$. When
$\mathbf{X}\in\mathbb{R}^{q_{n}}$ for $q_{n}$ depending on $n$, of primary
interest might be the behaviour of the empirical Gaussian copula covariance
matrix $\widehat{\mathbf{R}}_{n}$ given in (10). Proposition 4 (see D for a
proof) tells us how the dimension $q_{n}$ influences the consistency of the
estimator $\widehat{\mathbf{R}}_{n}$ in max norm
$||\mathbf{A}||_{\infty}=\max_{i,j}|\mathbf{A}_{ij}|$.

######  Proposition 4.

Let $q_{n}$ be a sequence of dimensions depending on $n$, and
$\mathbf{R}_{n}\in\mathbb{R}^{q_{n}\times q_{n}}$ corresponding Gaussian
copula covariance matrices. Assume that
$\sup_{n}\lambda_{\max}(\mathbf{R}_{n})<\infty$, where $\lambda_{\max}$
denotes the maximum eigenvalue, and let $\widehat{\mathbf{R}}_{n}$ be the
estimator given in (10) with $q=q_{n}$. Then, it holds that

$\left|\left|\widehat{\mathbf{R}}_{n}-\mathbf{R}_{n}\right|\right|_{\infty}=\mathcal{O}_{\text{p}}\left\\{\left(\ln(q_{n})/n\right)^{1/2}\right\\},$

for $n\to\infty$.

So, Proposition 4 ensures that $\widehat{\mathbf{R}}_{n}$ is a consistent
estimator as long as $\ln(q_{n})/n\to 0$ as $n\to\infty$, i.e., as long as we
are not in an ultra-high dimensional setting. In the context of the
penalization techniques, one can also assume that $q=q_{n}$ depends on $n$.
For the ridge regularization, this will lead to inconsistencies in high-
dimensional settings, see Section 3.1 in [50] (basically because sample
eigenvalues are known to be inconsistent when $q_{n}\to\infty$). Regarding the
$P_{\text{LT}}$ penalties, Theorem H.1. in [15] states that

$\left|\left|\widehat{\boldsymbol{\sigma}}_{n}-\boldsymbol{\sigma}\right|\right|_{2}=\mathcal{O}_{p}\left\\{\sqrt{\hskip
1.42271pt\widetilde{q}_{n}}\left(\ln(\ln(n))n^{-1/2}+a_{n}\right)\right\\},$

where $\widetilde{q}_{n}=q_{n}(q_{n}-1)/2$, and under the additional
conditions (next to those of Theorem 2) that
$\widetilde{q}_{n}^{2}\ln(\ln(n))n^{-1/2}\to 0$ and
$\widetilde{q}_{n}^{2}a_{n}\to 0$ for $n\to\infty$. Hence, $q_{n}$ is allowed
to increase with $n$, but consistency requires
$\widetilde{q}_{n}=\mathchoice{{\scriptstyle\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{{\scriptscriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptscriptstyle\mathcal{O}$}}(n^{1/4})$.
If, in addition (see Theorem H.2. in [15]),
$\widetilde{q}_{n}a_{n}\omega_{n}^{-1}$ and
$\sqrt{n}\omega_{n}(\widetilde{q}_{n}\ln(\ln(n)))^{-1}\to\infty$ for
$n\to\infty$, true zeroes are asymptotically identified with probability one
(but there might still be some false positives).

##  5 Simulation study

The aim of this section is to empirically study the finite sample performance
of the (non-)regularized plug-in estimators for the Gaussian copula based
dependence coefficients $\mathcal{D}_{1}$ and $\mathcal{D}_{2}$ discussed in
Section 3. The finite-sample distribution of the estimator of Section 3, is
investigated in Section 5.1, and the regularized estimators are considered in
Section 5.2.

###  5.1 Finite-sample performance of the estimator of Section 3

Theorem 1 gives an asymptotic normality result with explicit asymptotic
variance for the estimator
$\widehat{\mathcal{D}}_{r,n}=\mathcal{D}_{r}(\widehat{\mathbf{R}}_{n})$ for
$r\in\\{1,2\\}$, with $\widehat{\mathbf{R}}_{n}$ the matrix of sample normal
scores rank correlations given in (10). Let now $\widehat{\zeta}_{r,n}$ be the
plug-in estimator of the asymptotic standard deviation $\zeta_{r}$ obtained by
replacing $\mathbf{R}$ with $\widehat{\mathbf{R}}_{n}$. By simulating a sample
from a certain multivariate distribution having a Gaussian copula, we can
compute a realization of the actual sampling distribution of the studentized
estimator
$\sqrt{n}(\widehat{\mathcal{D}}_{r,n}-\mathcal{D}_{r})/\widehat{\zeta}_{r,n}$,
and several replications can be used to represent characteristics of the
entire distribution, which should approximately (for large $n$) be a standard
normal one according to Theorem 1.

We consider four settings which we can generate samples from:

Fig. 4: Standard normal Q–Q plots for $1000$ Monte Carlo runs of the
studentized plug-in estimator for $\mathcal{D}_{1}$ in four different settings
with sample sizes $n=50,200,1000,5000$. The median (“Med”) of the studentized
estimates is given in blue.

  * 1.

Setting 1: $k=2,d_{1}=d_{2}=2$, with standard normal marginals and a Gaussian
copula having an autoregressive AR(1) correlation matrix with $\rho=0.25$
(i.e., the $(i,j)$’th element of $\mathbf{R}$ equals $0.25^{|i-j|}$).

  * 2.

Setting 2: as Setting 1, but now with marginals

    * (a)

a $t$ distribution with $3$ degrees of freedom for $X_{11}$

    * (b)

an exponential distribution with mean $1$ for $X_{12}$

    * (c)

a beta distribution with parameters $2$ and $2$ for $X_{21}$

    * (d)

an $F$-distribution with degrees of freedom $2$ and $6$ for $X_{22}$.

  * 3.

Setting 3: similar to Setting 1, but with $\rho=0.8$.

  * 4.

Setting 4: $k=5,d_{1}=4,d_{2}=5,d_{3}=3,d_{4}=1,d_{5}=2$, with standard normal
marginals and a Gaussian copula having an equicorrelated correlation matrix
with $\rho=0.5$.

Each time, we draw $1000$ samples of sizes $n=50,200,1000,5000$ and make
standard normal Q–Q plots to assess the goodness-of-fit with a standard normal
distribution. See Fig. 4 for the results in case of $\mathcal{D}_{1}$. Similar
plots are obtained for $\mathcal{D}_{2}$. In each setting, there is a
qualitative fit with the standard normal distribution for larger sample sizes.
Changing the marginals has no influence (Setting 1 versus Setting 2). Small
correlations (Settings 1 and 2), yield a more pronounced lack-of-fit than
higher correlations (Setting 3). Increasing the total dimension to $15$
(Setting 4) results in a large positive bias for small sample sizes, calling
for regularization techniques.

Next, in Fig. 2, asymptotic standard deviations $\zeta_{1}$ and $\zeta_{2}$
were shown in the context of Example 3. Generating $N$ samples from, e.g., a
multivariate normal distribution with mean zero and covariance matrix given in
(18), we obtain $N$ estimates
$\widehat{\mathcal{D}}_{r,n}^{(1)},\dots,\widehat{\mathcal{D}}_{r,n}^{(N)}$,
whose sample standard deviation multiplied with $\sqrt{n}$, say
$\widehat{\zeta}_{r}$, can be seen as an approximation for $\zeta_{r}$ when
$n$ is large. Fig. 5 depicts $\widehat{\zeta}_{r}$ for different values of
$\rho_{1}$ and $\rho_{2}$ in case $n=10\hskip 1.42271pt000$ and $N=1000$. We
see the same patterns popping up as in Fig. 2, illustrating the asymptotic
variance formula given in Theorem 1 empirically in this particular setting.

Fig. 5: Empirical standard deviation (sample size $N=1000$)
$\widehat{\zeta}_{r}$ for $r\in\\{1,2\\}$ in the setting of Example 3 when
$n=10\hskip 1.42271pt000$.

###  5.2 Penalization techniques

We now turn our attention to the different covariance matrix penalization
techniques discussed in Section 4. We start with illustrating ridge
regularization, in particular the improvements it gives in plug-in estimation
of Gaussian copula based dependence coefficients between multiple random
vectors. We do this for increasing values of $q$, possibly larger than $n$.
Recall that the ridge estimator is rather easy to compute, and there are no
additional difficulties when $q>n$. This study thus also allows for an
impression on how the dependence coefficients behave with increasing $q$.
Afterwards, we go to sparsity inducing (group)-lasso type methods, and
investigate for two fixed values of $q$ and different sample sizes (with
$q<n$) their ability of recovering marginal independencies (interpretability)
on the one hand, and whether this improves the estimation of the dependence
coefficients (accuracy) on the other hand.

Ridge regularization

As explained in Section 4, the ridge estimator
$\widehat{\mathbf{R}}_{\text{R},n}$ of the Gaussian copula correlation matrix
tries to cope with biased eigenvalue dispersion of the empirical correlation
matrix which aggravates when $q$ is large compared to $n$. Moreover, it is
easy to implement (with a straightforward cross-validation procedure for
selecting the penalty parameter $\omega_{n}$), and guarantees a positive
definite outcome. One can expect that the performance of the estimator
$\mathcal{D}_{r}(\widehat{\mathbf{R}}_{\text{R},n})$ is better than the
performance of $\mathcal{D}_{r}(\widehat{\mathbf{R}}_{n})$ when $q$ is large
compared to $n$. To demonstrate this, we consider the following two designs:

  * 1.

Design 1: We let $q\in\\{4,6,8,\dots,98\\}$ and $k=q/2$ with
$d_{1}=d_{2}=\cdots=d_{k}=2$, e.g., if $q=30$, we are measuring the dependence
between $15$ random vectors of dimension $2$

  * 2.

Design 2: We let $q\in\\{4,6,8,\dots,98\\}$ and $k=2$ with $d_{1}=d_{2}=q/2$,
e.g., if $q=30$, we are measuring the dependence between $2$ random vectors of
dimension $15$.

In each design, we generate $1000$ samples of sizes $n=50,100,500$ from a
$\mathcal{N}_{q}(\mathbf{0}_{q},\mathbf{R})$ distribution (but, assume unknown
marginals), where $\mathbf{R}$ is the correlation matrix of an $\text{AR}(1)$
process with $\rho=0.5$, and compute the empirical mean squared error of
$\mathcal{D}_{r}(\widehat{\mathbf{R}}_{\text{R},n})$ and
$\mathcal{D}_{r}(\widehat{\mathbf{R}}_{n})$.

Fig. 6: Logarithm of the Monte Carlo mean squared error of the estimators
$\mathcal{D}_{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\scriptscriptstyle\bullet$}}}}}}(\widehat{\mathbf{R}}_{\text{R},n})$
(penalty) and
$\mathcal{D}_{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\scriptscriptstyle\bullet$}}}}}}(\widehat{\mathbf{R}}_{n})$
(no penalty) for
$\mathcal{D}_{t\ln(t)},\mathcal{D}_{(\sqrt{t}-1)^{2}},\mathcal{D}_{1}$ and
$\mathcal{D}_{2}$, based on $1000$ replications with sample sizes
$n=50,100,500$ as a function of $q$, in two different designs.

The penalty parameter $\omega_{n}$ is determined by $5$-fold cross-validation
(as described by equation (4) in [50]) on a grid of $50$ equidistant elements
in $[0.01,0.999]$. In addition, we do the same for two Gaussian copula-based
$\Phi$-dependence measures discussed in [12], known as the (normalized) mutual
information and Hellinger distance, respectively given by

$\mathcal{D}_{t\ln(t)}(\mathbf{R})=\left(1-\frac{\left|\mathbf{R}\right|}{\prod_{i=1}^{k}\left|\mathbf{R}_{ii}\right|}\right)^{1/2},\hskip
5.69046pt\text{and}\hskip
5.69046pt\mathcal{D}_{(\sqrt{t}-1)^{2}}(\mathbf{R})=1-\frac{2^{q/2}\left|\mathbf{R}\right|^{1/4}}{\left|\mathbf{I}_{q}+\mathbf{R}_{0}^{-1}\mathbf{R}\right|^{1/2}\prod_{i=1}^{k}\left|\mathbf{R}_{ii}\right|^{1/4}},$
(30)

where $\mathbf{R}_{0}$ is given in (4), and $|\cdot|$ denotes the determinant.
Note that the dependence coefficients in (30) depend on products of
eigenvalues (they are based on divergences of copula densities), while
$\mathcal{D}_{1}$ and $\mathcal{D}_{2}$ depend on sums of eigenvalues (arising
from the Bures-Wasserstein distance).

Fig. 6 shows plots of the logarithm of the Monte Carlo mean squared error of
$\mathcal{D}_{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\scriptscriptstyle\bullet$}}}}}}(\widehat{\mathbf{R}}_{\text{R},n})$
(penalty) and
$\mathcal{D}_{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\scriptscriptstyle\bullet$}}}}}}(\widehat{\mathbf{R}}_{n})$
(no penalty) for the four dependence measures for each sample size in both
designs as a function of $q$. In all cases, we see vast improvements when
using ridge regularization, especially when $n$ is small compared to $q$ (note
that when $n=50$ and $q>50$, the estimators
$\mathcal{D}_{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\scriptscriptstyle\bullet$}}}}}}(\widehat{\mathbf{R}}_{n})$
are not defined because of singularity of $\widehat{\mathbf{R}}_{n})$. Fig. 7
shows boxplots of the selected values for the penalization parameter
$\omega_{n}$. We clearly see compatibility with the theoretical property that
$\omega_{n}$ tends to one in probability when $n\to\infty$, and higher
dimensions require a larger correction of the eigenvalue dispersion, i.e., a
smaller $\omega_{n}$.

Fig. 7: Boxplots of selected penalty parameter $\omega_{n}$ via a $5$-fold
cross-validation search on a grid of $50$ equidistant elements in
$[0.01,0.999]$, for different values of $n$ and $q$.

Finally, Fig. 6 also indicates that it would be interesting to study the
behaviour of the dependence coefficients when $q\to\infty$, which can happen
in multiple ways. In the first design for example, $d_{i}$ remains fixed for
all $i=1,\dots,k$, but $k\to\infty$. Because $\Phi$-dependence measures
satisfy Axiom (A4) of A (see, e.g., [12]), we know that
$\mathcal{D}_{t\ln(t)},\mathcal{D}_{(\sqrt{t}-1)^{2}}\to 1$ when $k\to\infty$,
which is probably why the mean squared error of
$\widehat{\mathcal{D}}_{t\ln(t)},\widehat{\mathcal{D}}_{(\sqrt{t}-1)^{2}}$
first increases and afterwards decreases/becomes constant (in particular, the
variance tends to $0$ when $k\to\infty$, and the bias decreases/becomes
constant, see Fig. S1 and S2 of the Supplementary Material for plots of the
bias and variance). Regarding the optimal transport dependence measures
$\mathcal{D}_{1}$ and $\mathcal{D}_{2}$, we expect them (based on our
simulations) to converge to $0$ when $k\to\infty$, probably because the strong
normalization (denominator) grows faster than the numerator does. We see the
variance of $\widehat{\mathcal{D}}_{r}$ decreasing in $q$, and the biases
increasing in $q$.

In the second design, $k=2$ remains fixed, but $d_{1}=d_{2}\to\infty$. The
optimal transport measures seem (in our simulations) to converge to zero
again, while the $\Phi$-dependence measures remain constant. We leave a formal
study of the behaviour of these dependence coefficients when $q\to\infty$ for
further research.

(Adaptive/Group) lasso-type estimation

Recall the two different sparsity patterns discussed in Example 4 for a
$20\times 20$ correlation matrix. By performing lasso-type estimation, we hope
to recover zero entries (interpretability) on the one hand, and improve
accuracy on the other hand. In particular, for the latter, we desire better
performance of the plug-in estimator
$\widehat{\mathcal{D}}=\mathcal{D}(\widehat{\mathbf{R}}_{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.9}{$\scriptscriptstyle\bullet$}}}}},n})$
than of the non-penalized estimator $\mathcal{D}(\widehat{\mathbf{R}}_{n})$,
where $\mathcal{D}$ is a certain correlation matrix based dependence
coefficient.

For the simulation, we consider the following correlation structures

  * 1.

Scenario 1: We take $q=20$ with $k=7$, $d_{1}=\cdots=d_{6}=3$ and $d_{7}=2$
having the non-block sparse correlation structure given in the left plot of
Fig. 3, with $32.5\%$ zeroes.

  * 2.

Scenario 2: We take $q=20$ with $k=7$, $d_{1}=\cdots=d_{6}=3$ and $d_{7}=2$
such that $\mathbf{X}_{1},\dots,\mathbf{X}_{6}$ have independence within and
between each other, but all components of $\mathbf{X}_{7}$ are related and all
dependent on all components of $\mathbf{X}_{1},\dots,\mathbf{X}_{6}$, i.e.,
the block sparse correlation structure of Example 4, right plot in Fig. 3,
with $76.5\%$ zeroes.

  * 3.

Scenario 3: We take $q=70$ with $k=4$, $d_{1}=d_{2}=d_{3}=21$ and $d_{4}=7$
such that there is only dependence with and within $\mathbf{X}_{4}$, yielding
a block sparse correlation structure, with $79.7\%$ zeroes.

Each time, we generate $1000$ samples from a
$\mathcal{N}_{q}(\mathbf{0}_{q},\mathbf{R})$ distribution (marginals are again
assumed to be unknown), and compute the no penalty, lasso, adaptive lasso (of
[16] with $\omega_{n}=1$ and tuning parameter $\rho_{n}$), scad, and group
lasso estimator for $\mathbf{R}$. The considered sample sizes are
$n=50,100,500$ for Scenarios 1 and 2, and $n=100,500$ for Scenario 3. Based on
the $1000$ replications, we report on the average true positive rate (TPR),
which we want to be close to one, and false positive rate (FPR), which we want
to be close to zero. We also compute the empirical root mean squared error of
$||\widehat{\mathbf{R}}-\mathbf{R}||_{\text{F}}/q$ (where
$\widehat{\mathbf{R}}$ is the estimated correlation matrix in question), and
empirical mean squared error of $\widehat{\mathcal{D}}$, with $\mathcal{D}$
the mutual information $\mathcal{D}_{t\ln(t)}$, Hellinger distance
$\mathcal{D}_{(\sqrt{t}-1)^{2}}$, or one of the optimal transport dependence
measures $\mathcal{D}_{1}$ or $\mathcal{D}_{2}$.

For tuning $\omega_{n}$ (or $\rho_{n}$ in case of the adaptive lasso), we use
the BIC criterion (see, e.g., [17] for the case of a precision matrix in
Gaussian graphical models):

$\text{BIC}\left(\widehat{\boldsymbol{\Sigma}}_{\omega_{n}}\right)=-n\left[\ln\left|\widehat{\boldsymbol{\Sigma}}_{\omega_{n}}\right|+\text{tr}\left(\widehat{\boldsymbol{\Sigma}}_{\omega_{n}}^{-1}\widehat{\boldsymbol{\Sigma}}_{n}\right)\right]-\ln(n)\text{df}\left(\widehat{\boldsymbol{\Sigma}}_{\omega_{n}}\right),$
(31)

where $\widehat{\boldsymbol{\Sigma}}_{\omega_{n}}$ is the estimated candidate
covariance matrix using the penalty parameter $\omega_{n}$. We want to
maximize (31) in $\omega_{n}$. We do this over an equidistant grid of $50$
values in $[0.01;0.6]$. The number
$\text{df}(\widehat{\boldsymbol{\Sigma}}_{\omega_{n}})$ stands for the degrees
of freedom, and is estimated for the (adaptive) lasso and scad by the number
of non-zero entries in $\widehat{\boldsymbol{\Sigma}}_{\omega_{n}}$, not
taking the elements under the diagonal into account. For the group lasso, the
degrees of freedom can be estimated in a similar spirit as equation (23) in
[9]:

$\begin{split}\text{df}\left(\widehat{\boldsymbol{\Sigma}}_{\omega_{n}}\right)=\sum_{\begin{subarray}{c}i,m=1\\\
m>i\end{subarray}}^{k}\mathds{1}\left(\left|\left|\widehat{\boldsymbol{\Sigma}}_{\omega_{n},im}\right|\right|_{\text{F}}>0\right)&\left(1+\frac{\left|\left|\widehat{\boldsymbol{\Sigma}}_{\omega_{n},im}\right|\right|_{\text{F}}}{\left|\left|\widehat{\boldsymbol{\Sigma}}_{n,im}\right|\right|_{\text{F}}}\left(d_{i}d_{m}-1\right)\right)\\\
&+\sum_{i=1}^{k}\mathds{1}\left(\left|\left|\boldsymbol{\Delta}_{i}*\widehat{\boldsymbol{\Sigma}}_{\omega_{n},ii}\right|\right|_{\text{F}}>0\right)\left(1+\frac{\left|\left|\boldsymbol{\Delta}_{i}*\widehat{\boldsymbol{\Sigma}}_{\omega_{n},ii}\right|\right|_{\text{F}}}{\left|\left|\boldsymbol{\Delta}_{i}*\widehat{\boldsymbol{\Sigma}}_{n,ii}\right|\right|_{\text{F}}}\left(\frac{d_{i}(d_{i}-1)}{2}-1\right)\right)+q,\end{split}$

where $\widehat{\boldsymbol{\Sigma}}_{\omega_{n},im}$ is the $(i,m)$’th block
of $\widehat{\boldsymbol{\Sigma}}_{\omega_{n}}$, similarly for
$\widehat{\boldsymbol{\Sigma}}_{n,im}$, and
$\boldsymbol{\Delta}_{ii}\in\mathbb{R}^{d_{i}\times d_{i}}$ is a matrix with
ones as off-diagonal elements and zeroes on the diagonal. The results are
summarized in Table 1.

In Scenario 1 ($q=20$, no block sparsity, $32.5\%$ zeroes), there is no group
sparsity, which results in low TPR (and also low FPR, since only few entries
get shrunk to zero) for the group lasso estimator. The other penalization
techniques are, as expected, preferred for recovering zeroes. Clearly, any
type of considered penalization yields more accurate estimation of
$\mathbf{R}$ in Frobenius norm than in case no penalty is used. However, this
does not necessarily imply better estimation of the dependence coefficients,
especially for lasso and scad. For good estimation of these, it is important
not to lose sight of the $67.5\%$ non-zero entries, which is why the adaptive
lasso performs really well. Note that the mutual information is estimated very
well in the non-penalized case, but this is mainly due to the fact that the
true value equals $0.994$, which is very close to one, being an effect of the
dimension (recall also Fig. 6). All non-penalized mutual information estimates
are close to one because of a relatively large $q$, yielding low estimation
error.

In Scenario 2 ($q=20$, block sparsity, $76.5\%$ zeroes), all penalization
techniques perform well in identifying zeroes. For obtaining both high TPR and
low FPR, the group lasso performs slightly better than lasso and scad. Also,
when the focus is on estimating $\mathbf{R}$, penalization is clearly
beneficial, and the group lasso outperforms the other techniques. Regarding
the dependence coefficients, we see that, for the optimal transport measures,
lasso and scad give improvement compared to using no penalty, especially for
smaller sample sizes. The error of the $\Phi$-dependence

Table 1: True positive rate (TPR), false positive rate (FPR), and empirical
mean squared error for the estimated correlation matrix and corresponding
plug-in estimators of dependence coefficients based on $1000$ replications in
different scenarios, using different penalties. Scenario 1 ($q=20$, no block
sparsity, $32.5\%$ zeroes)
---
| n | TPR | FPR | $||\widehat{\mathbf{R}}-\mathbf{R}||_{\text{F}}/q$ | $\widehat{\mathcal{D}}_{t\ln(t)}$ | $\widehat{\mathcal{D}}_{(\sqrt{t}-1)^{2}}$ | $\widehat{\mathcal{D}}_{1}$ | $\widehat{\mathcal{D}}_{2}$ |
| $50$ | $\backslash$ | $\backslash$ | $0.137$ | $3.993\cdot 10^{-4}$ | $0.109$ | $0.058$ | $0.057$ |
no penalty | $100$ | $\backslash$ | $\backslash$ | $0.096$ | $2.898\cdot 10^{-4}$ | $0.044$ | $0.026$ | $0.026$ |
| $500$ | $\backslash$ | $\backslash$ | $0.043$ | $8.331\cdot 10^{-5}$ | $0.009$ | $0.005$ |  $0.005$ |
| $50$ | $0.831$ | $0.592$ | $0.105$ | $0.209$ | $0.327$ | $0.054$ | $0.052$ |
lasso | $100$ | $0.724$ | $0.355$ | $0.077$ | $0.005$ | $0.153$ | $0.032$ | $0.031$ |
| $500$ | $0.554$ | $0.121$ | $0.039$ | $0.001$ | $0.054$ | $0.014$ |  $0.014$ |
| $50$ | $0.818$ | $0.565$ | $0.122$ | $0.002$ | $0.075$ | $0.025$ | $0.025$ |
adaptive lasso | $100$ | $0.848$ | $0.480$ | $0.092$ | $2.640\cdot 10^{-4}$ | $0.030$ | $0.012$ | $0.013$ |
| $500$ | $0.932$ | $0.303$ | $0.044$ | $8.041\cdot 10^{-5}$ | $0.010$ | $0.003$ |  $0.004$ |
| $50$ | $0.831$ | $0.592$ | $0.105$ | $0.209$ | $0.327$ | $0.054$ | $0.053$ |
scad | $100$ | $0.744$ | $0.373$ | $0.078$ | $0.006$ | $0.165$ | $0.035$ | $0.034$ |
| $500$ | $0.432$ | $0.091$ | $0.039$ | $0.001$ | $0.054$ | $0.014$ |  $0.014$ |
| $50$ | $0.074$ | $0.038$ | $0.091$ | $0.045$ | $0.116$ | $0.024$ | $0.023$ |
group lasso | $100$ | $0.062$ | $0.026$ | $0.073$ | $0.004$ | $0.130$ | $0.029$ | $0.028$ |
| $500$ | $0.020$ | $0.006$ | $0.042$ | $0.001$ | $0.067$ | $0.018$ | $0.018$ |
Scenario 2 ($q=20$, block sparsity, $76.5\%$ zeroes)
| n | TPR | FPR | $||\widehat{\mathbf{R}}-\mathbf{R}||_{\text{F}}/q$ | $\widehat{\mathcal{D}}_{t\ln(t)}$ | $\widehat{\mathcal{D}}_{(\sqrt{t}-1)^{2}}$ | $\widehat{\mathcal{D}}_{1}$ | $\widehat{\mathcal{D}}_{2}$ |
| $50$ | $\backslash$ | $\backslash$ | $0.138$ | $0.006$ | $0.178$ | $0.065$ | $0.064$ |
no penalty | $100$ | $\backslash$ | $\backslash$ | $0.097$ | $0.004$ | $0.077$ | $0.030$ | $0.030$ |
| $500$ | $\backslash$ | $\backslash$ | $0.043$ | $0.001$ | $0.016$ | $0.006$ |  $0.006$ |
| $50$ | $0.940$ | $0.340$ | $0.058$ | $0.108$ | $0.246$ | $0.026$ | $0.026$ |
lasso | $100$ | $0.903$ | $0.135$ | $0.033$ | $0.020$ | $0.139$ | $0.013$ | $0.013$ |
| $500$ | $0.888$ | $0.048$ | $0.013$ | $0.004$ | $0.048$ | $0.005$ |  $0.005$ |
| $50$ | $0.858$ | $0.431$ | $0.101$ | $0.010$ | $0.103$ | $0.024$ | $0.024$ |
adaptive lasso | $100$ | $0.860$ | $0.329$ | $0.076$ | $0.004$ | $0.053$ | $0.013$ | $0.013$ |
| $500$ | $0.900$ | $0.164$ | $0.032$ | $0.001$ | $0.013$ | $0.003$ |  $0.003$ |
| $50$ | $0.940$ | $0.340$ | $0.058$ | $0.108$ | $0.246$ | $0.026$ | $0.026$ |
scad | $100$ | $0.906$ | $0.137$ | $0.033$ | $0.021$ | $0.142$ | $0.013$ | $0.013$ |
| $500$ | $0.846$ | $0.051$ | $0.016$ | $0.005$ | $0.062$ | $0.006$ |  $0.006$ |
| $50$ | $0.889$ | $0.011$ | $0.035$ | $0.021$ | $0.138$ | $0.013$ | $0.013$ |
group lasso | $100$ | $0.934$ | $0.012$ | $0.021$ | $0.011$ | $0.108$ | $0.010$ | $0.010$ |
| $500$ | $0.965$ | $0.013$ | $0.008$ | $0.003$ | $0.044$ | $0.004$ |  $0.004$ |

Scenario 3 ($q=70$, block sparsity, $79.7\%$ zeroes)
---
| n | TPR | FPR | $||\widehat{\mathbf{R}}-\mathbf{R}||_{\text{F}}/q$ | $\widehat{\mathcal{D}}_{t\ln(t)}$ | $\widehat{\mathcal{D}}_{(\sqrt{t}-1)^{2}}$ | $\widehat{\mathcal{D}}_{1}$ | $\widehat{\mathcal{D}}_{2}$ |
no penalty | $100$ | $\backslash$ | $\backslash$ | $0.099$ | $2.344\cdot 10^{-9}$ | $0.021$ | $0.138$ | $0.125$ |
$500$ | $\backslash$ | $\backslash$ | $0.044$ | $2.239\cdot 10^{-9}$ | $0.007$ | $0.023$ |  $0.022$ |
lasso | $100$ | $0.993$ | $0.757$ | $0.048$ | $0.532$ | $0.629$ | $0.073$ | $0.072$ |
$500$ | $0.874$ | $0.087$ | $0.015$ | $2.520\cdot 10^{-8}$ | $0.015$ | $0.007$ |  $0.008$ |
adaptive lasso | $100$ | $0.955$ | $0.740$ | $0.058$ | $0.001$ | $0.240$ | $0.033$ | $0.033$ |
$500$ | $0.838$ | $0.323$ | $0.038$ | $1.887\cdot 10^{-9}$ | $0.005$ | $0.012$ |  $0.011$ |
scad | $100$ | $0.993$ | $0.757$ | $0.048$ | $0.532$ | $0.629$ | $0.073$ | $0.072$ |
$500$ | $0.883$ | $0.092$ | $0.015$ | $3.969\cdot 10^{-8}$ | $0.019$ | $0.009$ |  $0.009$ |
group lasso | $100$ | $0.914$ | $0$ | $0.023$ | $4.331\cdot 10^{-8}$ | $0.013$ | $0.009$ | $0.009$ |
$500$ | $0.970$ | $0$ | $0.007$ | $8.943\cdot 10^{-9}$ | $0.009$ | $0.004$ |  $0.003$ |

estimates is rather low in case no penalty is used, which is again an effect
of the dimension. Aside from this, the group lasso or the adaptive lasso
(especially for larger sample sizes) gives the lowest error.

In Scenario 3 ($q=70$, block sparsity, $79.7\%$ zeroes), the group lasso is
again desirable for exploiting the sparsity structure, and even achieves zero
FPR. For accurate estimation of $\mathbf{R}$ in Frobenius norm, the group
lasso also performs best. Most accurate estimation of dependence (except for
the mutual information, where no penalty performs best because the true value
is again very close to one) is obtained by the group lasso when $n=100$, and
by the adaptive lasso when $n=500$ is large compared to $q$.

In conclusion, when the true correlation matrix is sparse, interpretability
can be enhanced by using a penalty that is able to completely shrink entries
to zero. The group lasso is preferred for obtaining both good TPR and FPR when
this sparsity is at the block level (Scenarios 2 and 3), and also performs
well in estimating dependence in such cases, particularly when $n$ is rather
small. The true $\mathbf{R}$ is estimated more accurately (compared to using
no penalty) in Frobenius norm when using any penalty and in any scenario, but
this does not necessarily result in better estimation of dependence.
Especially when there are still quite some non-zeroes (Scenario 1), but also
when $n$ is rather large, adaptive lasso is recommended for good accuracy of
the estimated dependence coefficients. The $\Phi$-dependence measures are
estimated with rather low error when using no penalty, but this is because
they attain their upper bound of $1$ rather quickly when the dimension
increases.

##  6 Real data applications

In Section 6.1, we look into an application of optimal transport dependence
measures between possibly more than two random vectors to sensory analysis. In
Section 6.2, we illustrate how these measures, together with the considered
penalization techniques, can be useful in finding clusters among speech signal
attributes used for detecting Parkinson’s disease.

###  6.1 Application to sensory analysis

Consider a caterer who wants to sell eight different smoothies, say
$S_{1},\dots,S_{8}$, on an event, and is looking for three employees willing
to take up this job. A total of $24$ people, say $P_{1},\dots,P_{24}$, show up
for this job opportunity, all equally qualified, and the caterer is looking
for a fair way to pick three candidates. Every candidate is asked to taste
each of the eight smoothies, and is given a sheet of paper in order to
position the different smoothies, knowing that the closer certain smoothies
are to each other, the more similar they are considered by the individual. For
example, according to candidate $P_{j}$ in Fig. 8, the smoothies
$S_{1},S_{2},S_{3}$ and $S_{4}$ are similar, but quite different from the
similarly tasting smoothies $S_{5},S_{6}$ and $S_{7}$, and none of them
resembles $S_{8}$.

Fig. 8: Example of smoothie similarity rating by candidate $P_{j}$ on a sheet
of paper and corresponding dataset $\mathcal{X}_{j}$.

As such, the caterer acquires $24$ datasets, say
$\boldsymbol{\mathcal{X}}_{j}\in\mathbb{R}^{8\times 2}$ for $j=1,\dots,24$,
containing the smoothies as rows and the $X$-$Y$ coordinates on the sheet of
paper for person $j$, denoted as a random vector $(X_{j},Y_{j})$, as columns,
representing the smoothie similarities of each candidate. So, the dataset
$\boldsymbol{\mathcal{X}}_{j}$ contains a sample from $(X_{j},Y_{j})$ of size
eight, denoted as $(x_{ij},y_{ij})$ for $i=1,\dots,8$. The data is available
in the R package SensoMineR ([30]). The criterion based on which three
employees are picked consists of finding the three individuals that have the
least similar spatial configurations, meaning three very diversified tastes
(in the hope of not selling only a few smoothies because of prepossessed
preferences by the sellers).

Typically, see, e.g., [32], the similarity between two configurations is
measured by the RV coefficient

$\text{RV}(\boldsymbol{\mathcal{X}}_{j_{1}},\boldsymbol{\mathcal{X}}_{j_{2}})=\frac{\text{tr}\left(\boldsymbol{\mathcal{X}}_{j_{1}}\boldsymbol{\mathcal{X}}_{j_{1}}^{\text{T}}\boldsymbol{\mathcal{X}}_{j_{2}}\boldsymbol{\mathcal{X}}_{j_{2}}^{\text{T}}\right)}{\sqrt{\text{tr}\left\\{\left(\boldsymbol{\mathcal{X}}_{j_{1}}\boldsymbol{\mathcal{X}}_{j_{1}}^{\text{T}}\right)^{2}\right\\}\text{tr}\left\\{\left(\boldsymbol{\mathcal{X}}_{j_{2}}\boldsymbol{\mathcal{X}}_{j_{2}}^{\text{T}}\right)^{2}\right\\}}},$
(32)

and for three configurations one can take, e.g., the average of all pairwise
RV coefficients. Yet, pairwise coefficients feel unnatural and it would be
better to compute a trivariate vector similarity. Note that (32) is actually
the RV coefficient between two random vectors $(X_{j_{1}},Y_{j_{1}})$ and
$(X_{j_{2}},Y_{j_{2}})$ of size two having joint, empirical covariance matrix

$\begin{pmatrix}\boldsymbol{\mathcal{X}}_{j_{1}}^{\text{T}}\boldsymbol{\mathcal{X}}_{j_{1}}&\boldsymbol{\mathcal{X}}_{j_{1}}^{\text{T}}\boldsymbol{\mathcal{X}}_{j_{2}}\vspace{0.2cm}\\\
\boldsymbol{\mathcal{X}}_{j_{2}}^{\text{T}}\boldsymbol{\mathcal{X}}_{j_{1}}&\boldsymbol{\mathcal{X}}_{j_{2}}^{\text{T}}\boldsymbol{\mathcal{X}}_{j_{2}}\end{pmatrix}\in\mathbb{R}^{4\times
4},$

and we get a similar (larger) block covariance matrix in $\mathbb{R}^{2m\times
2m}$ when taking $m$ individuals into account.

Since not restricted to two random vectors anymore, one can also opt for an
(estimated) optimal transport dependence coefficient $\mathcal{D}_{1}$ or
$\mathcal{D}_{2}$ between three vectors of size two for measuring the
similarity between three individual spatial configurations. Note that the
dispersion of the coordinates might differ among the individuals (some might
use the entire sheet, while others only use the right corner), but since we
use the normal scores rank correlations, we do not need any centering or
scaling of the data. Pairwise scatterplots of the normal scores of the $X$-$Y$
coordinates of the $24$ people are shown in Fig. S3 of the Supplementary
Material, indicating that dependencies are mainly correlation based, i.e., a
Gaussian copula model is suitable for modelling the dependencies. When zooming
in on candidates $P_{12},P_{13},P_{18}$ and $P_{20}$, we get the pairwise
scatterplots given in Fig. 9. From this, we expect for example that candidates
$P_{12}$ and $P_{13}$ have quite independent smoothies preferences, while
candidates $P_{18}$ and $P_{20}$ have rather strong correlations between their
coordinates.

Table 2: Arrangement of two and three candidates according to largest estimated similarity $\widehat{\mathcal{D}}_{r}$ for $r=1,2$. The first two and last two are shown. two based on $\widehat{\mathcal{D}}_{1}$ | two based on $\widehat{\mathcal{D}}_{2}$ | three based on $\widehat{\mathcal{D}}_{1}$ | three based on $\widehat{\mathcal{D}}_{2}$
---|---|---|---
1-2 3-4 5-6 7-8  $\widehat{\mathcal{D}}_{1}(P_{18},P_{20})$ | $0.561$ | $\widehat{\mathcal{D}}_{2}(P_{18},P_{20})$ | $0.561$ |  $\widehat{\mathcal{D}}_{1}(P_{15},P_{18},P_{20})$ | $0.585$ | $\widehat{\mathcal{D}}_{2}(P_{15},P_{18},P_{20})$ | $0.595$
$\widehat{\mathcal{D}}_{1}(P_{9},P_{23})$ | $0.542$ | $\widehat{\mathcal{D}}_{2}(P_{15},P_{20})$ | $0.521$ | $\widehat{\mathcal{D}}_{1}(P_{9},P_{10},P_{23})$ | $0.550$ | $\widehat{\mathcal{D}}_{2}(P_{10},P_{18},P_{23})$ | $0.561$
⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮
$\widehat{\mathcal{D}}_{1}(P_{12},P_{19})$ | $0.020$ |  $\widehat{\mathcal{D}}_{2}(P_{12},P_{19})$ | $0.016$ | $\widehat{\mathcal{D}}_{1}(P_{2},P_{3},P_{14})$ | $0.080$ | $\widehat{\mathcal{D}}_{2}(P_{2},P_{12},P_{19})$ | $0.075$
$\widehat{\mathcal{D}}_{1}(P_{12},P_{13})$ | $0.015$ | $\hskip 2.27626pt\widehat{\mathcal{D}}_{2}(P_{12},P_{13})$ | $0.013$ |  $\widehat{\mathcal{D}}_{1}(P_{12},P_{13},P_{21})$ | $0.069$ |  $\widehat{\mathcal{D}}_{2}(P_{12},P_{13},P_{21})$ | $0.074$

We denote $\widehat{\mathcal{D}}_{r}(P_{j_{1}},P_{j_{2}})$ and
$\widehat{\mathcal{D}}_{r}(P_{j_{1}},P_{j_{2}},P_{j_{3}})$ for the estimated
similarity between two candidates $P_{j_{1}},P_{j_{2}}$ or three candidates
$P_{j_{1}},P_{j_{2}},P_{j_{3}}$ for $r=1,2$. Recall that these are actually
dependencies between $2$ and $3$ random vectors of size $2$ respectively,
i.e., $q=4$ with $k=2,d_{1}=d_{2}=2$, or $q=6$ with $k=3,d_{1}=d_{2}=d_{3}=2$
respectively, estimated based on a sample of size $n=8$.

Table 2 shows the two strongest and two weakest couples or triplets according
to $\widehat{\mathcal{D}}_{1}$ or $\widehat{\mathcal{D}}_{2}$. The caterer
will definitely hire candidate $P_{12},P_{13}$ and $P_{21}$. In [32], they
cluster the candidates using the clustatis method, bringing forward three
classes of individuals, see their Fig. 5. We see that $P_{12}$ belongs to
class $2$, while $P_{13}$ belongs to class $3$ and $P_{21}$ to class $1$, also
indicating their diversified smoothie similarity pattern. The optimal
transport dependence measures give an unequivocal ordering of patterns based
on similarities that go beyond two random vectors.

Note that we can also switch the role of the smoothies and the individuals,
i.e., construct $8$ datasets in $\mathbb{R}^{24\times 2}$, and similarly look
at the dependence between smoothies. Doing so, both
$\widehat{\mathcal{D}}_{1}$ $(=0.0336)$ and $\widehat{\mathcal{D}}_{2}$
$(=0.0340)$ agree that $S_{6},S_{7}$ and $S_{8}$ are the least similar among
all possible triplets. They are respectively called Casino_PBC, Innocent_SB
and Carrefour_SB, and the biplots (based on various sensometrics methods) in
Fig. 2. of [32] indeed also reveal angels between these smoothies that are
close to $90^{\circ}$. So, $S_{6},S_{7}$ and $S_{8}$ would be a good choice if
the caterer wanted to limit his smoothie supply to three flavours that still
have a satisfactory amount of diversity.

Fig. 9: Pairwise scatterplots of normal scores of $X$-$Y$ coordinates of
candidates $P_{12},P_{13},P_{18}$ and $P_{20}$ of the smoothies dataset.

Fig. 10: Boxplot of $p$-values of pairwise goodness-of-fit tests for Gaussian
copula on $91$ dysphonia measures of the LSVT dataset.

###  6.2 Clustering dysphonia measures

In a second real data application, we study the LSVT voice rehabilitation
dataset of [47], freely accessible at the UCI Machine Learning Repository
(https://archive.ics.uci.edu/dataset/282/lsvt+voice+rehabilitation).
Parkinson’s disease frequently leads to vocal impairment, the extent of which
can be assessed using sustained vowel phonations. In particular, the sustained
vowel “ahh…” (denoted as /a/) is typically studied. Next, dysphonia measures
are used to extract clinically useful information from speech signals. The
dataset consists of $q=310$ such measurements on $n=126$ phonations. As
mentioned in Section C of [47], several of these dysphonia measures are
similar (algorithmic variations of the same basic ideas), e.g., there are many
jitter (a measure of cycle-to-cycle variation in frequency) and shimmer (a
measure of cycle-to-cycle variation in amplitude) variants. Hence, there is a
large amount of redundancy among the attributes, which could worsen the
performance of supervised learning (like, e.g., a classifier as studied in
[47]), and some kind of preliminary feature selection is recommended.

There are indeed many dysphonia measures that exhibit a (very) strong sample
normal scores rank correlation, see Fig. S4 of the Supplementary Material. In
order to eliminate strong detrimental redundancies (that closely approach
singularity), we iteratively search (pairwise) for attributes that have a
normal scores rank correlation larger than $0.8$ in absolute value, and each
time discard one of them. Remaining are $91$ dysphonia measures, whose (non-
singular) sample normal scores rank correlation matrix is given in the left
panel of Fig. S5 of the Supplementary Material.

For testing the adequacy of a Gaussian copula for modelling (at least
pairwise) dependencies, we consider the test based on the statistic $S_{n}$
given in equation (2) of [21], where we test the null hypothesis that the
copula between a pair of dysphonia measures is Gaussian, and compute the
$p$-value based on $1000$ bootstrap samples. In Fig. 10, a boxplot of
$p$-values of pairwise (that is between all $91\cdot 90/2=4095$ pairs of
variables) goodness-of-fit tests is shown. Based on the boxplot, we see that
is it reasonable to assume (at least pairwise) Gaussian dependencies.

We still have quite a large amount of attributes compared to the sample size
of $n=126$, and we still see some quite large correlations between several
first and last dysphonia measures (right upper corner of left panel in Fig.
S5), so further feature selection is desired. With this objective in mind, we
decide to cluster the remaining $91$ attributes in order to get an idea about
how they can be divided into properly separated groups that have rather strong
similarity within.

Fig. 11: Dendrogram (left) and redundancy (right) for the clustering of $91$
dysphonia measures.

In particular, we opt for agglomerative hierarchical variable clustering based
on similarity measures (see, e.g., [19] and references within for a survey),
where the main task is to measure the similarity, say
$\mathcal{D}(\mathbb{X}_{1},\mathbb{X}_{2})$ (or alternatively dissimilarity),
between two groups of variables $\mathbb{X}_{1}$ and $\mathbb{X}_{2}$. Denote
$X_{1},\dots,X_{91}$ for the $91$ dysphonia measures, which are initially
considered to be clusters on their own. In the first step, the two variables
$X_{i}$ and $X_{j}$ (for certain $i,j\in\\{1,\dots,91\\}$ with $i\neq j$) that
exhibit the largest similarity according to
$\mathcal{D}(\\{X_{i}\\},\\{X_{j}\\})$, i.e., taking
$\mathbb{X}_{1}=\\{X_{i}\\}$ and $\mathbb{X}_{2}=\\{X_{j}\\}$, are merged
together forming one single cluster. Next, all similarities between the
current clusters are again computed, and the two clusters showing the largest
similarity are merged. One keeps repeating this until only one big cluster
(composed of all the $91$ attributes) remains, yielding a total of $91$
partitions of the attribute space. In general, the main task is thus to
compute (for $\mathbb{X}_{1}=\\{X_{i_{1}},\dots,X_{i_{m}}\\}$ and
$\mathbb{X}_{2}=\\{X_{j_{1}},\dots,X_{j_{r}}\\}$)

$\mathcal{D}\left(\mathbb{X}_{1},\mathbb{X}_{2}\right)=\mathcal{D}\left(\\{X_{i_{1}},\dots,X_{i_{m}}\\},\\{X_{j_{1}},\dots,X_{j_{r}}\\}\right)$
(33)

for certain mutually exclusive
$i_{1},\dots,i_{m},j_{1},\dots,j_{r}\in\\{1,\dots,91\\}$. Typically, one
specifies a certain (estimated) bivariate dependence coefficient and a certain
link function (overlooking multivariate dependence structures) for computing
(33), or (estimated) multivariate concordance measures (between univariate
random variables) have also been studied (see, e.g., [19]). We decide to take
$\widehat{\mathcal{D}}_{1}((X_{i_{1}},\dots,X_{i_{m}}),(X_{j_{1}},\dots,X_{j_{r}}))$
for (33), and as such measure the similarity between clusters (in spite of
possible similarity within), which feels more natural and does not require a
link function. This $\widehat{\mathcal{D}}_{1}$ is taken to be
$\mathcal{D}_{1}(\widehat{\mathbf{R}}_{\text{R},n})$, where
$\widehat{\mathbf{R}}_{\text{R},n}$ is the ridge penalized estimated Gaussian
copula correlation matrix (23), keeping in mind that the total dimension $m+r$
in (33) might be large compared to $n$. The tuning parameter $\omega_{n}$ is
again determined by a $5$-fold cross-validation search on a grid of $50$
equidistant elements in $[0.01,0.999]$.

We obtain $91$ partitions, and one of these should be picked as preferable
clustering of the attributes. Since not restricted to two groups, we can also
use $\widehat{\mathcal{D}}_{1}$ for measuring the similarity
$\mathcal{D}(\mathbb{X}_{1},\dots,\mathbb{X}_{k})$ between $k$ clusters
$\mathbb{X}_{1},\dots,\mathbb{X}_{k}$ comprising a specific partition. This in
fact measures how separated the obtained clusters of that particular partition
are, i.e., it is a measure of redundancy among the clusters of attributes,
which we want to be rather small. The dendrogram of the clustering procedure
and redundancy $\mathcal{D}(\mathbb{X}_{1},\dots,\mathbb{X}_{k})$ as a
function of the number of clusters $k$ (for each of the $91$ obtained
partitions) are shown in Fig. 11. Based on this, we decide to look deeper into
the $36$ cluster partition (where the dotted line cuts the dendrogram), since
the redundancy is minimal here (in particular,
$\mathcal{D}(\mathbb{X}_{1},\dots,\mathbb{X}_{36})=0.07265$). Next, we
rearrange the attributes such that dysphonia measures that belong to the same
cluster follow each other in the dataset, and the cluster dimensions are

$(d_{1},\dots,d_{36})=(1,1,53,1,1,1,1,1,1,1,2,1,1,1,1,1,1,1,1,1,1,1,1,1,3,1,1,1,1,1,1,1,1,1,1,1),$

i.e., there is a big cluster of size $53$, two smaller clusters of size $2$
and $3$, and $33$ clusters of size $1$. The sample normal scores rank
correlation matrix after rearrangement, estimated without a penalty and with a
ridge penalty are respectively given in the middle and right panel of Fig. S5
in the Supplementary Material.

Denote now the $36$ clusters as random vectors
$\mathbf{X}_{1},\dots,\mathbf{X}_{36}$. We already know that
$\widehat{\mathcal{D}}_{1}(\mathbf{X}_{1},\dots,\mathbf{X}_{36})=0.07265$ when
using ridge penalization ($\omega_{n}=0.696$), which is already way closer to
$0$ than when using no penalization (then
$\widehat{\mathcal{D}}_{1}(\mathbf{X}_{1},\dots,\mathbf{X}_{36})=0.19013$).
Moreover, the rationale behind variable clustering is that clusters are truly
separated, i.e., there is (block) sparsity in the Gaussian copula correlation
matrix at levels corresponding to attributes belonging to different clusters.

Fig. 12: Sparsity patterns of estimated Gaussian copula correlation matrix of
$91$ dysphonia measures using different penalties.

Therefore, we also consider lasso, adaptive lasso, and group lasso estimation
of the Gaussian copula correlation matrix. For the lasso and group lasso, we
search for an optimal $\omega_{n}$ on an equidistant grid of $50$ values in
$[0.01;0.8]$, and for the adaptive lasso, we search for an optimal $\rho_{n}$
on an equidistant grid of $50$ values in $[0.01;0.6]$, each time aiming to
maximize the BIC (31).

The first panel in Fig. 12 shows the sparse structure of the estimated normal
scores rank correlation matrix of the $91$ clustered attributes when using the
lasso. To a certain extent, we recognize the diagonal blocks corresponding to
the within cluster correlations, and observe many zeroes between attributes
belonging to different clusters. The estimated inter-cluster dependence equals
$\widehat{\mathcal{D}}_{1}(\mathbf{X}_{1},\dots,\mathbf{X}_{36})=0.0014$.
Using the adaptive lasso (second panel), we get
$\widehat{\mathcal{D}}_{1}(\mathbf{X}_{1},\dots,\mathbf{X}_{36})=0.0141$, and
a group lasso estimator (third panel) results in
$\widehat{\mathcal{D}}_{1}(\mathbf{X}_{1},\dots,\mathbf{X}_{36})=0.0130$. Note
that the group lasso estimator did not shrink any of the diagonal blocks,
reflecting the stronger intra-cluster similarities. In the fourth panel of
Fig. 12, we visualize a group lasso estimate with a larger penalty parameter
(not determined through BIC). When forcing more sparsity by increasing
$\omega_{n}$, we see more off-diagonal blocks (and not diagonal blocks)
getting shrunk, and this leads to
$\widehat{\mathcal{D}}_{1}(\mathbf{X}_{1},\dots,\mathbf{X}_{36})=0.0001$. More
detailed images are given in Fig. S6 of the Supplementary Material.

Further feature selection can now be done via replacing each cluster by a
single (or maybe multiple) attribute that represents the cluster well (e.g.,
the attribute within the cluster that has largest similarity with the
cluster).

##  7 Discussion

In this paper, we proposed to use the $2$-Wasserstein distance between the
joint copula measure and the product of the marginal copula measures for
quantifying dependence between a finite, arbitrary amount of random vectors.
The obtained dependence coefficients satisfy several desirable properties and
especially have the powerful theoretical quality of detecting any departure
from independence. Examples illustrate that the choice of normalization
strongly influences the overall dependence quantification. Whether or not one
can explicitly calculate the infimum for the optimal transport map and/or the
supremum required for the normalization, depends on the specific form of the
copula, and a great deal of interesting work remains to be done (also when no
explicit copula can be assumed, i.e., when nonparametric estimators need to be
considered).

A Gaussian copula approach yields explicit formulas with a clear
interpretation. Using the sample matrix of normal scores rank correlation
coefficients results in an easily computable plug-in estimator for which we
obtained an asymptotic normality result with explicit asymptotic variance in
arbitrary dimensions. Expectedly, higher dimensions aggravate the finite
sample estimation performance. To cope with this, we studied rank-invariant
penalization techniques for estimating the Gaussian copula correlation matrix,
leading to estimators that are able to improve accuracy on the one hand, and
enhance interpretability by detecting marginal independencies on the other
hand. Such estimation challenges have enjoyed rather little attention, and
further research would definitely be worthwhile.

Another interesting theoretical challenge is to study the optimal transport
(and others as well) dependence coefficients when the dimension grows
unboundedly. In the Gaussian copula context, random matrix theory could
definitely be useful, and, as also touched upon in our simulations, various
behaviours can be expected depending on the nature of the dependence measure,
the normalization, whether letting $d_{i}\to\infty$ for some $i$, or
$k\to\infty$. Keeping the dimension fixed, our simulations illustrated the
asymptotic normality result and the benefits of using penalization techniques
when the sample size is rather small and/or when (group) sparsity is pursued.

Finally, in a first real data application, we illustrated the use of the
dependence coefficients in evaluating and comparing (possibly more than two)
consumer products or similarities in sensory analysis. Alongside, on a second
real dataset containing expressions on sustained vowel phonations, we
demonstrated how attributes can be hierarchically clustered via multivariate
similarities between random vectors (despite similarities within), disposing
of traditional link functions. Ridge penalization is preferred when the number
of attributes is large compared to the sample size, and (group)-lasso type
penalties can be used to reflect the homogeneity and separation of a partition
of, in general $k$, groups of variables.

Acknowledgments. The authors thank Dr Gilles Mordant, Georg-August-Universität
Göttingen, for scientific discussions during the startup phase of this
research. The authors gratefully acknowledge support from the Research Fund KU
Leuven [C16/20/002 project].

##  Appendix A Axioms for dependence measures between random vectors

In [11], a list of axioms is stated for a dependence measure
$\mathcal{D}^{d_{1}\dots,d_{k}}(\mathbf{X})=\mathcal{D}(\mathbf{X}_{1},\dots,\mathbf{X}_{k})$.
Up to some minor differences (small corrections and simplifications), the
axioms are given as follows.

  1. (A1)

For every permutation $\\{\pi(1),\dots,\pi(k)\\}$ of $\\{1,\dots,k\\}$:
$\mathcal{D}^{d_{1},\dots,d_{k}}(\mathbf{X})=\mathcal{D}(\mathbf{X}_{\pi(1)},\dots,\mathbf{X}_{\pi(k)})$:
and for every permu- tation $\\{\pi_{i}(1),\dots,\pi_{i}(d_{i})\\}$ of
$\\{1,\dots,d_{i}\\}$, $i\in\\{1,\dots,k\\}$:
$\mathcal{D}^{d_{1},\dots,d_{k}}(\mathbf{X})=\mathcal{D}(\mathbf{X}_{1},\dots,(X_{i\pi_{i}(1)},\dots,X_{i\pi_{i}(d_{i})}),\dots,\mathbf{X}_{k})$.

  2. (A2)

$0\leq\mathcal{D}^{d_{1},\dots,d_{k}}(\mathbf{X})\leq 1$.

  3. (A3)

$\mathcal{D}^{d_{1},\dots,d_{k}}(\mathbf{X})=0$ if and only if
$\mathbf{X}_{1},\dots,\mathbf{X}_{k}$ are mutually independent.

  4. (A4)

$\mathcal{D}(\mathbf{X}_{1},\dots,\mathbf{X}_{k},\mathbf{X}_{k+1})\geq\mathcal{D}(\mathbf{X}_{1},\dots,\mathbf{X}_{k})$
with equality if and only if $\mathbf{X}_{k+1}$ is independent of
$(\mathbf{X}_{1},\dots,\mathbf{X}_{k})$.

  5. (A5)

$\mathcal{D}^{d_{1},\dots,d_{k}}(\mathbf{X})$ is well-defined for any
$q$-dimensional random vector $\mathbf{X}$ (even if there is a singular part
in the distribution of $\mathbf{X}$).

  6. (A6)

$\mathcal{D}^{d_{1},\dots,d_{k}}(\mathbf{X})$ is a function of solely the
copula $C$ of $\mathbf{X}$ (which is equivalent to
$\mathcal{D}^{d_{1},\dots,d_{k}}(\mathbf{X})$ being invariant under strictly
increasing transformations of any of the components of $\mathbf{X}$).

  7. (A7)