inkscapepath=i/svg-inkscape/ svg/

# FastGCL: Fast Self-Supervised Learning on Graphs via
Contrastive Neighborhood Aggregation

Yuansheng Wang1 Wangbin Sun1 Kun Xu1 Zulun Zhu2 Liang Chen1 Zibin Zheng1 1Sun
Yat-sen University
2Nanyang Technological University {wangysh29, sunwb7,
<EMAIL_ADDRESS><EMAIL_ADDRESS>
{chenliang6<EMAIL_ADDRESS>
###### Abstract

Graph contrastive learning (GCL), as a popular approach to graph self-
supervised learning, has recently achieved a non-negligible effect. To achieve
superior performance, the majority of existing GCL methods elaborate on graph
data augmentation to construct appropriate contrastive pairs. However,
existing methods place more emphasis on the complex graph data augmentation
which requires extra time overhead, and pay less attention to developing
contrastive schemes specific to encoder characteristics. We argue that a
better contrastive scheme should be tailored to the characteristics of graph
neural networks (e.g., neighborhood aggregation) and propose a simple yet
effective method named FastGCL. Specifically, by constructing weighted-
aggregated and non-aggregated neighborhood information as positive and
negative samples respectively, FastGCL identifies the potential semantic
information of data without disturbing the graph topology and node attributes,
resulting in faster training and convergence speeds. Extensive experiments
have been conducted on node classification and graph classification tasks,
showing that FastGCL has competitive classification performance and
significant training speedup compared to existing state-of-the-art methods.

##  1 Introduction

Graph self-supervised learning (GSSL), which is capable of learning
representations of graph structured data without artificially crafted labels,
has numerous application scenarios in real systems Wu et al. (2021). Recently,
graph contrastive learning (GCL) Velickovic et al. (2019); Zhu et al. (2020);
You et al. (2020) has made a big splash on GSSL, successfully applying
contrastive learning from the field of computer vision to graph data. By
pulling the relevant pairs (positive) closer and pushing the irrelevant pairs
(negative) farther on the latent space, GCL methods achieve results comparable
to, and sometimes surpassing, supervised methods.

One of the critical stages in contrastive learning is data augmentation, which
aims to provide different views for contrastive pairs construction. Unlike
images, the augmentation of graph data has to take not only the semantic but
also the graph structural information into account. Existing GCL methods
utilize feature shuffling Velickovic et al. (2019), edge addition and removal
You et al. (2020), subgraph sampling Qiu et al. (2020), etc., for graph data
augmentation. The effectiveness of existing methods relies to their carefully
designed data augmentation strategies Lee et al. (2021).

However, these graph augmentation approaches mentioned above pose some
inevitable problems. For one thing, the time consumption is increased when
applying a complex augmentation algorithm. For another thing, it disturbs the
intrinsic structure of the original graph to a certain extent, accompanied
with the risk of semantic drift Sun et al. (2021), and it is difficult to
determine whether the augmented graph is positively correlated with the
original graph Lee et al. (2021). The majority of GCL methods follow the
InfoMax Linsker (1988) principle, and aim to maximize mutual information of
the corresponding graphs (or nodes) under different views of augmentation. To
correct the potential semantic drift introduced by these graph data
augmentations, encoders need to learn an extra nonlinear transformation to
rebuild the mapping between the augmented graph and the original graph, which
leads to the slow convergence of training.

Figure 1: t-SNE embeddings of nodes in the Amazon-Photo dataset. A randomly
initialized untrained SGC model can achieve discriminative representations
(right) relying only on the neighbor aggregation, compared to raw features
(left).

The limitations of existing augmentation strategies for GCL inspire us to look
back to the characteristics of the GNN encoder itself. In Figure 1 we
visualize the raw node features and the embedding of the randomly initialized
untrained SGC Wu et al. (2019) activation layer. It can be observed that
without the help of trained model parameters, a more discriminative graph
representation carrying semantics information can be obtained by neighbor
information aggregation alone. Thus, we argue that a better contrastive scheme
should be tailored to the characteristics of GNNs (e.g., neighborhood
aggregation).

Therefore, we propose a new solution for GCL, named Fast Graph Contrastive
Learnig (FastGCL). Unlike previous methods, FastGCL utilizes contrastive
neighborhood aggregation for self-supervised learning, without disturbing
graph topology or node attributes. Specifically, FastGCL generates different
representations by using a shared GNN encoder. To obtain the positive, it
weights each edge of the original graph for weighted neighborhood aggregation.
Edge weights are computed by learning interactions between node
representations of original graph output by the encoder. To obtain the
negative, it simply choose the graph with only self-loop as input of the
encoder. Our approach alleviates the time overhead of complex data
augmentation strategies, as well as the potential risk of semantic drift from
structural changes in graph data, which contributes to faster training and
convergence speed compared to previous methods.

Our contributions can be summarized as follows:

  * •

We demonstrate the effectiveness of the GNN neighborhood aggregation in the
unsupervised learning and suggest that a better graph contrastive scheme
should be tailored to the characteristics of GNNs.

  * •

We propose FastGCL, a fast and effective graph contrastive method in which the
contrastive neighborhood aggregation is employed without disturbing the graph
topology or node attributes.

  * •

Comprehensive experiments on diverse datasets demonstrate that FastGCL is able
to achieve competitive results with a remarkable training speedup compared to
the state-of-the-art GCL methods.

##  2 Related Work

Figure 2: Overview of FastGCL. It mainly consists of a shared GNN encoder and
a learnable edge-weighting module, and the latter takes the node
representation $\mathbf{H}^{\alpha}$ of the original graph $\mathcal{G}$ as
input (Red line). The positive graph $\mathcal{G}^{\rho}$ is generated by the
edge-weighting module, while the $\mathcal{G}^{\eta}$ with only self-loop
(mathematically, $\mathbf{A}=\mathbf{I}$) is straightforwardly considered as
the negative graph. All learnable parameters are optimized via contrastive
loss.

Inspired by the field of computer vision, GCL methods have been explored for
self-supervised learning on graphs with remarkable success. The most
substantial challenge for GCL is how to develop contrastive pairs on abstract
and irregular graph data. The general approach is performing data augmentation
to generate different graph views. If the augmented graphs (or nodes)
originate from the same graph (or node), they are constructed as positive
sample pairs (semantic correlation), otherwise they are negative sample pairs.
According to the strategies employed in current studies, we can categorize
them into the following three groups:

#### Uniform Graph Data Augmentation.

This strategy focuses on the structure and attributes of graph data and
implements a non-discriminatory augmentation for all study objects (i.e.,
graphs or nodes). For example, DGI Velickovic et al. (2019) generates negative
graphs by corruption functions (e.g., random shuffling of features in row-
wise). GRACE Zhu et al. (2020) facilitates the optimization of contrastive
objectives from two augmentation views through random edge removal and feature
masking. GCC Qiu et al. (2020) uses ego-net sampling to contrst the L-hop
neighbor subgraph of nodes. GraphCL You et al. (2020) parameterizes four graph
augmentation strategies (i.e., node dropping, edge perturbation, attribute
masking, and subgraph sampling) to select best combination of augmentation.

#### Guided Graph Data Augmentation.

Uniform graph data augmentation may perturb important nodes or edges and then
affect the positive correlation of the same object in different views.
Therefore, some guided approaches emerge to identify the importance of nodes
or edges through domain knowledge or learnable strategies to keep the
semantics inherent in the augmented views. GCA Zhu et al. (2021) proposes
adaptive data augmentation, specifically highlighting the important links
based on node centrality, and adding noise to unimportant nodes. MoCL Sun et
al. (2021) combines domain knowledge for molecular graph augmentation, without
altering the molecular properties too much. AD-GCL Suresh et al. (2021), which
extends the graph information bottleneck Wu et al. (2020) to unsupervised
tasks, guides the process of graph data augmentation by a learnable
perturbation model for edge removals.

#### No Negative Sampling or Augmentation.

Instead of using negative sampling or data augmentation, some works looked for
other information to construct contrastive pairs, which inspired our work to
consider the characteristics of the encoder. GMI Peng et al. (2020), InfoGraph
Sun et al. (2019) do not require explicit data augmentation, the former
maximizes the mutual information of edge and node features on the original
graph, and the latter constructs contrastive pairs in a batch-wise fashion.
BGRL Thakoor et al. (2021) is based on the BYOL Grill et al. (2020) framework
without the use of negative samples, which avoids the sampling bias Bielak et
al. (2021). AFGRL Lee et al. (2021) discovers nodes that can be treated as
positive samples by considering the local structural information and global
semantics of the graph, without creating an augmented graph view and negative
sampling.

Most of the aforementioned methods, however, disturb the structural and
attribute information when providing different views for graphs. In contrast,
our approach tends to preserve topology and attributes. We use a learnable
edge weighting module to provide the positive graph view and preserve semantic
information to the maximum extent. The design for the negative sample is
straightforwardly based on GNNs’ aggregation mechanism, which avoids a large
amount of negative samples and sampling bias. Thus, by considering the
difference of neighbor information aggregation as a contrastive objective, our
augmentation strategy is more suitable for GNNs.

##  3 Proposed Method

###  3.1 Notations and Preliminaries

Let $\mathcal{G}=(\mathcal{V},\mathcal{E})$ denote an undirected graph, where
$\mathcal{V}=\\{v_{1},v_{2},\cdots,v_{n}\\}$ represents all nodes,
$N=|\mathcal{V}|$, and $\mathcal{E}$ is the set of all edges,
$M=|\mathcal{E}|$. The matrix $\mathbf{X}\in\mathbb{R}^{N\times F}$ is
composed of $F$-dimensional attributes $x_{v}$ of node $v$ on graph
$\mathcal{G}$. The symmetric matrix $\mathbf{A}\in\mathbb{R}^{N\times N}$
denotes the adjacency matrix where $A_{ij}=1$, if
$(v_{i},v_{j})\in\mathcal{E}$, otherwise $A_{ij}=0$. In addition, we use the
matrix $\mathbf{E}\in\mathbb{R}^{N\times N}$ to denote the weights on the
corresponding edges of the adjacency matrix, which are floating point numbers
between 0 and 1. $\mathcal{N}_{v}$ is the set of neighbors of the node $v$ and
$\mathbf{D}$ is the degree matrix ${D}_{ii}=\sum^{N}_{j=0}A_{ij}$. We denote
the symmetric renormalized adjacency matrix $\mathbf{\hat{A}}$ by
$\mathbf{\hat{A}}=\mathbf{\tilde{D}^{1/2}\tilde{A}\tilde{D}^{1/2}}$, where
$\mathbf{\tilde{A}}=\mathbf{A}+\mathbf{I}$ is the adjacency matrix with self-
loop and $\mathbf{I}$ is the identity matrix.

In our work, we use GNN as our encoder. It typically uses a recursive
neighborhood aggregation protocol where the representation $\mathbf{h}_{v}$ of
node $v$ is obtained by aggregating its neighborhood representations by
$agg(\cdot)$ at each iteration, which is then updated by a nonlinear function
$upd(\cdot)$. The representation of a node can capture the semantic and
structural information of its $k$-hop neighbors after $k^{th}$ iterations of
aggregation. The process can be expressed as,

$\mathbf{h}_{v}^{(k)}=\textit{upd}^{(k)}\left(\mathbf{h}_{v}^{(k-1)},\textit{
agg}^{(k)}\left(\left\\{\mathbf{h}_{u}^{(k-1)}\mid
u\in\mathcal{N}_{v}\right\\}\right)\right)$ (1)

where the $\mathbf{x}_{v}$ is as the initialization $\mathbf{h}_{v}^{(0)}$ of
the node representation.

Our goal is to train an encoder $f_{\theta}(\cdot)$, given a graph
$\mathcal{G}$ along with $\mathbf{X}$ and $\mathbf{A}$, to learn a high-level
node representation $\mathbf{H}=f_{\theta}(\mathbf{A},\mathbf{X})$ which can
be generalized to various downstream tasks, without using any label
information of categories.

###  3.2 FastGCL

The overall framework of FastGCL is shown in Figure 2. The generation of graph
views and the objective function are two vital components of GCL methods,
which are introduced in this section.

#### Graph View Generation.

Standard contrastive learning approaches implicitly assume that the positive
and the negative belong to different categories in downstream tasks Saunshi et
al. (2019). Zhu et al. (2021) states that data augmentation strategies should
preserve the intrinsic structure and attributes of the graph. Our design
adheres to the above constraints. Considering the representation of original
graph $\mathcal{G}$ as the anchor, we need to construct the positive and the
negative under two other graph views (i.e., $\mathcal{G}^{\rho}$ and
$\mathcal{G}^{\eta}$) for contrastive learning. This process can be described
as follows.

Firstly, with a GNN encoder $f_{\theta}(\cdot)$ , all node representations of
the original graph $\mathbf{H}^{\alpha}$ can be obtained. Mathematically,

$\mathbf{H}^{\alpha}=f_{\theta}(\mathbf{A},\mathbf{X})$ (2)

$\mathbf{H}^{\alpha}$ aggregates rich information about each node and its
neighbors. We take $\mathbf{H}^{\alpha}$ as input of an edge-weighting module
to generate the weights of the corresponding edges. This process involves the
help of a MLP $g_{\varphi}(\cdot)$, which projects the node representations to
a new latent space, denoted as,

$\mathbf{Z}=g_{\varphi}(\mathbf{H}^{\alpha})$ (3)
$e_{u,v}=\left\\{\begin{matrix}\sigma(\mathbf{z}_{u}\cdot\mathbf{z}_{v}^{\top}),&\textit{
if }(u,v)\in\mathcal{E}\\\ \\\ 0,&\textit{ otherwise. }\end{matrix}\right.$
(4)

where $\sigma$ is sigmoid function and $e_{u,v}$ represents the weight between
node $u$ and node $v$, if there is an edge between them. Then these weights
are performed a element-wise product with the adjacency matrix $\mathbf{A}$ of
the original graph. The positive representation $\mathbf{H}^{\rho}$ is then
obtained by the shared GNN encoder $f_{\theta}(\cdot)$, represented as,

$\mathbf{H}^{\rho}=f_{\theta}(\mathbf{A}\odot\mathbf{E},\mathbf{X})$ (5)

This is equivalent to the fact that we bring these learned edge weights to the
GNN for weigted neighborhood aggregation, instead of equal aggregation for
Equation 2. For the negative, we simply choose the identity matrix
$\mathbf{I}$ as the shared GNN encoder’s input,

$\mathbf{H}^{\eta}=f_{\theta}(\mathbf{I},\mathbf{X})$ (6)

By taking the identity matrix $\mathbf{I}$ as input, the GNN will degenerate
into an MLP, which can only aggregate information from the nodes themselves.
It is suitable to serve as an uncorrelated negative representation.

With this design, for the positive sample, it does not lose neighbor
information during aggregation, but only affects how much or how little
neighbor information each node receives. This weighted aggregation of
neighbors brings a difference in representation that can be seen as an
augmentation, allowing the encoder to better extract positively relevant
features. For the negative sample, it loses all information about its
neighbors, in which case it is likely to be predicted as another category by
the downstream task. As we know that the information of the graph data is
mainly determined by the composition of the node and its neighbors, see
Equation 1.

#### Objective Function.

Inspired by noise-contrastive estimation (NCE) Gutmann and Hyvärinen (2010),
to pull the positive closer and push the negative farther, our self-supervised
learning loss is designed as,

$\mathcal{L}_{\mathrm{ssl}}=-\frac{1}{N}\sum_{v\in\mathcal{V}}\left[\log\sigma(sim^{\rho}_{v})+\log(1-\sigma(sim^{\eta}_{v}))\right]$
(7)

where $sim$ is the similarity score between two representations, and the
positive and negative samples similarity scores of node $v$ can be expressed
as,

$sim^{\rho}_{v}=\ell(\mathbf{h}_{v}^{\alpha},\mathbf{h}_{v}^{\rho}),\,\,sim^{\eta}_{v}=\ell(\mathbf{h}_{v}^{\alpha},\mathbf{h}_{v}^{\eta})$
(8)

where $\ell$ is the cosine similarity function. For the node classification
task, our contrastive objective is the corresponding nodes in different views.
For the graph classification task, our aim is to learn discriminative and
semantic graph representations, so our contrastive objective becomes the graph
representations of different views. The representation of graph $\mathcal{G}$,
$\mathbf{h}_{\mathcal{G}}$, can be obtained by a pooling operation that
combines the representations of all nodes in $K$ iterations. It is formally
defined as,

$\mathbf{h}_{\mathcal{G}}=\textit{pool}\left(\left\\{\mathbf{h}_{v}^{(K)}\mid
v\in\mathcal{V}\right\\}\right)$ (9)

To prevent the positive graph from becoming almost identical to the original
graph, which would lead to an indistinguishable representation compared with
the anchor, a regularization term is applied. Mathematically,

$\mathcal{L}_{\mathrm{norm}}=-\frac{1}{M}\sum_{i=1}^{N}\sum_{j=1}^{N}\log\left(1+\exp\left(\mathbf{A}-\mathbf{E}\right)\right)$
(10)

which makes the $\mathbf{E}$ less close to the $\mathbf{A}$ to a better
learning. In conjunction with Equation 7, we propose an objective function
that is to minimize the loss $\mathcal{L}$ of the entire network, i.e.,

$\mathrm{\mathop{min}\limits_{\theta,\varphi}}\,\mathcal{L}=\mathcal{L}_{\mathrm{ssl}}+\lambda\mathcal{L}_{\mathrm{norm}}$
(11)

where $\lambda$ is a hyper-parameter to regulate the strength of the penalty.

In summary, 1) FastGCL constructs contrastive pairs for the shared GNN encoder
without disturbing the topology and attributes of the graph and learns the
representation by contrasting the neighborhood aggregations. 2) The positive
samples of FastGCL are obtained by a learnable edge-weighting module, which
can guide the weighted graph to maintain a certain amount of information
associated with downstream tasks. The negative examples, on the other hand,
are straightforward in design and do not require the involvement of complex
augmentation functions.

##  4 Experiments

Table 1: Node classification performance (Accuracy). R.I.U. indicates randomly initialized untrained model. Method | WikiCS | Computers | Photo | Co.CS | Co.Physics
---|---|---|---|---|---
GCN Kipf and Welling (2017) | 78.85 ± 0.09 | 86.42 ± 0.09 | 92.71 ± 0.03 | 93.07 ± 0.02 | 95.90 ± 0.01
GAT Veličković et al. (2017) | 78.93 ± 0.05 | 88.30 ± 0.04 | 92.10 ± 0.04 | 92.71 ± 0.04 | 95.91 ± 0.01
SGC Wu et al. (2019) | 78.63 ± 0.05 | 87.36 ± 0.04 | 92.95 ± 0.04 | 93.03 ± 0.02 | OOM
GCN R.I.U. | 77.97 ± 0.09 | 88.16 ± 0.12 | 92.72 ± 0.07 | 93.14 ± 0.04 | 95.72 ± 0.03
GAT R.I.U. | 78.03 ± 0.05 | 85.72 ± 1.21 | 91.69 ± 0.38 | 91.85 ± 0.14 | 95.21 ± 0.04
SGC R.I.U. | 77.44 ± 0.10 | 88.30 ± 0.13 | 92.83 ± 0.07 | 93.17 ± 0.03 | 95.71 ± 0.03
Raw feature Zhu et al. (2021) | 71.98 ± 0.00 | 73.81 ± 0.00 | 78.53 ± 0.00 | 90.37 ± 0.00 | 93.58 ± 0.00
DGI Velickovic et al. (2019) | 77.31 ± 0.10 | 87.09 ± 0.14 | 92.19 ± 0.12 | 92.68 ± 0.06 | OOM
GRACE Zhu et al. (2020) | 78.50 ± 0.19 | 89.53 ± 0.35 | 92.78 ± 0.45 | 91.12 ± 0.02 | OOM
GCA Zhu et al. (2021) | 78.35 ± 0.05 | 87.85 ± 0.31 | 92.53 ± 0.16 | 93.10 ± 0.01 | 95.73 ± 0.03
BGRL Thakoor et al. (2021) |  79.98 ± 0.10 |  90.34 ± 0.19 | 93.17 ± 0.30 |  93.31 ± 0.13 |  95.73 ± 0.05
AFGRL Lee et al. (2021) | 77.62 ± 0.49 | 89.88 ± 0.33 |  93.22 ± 0.28 | 93.27 ± 0.17 | 95.69 ± 0.10
FastGCL (Ours) | 79.20 ± 0.07 | 89.35 ± 0.09 | 92.91 ± 0.07 | 92.71 ± 0.07 | 95.53 ± 0.02
Table 2: Graph classification performance (Accuracy). Method | NCI1 | PROTEINS | DD | COLLAB | RDT-B | RDT-M5K
---|---|---|---|---|---|---
GIN Xu et al. (2019) | 76.37 ± 1.93 | 73.69 ± 1.79 | 74.53 ± 1.57 | 73.62 ± 0.99 | 87.35 ± 2.61 | 53.29 ± 1.14
GIN R.I.U. | 62.35 ± 1.00 | 63.67 ± 0.67 | 60.00 ± 1.37 | 63.61 ± 1.28 | 63.12 ± 1.34 | 39.23 ± 1.13
InfoGraph Sun et al. (2019) | 69.42 ± 0.58 | 73.03 ± 0.26 | 69.57 ± 0.95 | 68.78 ± 0.47 | 75.58 ± 0.71 | 45.59 ± 0.84
GCC Qiu et al. (2020) | 62.41 ± 0.13 | 72.61 ± 0.23 | 74.09 ± 0.15 |  71.16 ± 0.04 | 77.97 ± 0.24 | 49.31 ± 0.20
GraphCL You et al. (2020) | 70.89 ± 1.25 | 73.02 ± 0.51 | 70.00 ± 1.17 | 69.40 ± 0.72 | 76.11 ± 2.38 | 46.73 ± 0.39
AD-GCL Suresh et al. (2021) | 71.74 ± 0.42 | 73.16 ± 0.51 | 75.00 ± 0.47 | 69.85 ± 0.28 | 78.49 ± 1.04 | 46.36 ± 0.78
FastGCL |  72.68 ± 0.70 |  75.77 ± 0.28 |  76.61 ± 0.44 | 69.34 ± 0.40 |  82.75 ± 2.70 |  51.50 ± 0.91

###  4.1 Experimental Settings

#### Datasets.

To evaluate the performance of our proposed method on different graph tasks,
we conduct experiments on diverse domain datasets. The performance of node
classification is evaluated on 5 medium-scale and widely-used datasets Mernyei
and Cangea (2020); Shchur et al. (2018). The performance of graph
classification is evaluated on 6 datasets with or without node features,
involving molecules, proteins and social networks Morris et al. (2020).

#### Evaluation Protocol.

We train all models in an unsupervised manner. For node classification, we use
a standard linear evaluation protocol Velickovic et al. (2019) to evaluate the
quality of the embeddings learned by the model. Specifically, the parameters
and gradients of the model are frozen and the embeddings are used as input for
classification using a logistic regression with the $\ell_{2}$ regularization.
In addition to using the public 20 divisions on the WikiCS dataset, we
randomly divide the training set, validation set, and test set 20 times
according to 10%, 10%, and 80% on the other datasets Thakoor et al. (2021),
and report the best round on average performance of the validation set as the
evaluation result. For graph classification, Suresh et al. (2021) specifically
states that downstream classifiers can significantly affect the evaluation
results of unsupervised learning. We therefore uniformly use the standard
linear evaluation protocol instead of SVM to evaluate all models and use ten-
fold cross-validation to evaluate the performance of graph classification.

#### Implementation Details.

We use GCN Kipf and Welling (2017) with edge weights as an encoder for the
node-classification task, and the node-wise expression is as follows.

$\mathbf{h}_{v}^{(k)}=\boldsymbol{\Theta}^{\top}\sum_{u\in\mathcal{N}_{v}\cup\\{v\\}}\frac{e_{u,v}}{\sqrt{\hat{d}_{u}\hat{d}_{v}}}\mathbf{h}_{u}^{(k-1)}$
(12)

GIN Xu et al. (2019) with edge weights is applied as an encoder for the graph
classification task. It can be formally denoted as,

$\mathbf{h}_{v}^{(k)}=\textit{MLP}_{\Theta}\left((1+\epsilon)\cdot\mathbf{h}_{v}^{(k-1)}+\sum_{u\in\mathcal{N}_{v}}e_{u,v}\cdot\mathbf{h}_{u}^{(k-1)}\right)$
(13)

The sum pooling is the readout function for graph representation and learn
layer-wise graph representation,

$\mathbf{h}_{\mathcal{G}}=\mathop{\parallel}\limits_{k=1}^{K}\sum_{v\in\mathcal{V}}\mathbf{h}_{v}^{(k)}$
(14)

where $\parallel$ denotes the concatenation operation and $K$ is the number of
layers. For graph data without node features, we simply make a 1-dimensional
unit vector as node features. These settings above are kept consistent with
all baseline models. We implement our proposed method based on the PyTorch and
PyG111https://www.pyg.org/ code libraries and evaluate embeddings more quickly
using a logistic regression classifier on CUDA. The Adam algorithm with weight
decay was used as the optimizer for our model. Hyperparameters such as
learning rate, regular strength $\lambda$, embedding dimension size, and
number of layers of the model were fully examined with grid search.

#### Baselines.

We compare the representative GCL methods mentioned in Section 2. Since the
encoders of the existing methods are basically GCN or GIN, the difference
between their encoders may only be the number of layers, dimension size,
regularization, and implementation. Therefore, we additionally evaluate the
performance of randomly initialized untrained GNNs for node classification
(i.e., GCN, GAT, SGC) and graph classification (i.e., GIN), respectively, to
explore the effect brought by GNN neighbor aggregation. Notably, since our
evaluation using logistic regression on CUDA would be a little worse than the
previously officially reported results for the graph classification task using
SVM, we re-report our experimental results. On the node classification task,
there is not much difference with the previous publicly available results, so
here we directly report them, except for DGI (could be better) and GRACE on
WikiCS (slightly worse). Here we report our reproduced results of them. All
our experiments were run ten times to avoid bias from random initialization
under the same random division.

###  4.2 Experimental Analysis

#### Overall Performance.

The experimental results of node classification and graph classification are
summarized in Table 1 and Table 2, respectively. We can observe from them
that: 1) Our proposed method FastGCL achieves competitive classification
performance on node classification compared to recent GCL method and
outperforms supervised GNNs on some datasets (e.g., WikiCS, Computers). For
the graph classification task, our model exceeds the state-of-the-art methods
on most of datasets. There is a significant improvement on datasets without
node features (i.e., RDT-B, RDT-M5K), which is closer to the performance of
supervised GIN compared to previous methods. 2) Randomly initialized untrained
GNNs also has relatively good performances, even better than DGI, a pioneering
work of GCL, on the node classification task. In comparison with the
evaluation performances of raw feature, the node representation obtained by
neighbor aggregation of GNNs helps the downstream classifier achieve better
performance, especially on Computers. 3) Our method improves on the basis of
randomly initialized untrained GNNs, indicating that our contrastive scheme is
effective. This improvement is obvious on the WikiCS that contains dense node
features. Moreover, on the graph classification datasets without rich node
features, our contrastive scheme drives GNN encoder to learn richer potential
information from the graph structure, obtaining a significant improvement in
classification performance compared to randomly initialized untrained GIN.

#### Training Efficiency.

In Figure 3, we plot the learning curve of all models on the test set of
WikiCS on the node classification task. It can be noticed that FastGCL
converges quickly, with less than 100 epochs. While other models need longer
training epochs for unsupervised learning. In particular, BGRL, which uses the
BYOL framework, maintained its upward trend after 300 epochs, but does not
obtain its best results until training exceeded 9,000 epochs. Within 2000
training epochs, FastGCL outperforms all the comparison methods. The fast
convergence of FastGCL is due to the absence of structure-changed
perturbations, while other models need to implicitly learn a nonlinear mapping
to correct the potential semantic drift brought by data augmentation. Besides,
it benefits from the contrastive learning designed for GNN characteristics,
contrasting weighted and unweighted neighborhood aggregation to quickly
optimize the parameters of the GNN aggregation function.

Figure 3: Learning curve of all node classification baselines on WikiCS.
Figure 4: Performance over training time on WikiCS and RDT-M5K.

Figure 4 shows the training time performance of all models on WikiCS and
RDT-M5K222Since GCC is a general pre-trained graph model that is not specific
to a specific dataset, we have not compared it., respectively, relative to
that of FastGCL. It can be noticed that FastGCL is faster than all previous
methods while achieving competitive classification accuracy. The methods using
complex data augmentation need more time to construct contrastive pairs.
Meanwhile, the method that does not use data augmentation, AF-GCL, is
obviously more time consuming because of the clustering operation used to
identify potential positive sample nodes.

#### Sensitivity Analysis.

We performed sensitivity analysis on several key parameters of FastGCL on
WikiCS, namely the embedding dimension size $d$ and the number of layers $l$
of GNN, as well as the regular strength $\lambda$. The effects of $d$ and $l$
on the randomly initialized untrained GCN and FastGCL are demonstrated in
Figure 5 in (a) and (b), respectively. As can be seen, the FastGCL trend is
mostly consistent with the randomly initialized untrained GCN, and it can be
inferred that the wider dimension and the shallower layers the GNN encoder
has, the better it performs in the downstream task evaluation. The Figure 5(c)
demonstrates the effect of the size of $\lambda$ on the FastGCL learning
stage. Too small $\lambda$ can cause the collapse of learning while too large
$\lambda$ can lead to a slower convergence. A suitable $\lambda$ can make our
contrastive learning more stable.

Figure 5: Effects of key parameters on WikiCS.

#### Ablation Study.

An ablation study has been conducted to explore the effects of our edge-
weighting method with the same hyper-parameters for fair comparison. The
results can find in Table 3. It can be concluded that the edge-weighting
method we designed is effective in guiding the model to discriminate the
importance of the graph topology and highlight critical connections compared
to random weighting.

Table 3: Performance of different edge-weighting methods. Datatset | Random Weighting | FastGCL
---|---|---
WikiCS | 79.11 ± 0.08 | 79.20 ± 0.07
Computers | 89.08 ± 0.09 | 89.35 ± 0.09
NCI1 | 68.18 ± 1.09 | 72.68 ± 0.70
RDT-B | 76.76 ± 3.24 | 82.75 ± 2.70

##  5 Conclusion

In this paper, we propose a simple yet effective GCL method termed FastGCL.
Specifically, it performs contrastive learning by contrasting weighted-
aggregated and non-aggregated neighborhood information, rather than disturbing
the topology and attributes of graph data to generate augmented views. This
aggregation-based contrastive scheme provides a new perspective for GCL
research. Extensive experiments demonstrate that FastGCL has an advantage in
training time and can achieve competitive performance in downstream tasks.

## References

  * Bielak et al. [2021] Piotr Bielak, Tomasz Kajdanowicz, and Nitesh V. Chawla.  Graph barlow twins: A self-supervised representation learning framework for graphs.  arXiv preprint arXiv: Arxiv-2106.02466, 2021.
  * Grill et al. [2020] Jean-Bastien Grill, Florian Strub, Florent Altché, Corentin Tallec, Pierre H. Richemond, Elena Buchatskaya, Carl Doersch, Bernardo Ávila Pires, Zhaohan Guo, Mohammad Gheshlaghi Azar, Bilal Piot, Koray Kavukcuoglu, Rémi Munos, and Michal Valko.  Bootstrap your own latent - A new approach to self-supervised learning.  In NeurIPS, 2020.
  * Gutmann and Hyvärinen [2010] Michael Gutmann and Aapo Hyvärinen.  Noise-contrastive estimation: A new estimation principle for unnormalized statistical models.  In AISTATS, volume 9 of JMLR Proceedings, pages 297–304. JMLR.org, 2010.
  * Kipf and Welling [2017] Thomas N. Kipf and Max Welling.  Semi-supervised classification with graph convolutional networks.  In ICLR (Poster). OpenReview.net, 2017.
  * Lee et al. [2021] Namkyeong Lee, Junseok Lee, and Chanyoung Park.  Augmentation-free self-supervised learning on graphs.  arXiv preprint arXiv: Arxiv-2112.02472, 2021.
  * Linsker [1988] Ralph Linsker.  An application of the principle of maximum information preservation to linear systems.  In NIPS, pages 186–194. Morgan Kaufmann, 1988.
  * Mernyei and Cangea [2020] Péter Mernyei and Cătălina Cangea.  Wiki-cs: A wikipedia-based benchmark for graph neural networks.  arXiv preprint arXiv: Arxiv-2007.02901, 2020.
  * Morris et al. [2020] Christopher Morris, Nils M. Kriege, Franka Bause, Kristian Kersting, Petra Mutzel, and Marion Neumann.  Tudataset: A collection of benchmark datasets for learning with graphs.  arXiv preprint arXiv: Arxiv-2007.08663, 2020.
  * Peng et al. [2020] Zhen Peng, Wenbing Huang, Minnan Luo, Qinghua Zheng, Yu Rong, Tingyang Xu, and Junzhou Huang.  Graph representation learning via graphical mutual information maximization.  In WWW, pages 259–270. ACM / IW3C2, 2020.
  * Qiu et al. [2020] Jiezhong Qiu, Qibin Chen, Yuxiao Dong, Jing Zhang, Hongxia Yang, Ming Ding, Kuansan Wang, and Jie Tang.  GCC: graph contrastive coding for graph neural network pre-training.  In KDD, pages 1150–1160. ACM, 2020.
  * Saunshi et al. [2019] Nikunj Saunshi, Orestis Plevrakis, Sanjeev Arora, Mikhail Khodak, and Hrishikesh Khandeparkar.  A theoretical analysis of contrastive unsupervised representation learning.  In ICML, volume 97 of Proceedings of Machine Learning Research, pages 5628–5637. PMLR, 2019.
  * Shchur et al. [2018] Oleksandr Shchur, Maximilian Mumme, Aleksandar Bojchevski, and Stephan Günnemann.  Pitfalls of graph neural network evaluation.  arXiv preprint arXiv: Arxiv-1811.05868, 2018.
  * Sun et al. [2019] Fan-Yun Sun, Jordan Hoffmann, Vikas Verma, and Jian Tang.  Infograph: Unsupervised and semi-supervised graph-level representation learning via mutual information maximization.  arXiv preprint arXiv: Arxiv-1908.01000, 2019.
  * Sun et al. [2021] Mengying Sun, Jing Xing, Huijun Wang, Bin Chen, and Jiayu Zhou.  Mocl: Contrastive learning on molecular graphs with multi-level domain knowledge.  arXiv preprint arXiv: Arxiv-2106.04509, 2021.
  * Suresh et al. [2021] Susheel Suresh, Pan Li, Cong Hao, and Jennifer Neville.  Adversarial graph augmentation to improve graph contrastive learning.  arXiv preprint arXiv: Arxiv-2106.05819, 2021.
  * Thakoor et al. [2021] Shantanu Thakoor, Corentin Tallec, Mohammad Gheshlaghi Azar, Mehdi Azabou, Eva L. Dyer, Rémi Munos, Petar Veličković, and Michal Valko.  Large-scale representation learning on graphs via bootstrapping.  arXiv preprint arXiv: Arxiv-2102.06514, 2021.
  * Velickovic et al. [2019] Petar Velickovic, William Fedus, William L. Hamilton, Pietro Liò, Yoshua Bengio, and R. Devon Hjelm.  Deep graph infomax.  In ICLR (Poster). OpenReview.net, 2019.
  * Veličković et al. [2017] Petar Veličković, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pietro Liò, and Yoshua Bengio.  Graph attention networks.  arXiv preprint arXiv: Arxiv-1710.10903, 2017.
  * Wu et al. [2019] Felix Wu, Amauri H. Souza Jr., Tianyi Zhang, Christopher Fifty, Tao Yu, and Kilian Q. Weinberger.  Simplifying graph convolutional networks.  In ICML, volume 97 of Proceedings of Machine Learning Research, pages 6861–6871. PMLR, 2019.
  * Wu et al. [2020] Tailin Wu, Hongyu Ren, Pan Li, and Jure Leskovec.  Graph information bottleneck.  In NeurIPS, 2020.
  * Wu et al. [2021] Lirong Wu, Haitao Lin, Zhangyang Gao, Cheng Tan, and Stan. Z. Li.  Self-supervised learning on graphs: Contrastive, generative,or predictive.  arXiv preprint arXiv: Arxiv-2105.07342, 2021.
  * Xu et al. [2019] Keyulu Xu, Weihua Hu, Jure Leskovec, and Stefanie Jegelka.  How powerful are graph neural networks?  In ICLR. OpenReview.net, 2019.
  * You et al. [2020] Yuning You, Tianlong Chen, Yongduo Sui, Ting Chen, Zhangyang Wang, and Yang Shen.  Graph contrastive learning with augmentations.  In NeurIPS, 2020.
  * Zhu et al. [2020] Yanqiao Zhu, Yichen Xu, Feng Yu, Qiang Liu, Shu Wu, and Liang Wang.  Deep graph contrastive representation learning.  arXiv preprint arXiv: Arxiv-2006.04131, 2020.
  * Zhu et al. [2021] Yanqiao Zhu, Yichen Xu, Feng Yu, Qiang Liu, Shu Wu, and Liang Wang.  Graph contrastive learning with adaptive augmentation.  In WWW, pages 2069–2080. ACM / IW3C2, 2021.


11institutetext: Marc Hallin and Dimitri Konen 22institutetext: Université
libre de Bruxelles, Belgium, and University of Warwick, United Kingdom
22email<EMAIL_ADDRESS><EMAIL_ADDRESS>
# Multivariate Quantiles: Geometric and Measure-Transportation-Based Contours

Marc Hallin and Dimitri Konen

###### Abstract

Quantiles are a fundamental concept in probability and theoretical statistics
and a daily tool in their applications. While the univariate concept of
quantiles is quite clear and well understood, its multivariate extension is
more problematic. After half a century of continued efforts and many
proposals, two concepts, essentially, are emerging: the so-called (relabeled)
geometric quantiles, extending the characterization of univariate quantiles as
minimizers of an L1 loss function involving the check functions, and the more
recent center-outward quantiles based on measure transportation ideas. These
two concepts yield distinct families of quantile regions and quantile
contours. Our objective here is to present a comparison of their main
theoretical properties and a numerical investigation of their differences.

##  1 Introduction

Quantiles are a fundamental concept in probability and theoretical statistics
and a familiar tool in their applications, ranging from descriptive statistics
and data analysis to statistical inference. The univariate concept of a
quantile function (the inverse of a distribution function) is well understood
and well studied in dimension $d=1$. Despite half a century of efforts,
however, the multivariate extension (dimension $d\geq 2$) of this essential
notion remains quite problematic—the obvious difficulty being the absence of a
canonical ordering of the real space for dimension two and higher.

Various proposals have been made in the literature, none of which is fully
agreed upon. Two concepts, essentially, are emerging from these many attempts:
the so-called (relabeled) geometric quantiles (Chaudhuri 1996), extending the
characterization of univariate quantiles as minimizers of an L1 loss function
involving the check functions, and the more recent center-outward quantiles
(Hallin et al. 2021) based on measure transportation ideas.

The objective of this note is a brief comparison of some of the properties of
the quantile regions and quantile contours associated with these two concepts.

##  2 Quantiles: from univariate to multivariate

Throughout, the inner product of two vectors $x,y\in{\mathbb{R}}^{d}$ is
denoted as $x^{\prime}y$; $\mathbb{B}_{d}$ and $\mathbb{S}^{d-1}$ stand for
the open unit ball and the $(d-1)$-dimensional unit sphere, respectively, in
${\mathbb{R}}^{d}$.

###  2.1 Definitions and basic properties

Let $X$ denote a Lebesgue-absolutely continuous real-valued random variable
with distribution $\rm P$ over ${\mathbb{R}}$ and distribution function

$F_{\rm P}:x\mapsto F_{\rm P}(x)\coloneqq{\rm P}(X\leq x).$

To simplify the exposition, assume that $X$ admits a nonvanishing density
$f\coloneqq{\rm dP/d\lambda}$ with respect to the Lebesgue measure $\lambda$
so that, in particular, $F_{\rm P}$ is continuous and strictly monotone
increasing, hence invertible. The quantile function of $\rm P$ (equivalently,
the quantile function of $X\sim\rm P$) then is defined as the inverse $Q_{\rm
P}\coloneqq F_{\rm P}^{-1}$ of the distribution function $F_{\rm P}$ and the
quantile of order $\alpha\in(0,1)$ of $\rm P$ (of $X\sim{\rm P}$) is the value
$Q_{\rm P}(\alpha)$ of $Q_{\rm P}$ at $\alpha$.

Along with quantiles of given order $\alpha\in(0,1)$, it is often convenient
to consider central quantile regions and contours. In dimension $d=1$, call
quantile region of order $\tau\in~{}\\![0,1)$ of $\rm P$ (of $X\sim{\rm P}$)
the closed interval

$\mathscr{R}_{\rm P}(\tau)\coloneqq\Big{[}Q_{\rm
P}\Big{(}\frac{1-\tau}{2}\Big{)},\,Q_{\rm
P}\Big{(}\frac{1+\tau}{2}\Big{)}\Big{]}$

(an interquantile interval) and quantile contour of order $\tau$ of $\rm P$
(of $X\sim{\rm P}$) the boundary

${\mathscr{C}_{\rm P}(\tau)}\coloneqq\Big{\\{}Q_{\rm
P}\Big{(}\frac{1-\tau}{2}\Big{)},\,Q_{\rm
P}\Big{(}\frac{1+\tau}{2}\Big{)}\Big{\\}}$

of $\mathscr{R}_{\rm P}(\tau)$.

The quantile region of order 0 is the degenerate interval $\mathscr{R}_{\rm
P}({0})=\\{Q_{\rm P}(1/2)\\}$ and has Lebesgue measure zero: call it the
median region. Similarly, the quantile region of order 1/2 is the traditional
interquartile interval: $\mathscr{R}_{\rm P}({1/2})=[Q_{\rm P}(1/4),Q_{\rm
P}(3/4)]$. Contrary to the quantile halflines $(-\infty,Q_{\rm P}(\alpha)]$,
which very much depend on the orientation of the real line (an orientation
that is no longer meaningful in ${\mathbb{R}}^{d}$ for $d\geq~{}\\!2$), the
concepts of quantile regions and contours of order $\tau$ are “center-outward”
concepts (the “center” being the median $Q_{\rm P}(1/2)$) and are invariant
under a change of the orientation of ${\mathbb{R}}$ (an orthogonal
transformation). The collection $\left\\{\mathscr{R}_{\rm
P}(\tau):\,\tau\in[0,1)\right\\}$ of quantile regions is stricly nested (in
the sense of set inclusion) as $\tau$ ranges from zero to one, with
$\mathscr{R}_{\rm P}({0})=\bigcap_{\tau\in[0,1)}{\mathscr{R}_{\rm P}}(\tau)$.
It is easily seen that ${\mathscr{C}_{\rm P}(|2F_{\rm P}(x)-1|)}$ is the
unique quantile contour running through $x\in\mathbb{R}$; actually,

${\mathscr{R}_{\rm P}(\tau)}=\left\\{x:|2F_{\rm P}(x)-1|\leq\tau\right\\}\
\text{ and }\ {\mathscr{C}_{\rm P}(\tau)}=\left\\{x:|2F_{\rm
P}(x)-1|=\tau\right\\}\quad\tau\in[0,1).$ (1)

An essential property of a quantile function is that the $\rm P$-probability
content of a quantile region $\mathscr{R}_{\rm P}(\tau)$ is $\tau$
irrespective of $\rm P$: namely, ${\rm P}(\mathscr{R}_{\rm P}(\tau))=\tau$ for
all $\rm P$. What would be the relevance, indeed, of an interquartile region
${\mathscr{R}_{\rm P}(1/2)}$ with probability content 0.5 under ${\rm
P=P}_{1}$ (as should be) but probability content 0.6 under ${\rm P=P}_{2}$,
and probability content 0.7 under ${\rm P=P}_{3}$? This property, in dimension
$d=1$, is equivalent to the property of $F_{\rm P}(X)$ being uniform over
$[0,1)$ for $X\sim{\rm P}$.

When the center-outward quantile regions $\mathscr{R}_{\rm P}(\tau)$ and
contours $\mathscr{C}_{\rm P}(\tau)$ are the main points of interest, it is
natural to replace $F_{\rm P}$ and $Q_{\rm P}$ by center-outward counterparts:
call $F^{\scriptscriptstyle\pm}_{\rm P}\coloneqq 2F_{\rm P}-1$ and
$Q^{\scriptscriptstyle\pm}_{\rm P}\coloneqq(F^{\scriptscriptstyle\pm}_{\rm
P})^{-1}$ the center-outward distribution and quantile functions of $\rm P$
(or $X\sim{\rm P}$), respectively. Just like $F_{\rm P}$ and $Q_{\rm P}$, the
maps $F^{\scriptscriptstyle\pm}_{\rm P}$ and $Q^{\scriptscriptstyle\pm}_{\rm
P}$ are monotone increasing and carry the same information about $\rm P$,
which they fully characterize. Instead of being uniformly distributed over
$[0,1)$, $F^{\scriptscriptstyle\pm}_{\rm P}(X)$, when $X\sim{\rm P}$, is
uniform over $(-1,1)$, which is the open unit ball $\mathbb{B}_{1}$ in
${\mathbb{R}}$. In fact, denoting by ${\rm U}_{1}$ the uniform distribution
over $\mathbb{B}_{1}$, we have $F^{\scriptscriptstyle\pm}_{\rm P}(X)\sim
U_{1}$ when $X\sim{\rm P}$, hence $Q^{\scriptscriptstyle\pm}_{\rm
P}(U)\sim{\rm P}$ if $U\sim{\rm U}_{1}$. In the terminology of measure
transportation, we say that the transport $F^{\scriptscriptstyle\pm}_{\rm P}$
is pushing $\rm P$ forward to ${\rm U}_{1}$ and that its inverse
$Q^{\scriptscriptstyle\pm}_{\rm P}$ is pushing ${\rm U}_{1}$ forward to ${\rm
P}$, which we conveniently denote as

$F^{\scriptscriptstyle\pm}_{\rm P}\\#{\rm P}={\rm U}_{1}\quad\text{and }\quad
Q^{\scriptscriptstyle\pm}_{\rm P}\\#{\rm U}_{1}={\rm P},$

respectively. The quantile region $\mathscr{R}_{\rm P}(\tau)$ and quantile
contour $\mathscr{C}_{\rm P}(\tau)$ then take the simple forms (equivalent to
(1))

$\mathscr{R}_{\rm P}(\tau)=Q^{\scriptscriptstyle\pm}_{\rm
P}(\tau\mathbb{B}_{1})\quad\text{and}\quad\mathscr{C}_{\rm
P}(\tau)=Q^{\scriptscriptstyle\pm}_{\rm
P}(\tau\mathbb{S}^{0}),\quad\tau\in[0,1),$ (2)

respectively, where $\mathbb{S}^{0}\coloneqq\\{-1,1\\}$ stands for the
($0$-dimensional) unit sphere in ${\mathbb{R}}$.

###  2.2 Characterization of univariate quantiles as minimizers of expected
check functions

When $\rm P$ is a probability measure over ${\mathbb{R}}$ and $X\sim{\rm P}$,
and given a level $\alpha\in[0,1)$, a quantile of order $\alpha$ of $X$, by
definition, is any $x\in{\mathbb{R}}$ satisfying

${\rm P}(X<x)\leq\alpha\leq{\rm P}(X\leq x).$

It turns out that such $x$’s are exactly the minimizers (with respect to
$\theta\in{\mathbb{R}}$) of the objective function

$\theta\mapsto M_{\alpha}^{\rm P}(\theta):={\rm E}_{\rm
P}[\rho_{\alpha}(\theta-X)-\rho_{\alpha}(X)]$ (3)

where, denoting by $z^{+}$ and $z^{-}$ the positive and negative parts,
respectively, of $z\in\mathbb{R}$,

$\rho_{\alpha}:z\mapsto\rho_{\alpha}(z)\coloneqq|z|+(2\alpha-1)z=2\big{\\{}\alpha
z^{+}+(1-\alpha)z^{-}\big{\\}},\quad z\in{\mathbb{R}},\ \alpha\in[0,1)$

stands for the so-called check function. Noting that $M_{\alpha}^{\rm
P}(\theta)$ rewrites as

$M_{\alpha}^{\rm P}(\theta)={\rm E}_{\rm
P}\big{[}|\theta-X|-|X|\big{]}-(2\alpha-1)\theta$ (4)

for all $\theta\in{\mathbb{R}}$, it follows from the triangular inequality
that $\big{|}|X-\theta|-|X|\big{|}\leq|\theta|$, which is integrable. The
objective function $M_{\alpha}^{\rm P}$, hence the correspondence between
quantiles and the minimizing procedure in (3), thus remains well defined for
an arbitrary probability measure ${\rm P}$, without any moment assumption. In
this general setting, we have, even in the absence of a finite moment of order
one,

$Q_{\rm P}(\alpha)={\rm argmin}_{\theta\in{\mathbb{R}}}M_{\alpha}^{\rm
P}(\theta),\qquad\alpha\in[0,1).$

For $\alpha=0.5$, $M^{\rm P}_{\alpha}(\theta)={\rm E}_{\rm
P}[|\theta-X|-|X|]$; in case ${\rm E_{P}}(X)<\infty$, thus, we recover the
familiar characterization of the median $Q_{\rm P}(1/2)$ as an ${\rm L}_{1}$
location parameter, the minimizer of the expected absolute deviation ${\rm
E}_{\rm P}[|\theta-X|]$.

Recalling that $F_{\rm P}^{\scriptscriptstyle\pm}\coloneqq 2F_{\rm P}-1$, we
see that the natural labeling for center-outward quantiles is provided by
$\beta:=2\alpha-1\in(-1,1)$. Therefore, the center-outward version of the
characterization of quantiles as minimizers is obtained by relabeling
$2\alpha-1$ as $\beta$ in (3). This leads to the characterization of center-
outward quantiles $Q_{\rm P}^{\scriptscriptstyle\pm}(\beta)$ as minimizers
with respect to $\theta\in{\mathbb{R}}$ of the objective function

$\theta\mapsto O_{\beta}^{\rm P}(\theta)\coloneqq{\rm E}_{\rm
P}\big{[}|\theta-X|-|X|\big{]}-\beta\theta,\qquad\beta\in(-1,1).$ (5)

We then have (still, without any moment assumption)

$Q_{\rm P}^{\scriptscriptstyle\pm}(\beta)={\rm
argmin}_{\theta\in{\mathbb{R}}}O_{\beta}^{\rm
P}(\theta),\qquad\beta\in(-1,1).$

The value $\beta=0$ yields the median $Q^{\scriptscriptstyle\pm}_{\rm
P}(0)=Q_{\rm P}(1/2)$ of $\rm P$, while positive (resp. negative) values of
$\beta$ are associated to quantiles sitting to the right-hand (resp. left-
hand) side of the median.

###  2.3 Characterization of univariate quantile functions as monotone
transports

Consider a probability measure ${\rm P}$ with non-vanishing density $f$ with
respect to the Lebesgue measure on ${\mathbb{R}}$. An alternative to the
characterization, in dimension $d=1$, of $Q^{\scriptscriptstyle\pm}_{\rm P}$
(hence of $Q_{\rm P}$ and $F^{\scriptscriptstyle\pm}_{\rm P}$) as minimizers
of expected check functions is obtained by noting that
$Q^{\scriptscriptstyle\pm}_{\rm P}$ is monotone increasing and pushes ${\rm
U}_{1}$ forward to $\rm P$. This fully characterizes
$Q^{\scriptscriptstyle\pm}_{\rm P}$ up to a set of values with Lebesgue
measure zero.111As we shall see, the existence and (almost everywhere)
uniqueness, in dimension $d=1$, of a monotone non-decreasing mapping (that is,
the derivative, hence the gradient, of a convex function) pushing ${\rm
U}_{1}$ forward to $\rm P$ is a particular case of a famous and more general
measure transportation result by McCann (1995). The center-outward quantile
function $Q^{\scriptscriptstyle\pm}_{\rm P}$, thus, can be characterized as
the almost everywhere unique monotone map from the open unit ball
$\mathbb{B}_{1}=(-1,1)$ to ${\mathbb{R}}$ such that
$Q^{\scriptscriptstyle\pm}_{\rm P}\\#{\rm U}_{1}={\rm P}$.

##  3 Multivariate quantile functions, regions, and contours

Extending to dimension $d>1$ the characterizations of quantile functions
developed in Section 2.2 leads to the concept of geometric quantiles first
proposed by Chaudhuri (1996) while the measure-transportation-based
characterization developed in Section 2.3 leads to the concept of multivariate
center-outward quantiles proposed in Hallin et al. (2021).

###  3.1 Geometric quantiles and geometric quantile contours

Consider a probability measure ${\rm P}$ supported on ${\mathbb{R}}^{d}$, with
$d\geq 2$. Extending toa $d$-dimensional framework the characterization of
univariate quantiles as minimizers of expected check functions is obtained by
replacing, in (5), the absolute values with Euclidean norms and
$\beta\in(-1,1)=\mathbb{B}_{1}$ with ${\bf v}\in\mathbb{B}_{d}$. Letting ${\bf
v}\eqqcolon\tau{\bf u}\in\mathbb{B}_{d}$with $\tau=\|{\bf v}\|\in(0,1)$ and
${\bf u}={\bf v}/\|{\bf v}\|\in\mathbb{S}^{d-1}$, define the geometric
quantile of order $\tau>0$ in direction ${\bf u}={\bf v}/\|{\bf v}\|$ of $\rm
P$ (of ${\bf X}\sim{\rm P}$) as an arbitrary minimizer, over ${\bf
z}\in{\mathbb{R}}^{d}$, of the objective function

${\bf z}\mapsto O_{\tau,{\bf u}}^{\rm P}({\bf z}):={\rm E}_{\rm
P}\big{[}\|{\bf z}-{\bf X}\|-\|{\bf X}\|\big{]}-\tau{\bf u}^{\prime}{\bf z}$
(6)

where ${\bf X}\sim{\rm P}$; throughout, boldface is used to stress that the
variables take values in ${\mathbb{R}}^{d}$. The objective function
$O_{\tau,{\bf u}}^{\rm P}$ being convex, it can be shown that the geometric
quantiles of ${\rm P}$ are unique as soon as $\rm P$ is not supported on a
single line of ${\mathbb{R}}^{d}$—a situation which essentially boils down to
a univariate setting: see Theorem 1 in Paindaveine and Virta (2020) for a
discussion. For $\tau=0$, define the geometric median as the minimizer of
${\rm E}_{\rm P}\big{[}\|{\bf z}-{\bf X}\|-\|{\bf X}\|\big{]}$ which, still in
view of the triangular inequality, exists without any moment assumption and is
better known as the Fréchet median.

When uniqueness holds ($\rm P$ not supported on a line), we denote by ${\bf
Q}^{\scriptscriptstyle{\text{\rm g}}}_{\rm P}(\tau{\bf u})$ the geometric
quantile of order $\tau\in[0,1)$ in direction $\bf u$ of ${\rm P}$ (of ${\bf
X}\sim{\rm P}$) and by ${\bf Q}^{\scriptscriptstyle{\text{\rm g}}}_{\rm
P}({\boldsymbol{0}})$ the geometric or Fréchet median, which can be
interpreted as the geometric quantile of order $\tau=0$ in any direction $\bf
u$. Here again, boldface is used to stress that ${\bf
Q}^{\scriptscriptstyle\pm}_{\rm P}$ takes values in ${\mathbb{R}}^{d}$.

Strongly related to the gradient of $O_{\tau,{\bf u}}^{\rm P}$, we define the
geometric distribution function ${\bf F}^{\scriptscriptstyle{\text{\rm
g}}}_{\rm P}$ of ${\rm P}$ as the mapping

${\bf z}\mapsto{\bf F}^{\scriptscriptstyle{\text{\rm g}}}_{\rm P}({\bf
z})\coloneqq{\rm E}_{\rm P}\Big{[}\frac{{\bf z}-{\bf X}}{\|{\bf z}-{\bf
X}\|}{\mathbb{1}}[{\bf X}\neq{\bf z}]\Big{]},\qquad{\bf
z}\in{\mathbb{R}}^{d},$ (7)

where ${\mathbb{1}}[A]$ stands for the indicator function of $A$. While the
definition of ${\bf Q}^{\scriptscriptstyle{\text{\rm g}}}_{\rm P}$ is
motivated by the characterization of univariate center-outward quantiles as
minimizers of an objective function, let us stress that the geometric
distribution function in (7) is a natural analogue of the univariate center-
outward distribution as well. Indeed, recalling that, for $d=1$, $F_{\rm
P}^{\scriptscriptstyle\pm}\coloneqq 2F_{\rm P}-1$ and observing that
$2{\mathbb{1}}[X\leq z]-1={\rm sign}(z-X)$, one obtains

$F_{\rm P}^{\scriptscriptstyle\pm}(z)={\rm E}_{\rm P}[{\rm sign}(z-X)],\quad
z\in{\mathbb{R}},$ (8)

to which ${\bf F}^{\scriptscriptstyle{\text{\rm g}}}_{\rm P}$ in (7) reduces
when $d=1$. The direction $({\bf z}-{\bf X})/\|{\bf z}-{\bf X}\|$, thus, can
be interpreted as a multivariate sign for ${\bf z}-{\bf
X}\neq{\boldsymbol{0}}$.

Provided that ${\rm P}$ is non-atomic and is not supported on a single line of
${\mathbb{R}}^{d}$, Theorem 6.2 in Konen and Paindaveine (2022) entails that
${\bf F}^{\scriptscriptstyle{\text{\rm g}}}_{\rm P}$ is a homeomorphism
between ${\mathbb{R}}^{d}$ and the open unit ball $\mathbb{B}_{d}$ of
${\mathbb{R}}^{d}$, with inverse ${\bf Q}^{\scriptscriptstyle{\text{\rm
g}}}_{\rm P}$. Consequently, one also can define ${\bf
Q}^{\scriptscriptstyle{\text{\rm g}}}_{\rm P}$ as $\left({\bf
F}^{\scriptscriptstyle{\text{\rm g}}}_{\rm P}\right)^{-1}$ with ${\bf
F}^{\scriptscriptstyle{\text{\rm g}}}_{\rm P}$ defined in (7), or define ${\bf
F}^{\scriptscriptstyle{\text{\rm g}}}_{\rm P}$ as $\left({\bf
Q}^{\scriptscriptstyle{\text{\rm g}}}_{\rm P}\right)^{-1}$ with ${\bf
Q}^{\scriptscriptstyle{\text{\rm g}}}_{\rm P}$ defined as the minimizer of
$O_{\tau,{\bf u}}^{\rm P}$ in (6). Further results about the regularity of
${\bf F}^{\scriptscriptstyle{\text{\rm g}}}_{\rm P}$ and ${\bf
Q}^{\scriptscriptstyle{\text{\rm g}}}_{\rm P}$ are provided in Konen (2022).

With this, one can define geometric quantile regions and contours, which are
multivariate extensions of the univariate center-outward quantile regions and
contours, by letting

$\boldsymbol{\mathscr{R}}_{{\rm P}}^{\scriptscriptstyle{\text{\rm
g}}}(\tau):={\bf Q}^{\scriptscriptstyle{\text{\rm g}}}_{\rm
P}(\tau\mathbb{B}_{d})\quad\textrm{and}\quad\boldsymbol{\mathscr{C}}_{\rm
P}^{\scriptscriptstyle{\text{\rm g}}}(\tau):={\bf
Q}^{\scriptscriptstyle{\text{\rm g}}}_{\rm
P}(\tau\mathbb{S}^{d-1}),\quad\tau\in[0,1).$ (9)

For $\tau=0$, $\boldsymbol{\mathscr{R}}_{\rm P}^{\scriptscriptstyle{\text{\rm
g}}}(0)$ and $\boldsymbol{\mathscr{C}}_{\rm P}^{\scriptscriptstyle{\text{\rm
g}}}(0)$ reduce to the singlepoint set consisting of the geometric or Fréchet
median; as in the univariate case, the geometric median is an ${\rm L}_{1}$
location parameter that minimizes, when finite first-order moments exist, the
expected absolute deviation (absolute deviations, here, are to be understood
as Euclidean norms of deviations).

###  3.2 Measure-transportation-based quantiles and quantile contours

Denote by ${\rm U}_{d}$ the spherical uniform over the open unit ball
$\mathbb{B}_{d}$, that is, the product of a uniform distribution over the unit
sphere $\mathbb{S}^{d-1}$ and a (univariate) uniform over the distances from
the origin; in other words, ${\rm U}_{d}$ is the probability measure with
density $f_{{\rm U}_{d}}({\bf x})$ proportional to $1/\|{\bf
x}\|^{d-1}{\mathbb{1}}[{\boldsymbol{0}}<\|{\bf x}\|\leq 1]$ for
${\boldsymbol{0}}<\|{\bf x}\|$ and $f_{{\rm U}_{d}}({\boldsymbol{0}})=\infty$.

A celebrated result by McCann (1995) implies that, for any distribution $\rm
P$ over ${{\mathbb{R}}}^{d}$, there exists an almost everywhere unique
gradient of a convex function pushing $\rm P$ forward to the spherical uniform
${\rm U}_{d}$ over the (open) unit ball $\mathbb{B}_{d}$ of
${\mathbb{R}}^{d}$. In dimension $d=1$, a gradient of convex function is a
monotone non-decreasing function and the unique monotone non-decreasing
function pushing $\rm P$ forward to ${\rm U}_{1}$ is the center-outward
distribution function $F^{\scriptscriptstyle\pm}_{\rm P}\coloneqq 2F_{\rm
P}-1$. For arbitrary dimension $d\geq~{}\\!1$, this leads Hallin et al. (2021)
to define the center-outward distribution function of $\rm P$ (or ${\bf
X}\sim{\rm P}$) as the almost everywhere unique gradient of a convex function
${\bf F}^{\scriptscriptstyle\pm}_{\rm P}$ pushing $\rm P$ to the spherical
uniform ${\rm U}_{d}$.

Assuming that $\rm P$ belongs to the class $\mathcal{P}_{d}^{+}$ of Lebesgue-
absolutely continuous distributions with density bounded away from zero and
infinity on any compact222Formally, for any compact subset
$K\subset{\mathbb{R}}^{d}$, there exist constants $0<c^{-}_{K}\leq
c_{K}^{+}<\infty$ such that $c_{K}^{-}\leq{\rm dP}/{\rm d\lambda}({\bf x})\leq
c_{K}^{+}$ for all ${\bf x}\in K$. This assumption is maintained here for
clarity of exposition, but it can be relaxed: see del Barrio et al. (2020) and
del Barrio and Gonzáles-Sanz (2023). subset of ${\mathbb{R}}^{d}$, Figalli
(2018) showed that ${\bf F}^{\scriptscriptstyle\pm}_{\rm P}$ is a
homeomorphism between $\mathbb{R}^{d}\\!\setminus\\!({\bf
F}^{\scriptscriptstyle\pm}_{\rm P})^{-1}({\bf 0})$ and the punctured unit ball
$\mathbb{B}_{d}\\!\setminus\\!\\{{\bf 0}\\}$. The restriction of ${\bf
F}^{\scriptscriptstyle\pm}_{\rm P}$ to $\mathbb{R}^{d}\\!\setminus\\!({\bf
F}^{\scriptscriptstyle\pm}_{\rm P})^{-1}({\bf 0})$ is thus continuously
invertible. We call its (continuous) inverse ${\bf
Q}^{\scriptscriptstyle\pm}_{\rm P}\coloneqq({\bf
F}^{\scriptscriptstyle\pm}_{\rm P})^{-1}$ the center-outward quantile function
of $\rm P$ (or ${\bf X}\sim{\rm P}$). This yields center-outward quantile
regions and contours, of order $\tau$ which are multivariate extensions of the
univariate center-outward quantile regions and contours defined in (2), of the
form

$\boldsymbol{\mathscr{R}}_{\rm P}^{\scriptscriptstyle\pm}(\tau)\coloneqq{\bf
Q}^{\scriptscriptstyle\pm}_{\rm
P}(\tau\mathbb{B}_{d})\quad\text{and}\quad\boldsymbol{\mathscr{C}}_{\rm
P}^{\scriptscriptstyle\pm}(\tau)\coloneqq{\bf Q}^{\scriptscriptstyle\pm}_{\rm
P}(\tau\mathbb{S}^{d-1}),\qquad\tau\in(0,1),$ (10)

respectively. As for $\tau=0$, define $\boldsymbol{\mathscr{R}}_{\rm
P}^{\scriptscriptstyle\pm}(0)\coloneqq\bigcap_{0<\tau<1}\boldsymbol{\mathscr{R}}_{\rm
P}^{\scriptscriptstyle\pm}(\tau)$ as the center-outward median region, which
is shown to be convex and compact, with Hausdorff dimension at most $d-1$
(Figalli, 2018). Note the parallel between (10) and (9).

##  4 Main properties

###  4.1 Geometric quantile regions and contours

Consider a probability measure ${\rm P}$ which is non-atomic and is not
supported on a single line of ${\mathbb{R}}^{d}$; in particular, ${\bf
F}^{\scriptscriptstyle{\text{\rm g}}}_{\rm P}$ and ${\bf
Q}^{\scriptscriptstyle{\text{\rm g}}}_{\rm P}$ are homeomorphisms with ${\bf
Q}^{\scriptscriptstyle{\text{\rm g}}}_{\rm P}=({\bf
F}^{\scriptscriptstyle{\text{\rm g}}}_{\rm P})^{-1}$. Then, geometric quantile
regions and contours have the following topological properties:

  1. (i)

the regions $\boldsymbol{\mathscr{R}}_{\rm P}^{\scriptscriptstyle{\text{\rm
g}}}(\tau)$ are compact and arc-connected for all $\tau\in[0,1)$;

  2. (ii)

the regions $\boldsymbol{\mathscr{R}}_{\rm P}^{\scriptscriptstyle{\text{\rm
g}}}(\tau)$ are stricly nested as $\tau$ increases: for any
$0\leq\tau_{1}<\tau_{2}<1$, we have $\boldsymbol{\mathscr{R}}_{\rm
P}^{\scriptscriptstyle{\text{\rm
g}}}(\tau_{1})\subsetneq\boldsymbol{\mathscr{R}}_{\rm
P}^{\scriptscriptstyle{\text{\rm g}}}(\tau_{2})$;

  3. (iii)

if ${\rm P}$ has a bounded density $f$, then the contour
$\boldsymbol{\mathscr{C}}_{\rm P}^{\scriptscriptstyle{\text{\rm g}}}(\tau)$ is
a $(d-1)$-dimensional manifold of class $\mathcal{C}^{d-1}$ for all
$\tau\in(0,1)$; if, in addition, $f$ is of class $\mathcal{C}^{k,\alpha}$ for
some $k\in\mathbb{N}$ and $\alpha\in(0,1)$, then
$\boldsymbol{\mathscr{C}}_{\rm P}^{\scriptscriptstyle{\text{\rm g}}}(\tau)$ is
of class $\mathcal{C}^{d+k}$ for all $\tau\in(0,1)$;

  4. (iv)

for $d=1$, ${\bf F}^{\scriptscriptstyle{\text{\rm g}}}_{\rm P}$, ${\bf
Q}^{\scriptscriptstyle{\text{\rm g}}}_{\rm P}$, $\boldsymbol{\mathscr{R}}_{\rm
P}^{\scriptscriptstyle{\text{\rm g}}}$, and $\boldsymbol{\mathscr{C}}_{\rm
P}^{\scriptscriptstyle{\text{\rm g}}}$ reduce to
$F^{\scriptscriptstyle\pm}_{\rm P}$, $Q^{\scriptscriptstyle\pm}_{\rm P}$,
${\mathscr{R}}_{\rm P}^{\scriptscriptstyle\pm}$, and ${\mathscr{C}}_{\rm
P}^{\scriptscriptstyle\pm}$.

While properties (i)-(ii) are straightforward consequences of the continuity
of ${\bf Q}^{\scriptscriptstyle{\text{\rm g}}}_{\rm P}$, property (iii)
requires more advanced mathematical tools. That property was established in
Konen (2022), where it is shown that geometric distribution functions satisfy
a linear partial differential equation, which was the key to te derivation of
the regularity properties of geometric contours. Among other things, this PDE
implies that geometric distribution functions fully characterize arbitrary
probability measures (no moment or non-atomicity assumptions required),
namely,

  1. (v)

if ${\rm P_{1}}$ and ${\rm P_{2}}$ are arbitrary probability measures such
that ${\bf F}_{\rm P_{1}}^{\scriptscriptstyle{\text{\rm g}}}={\bf F}_{\rm
P_{2}}^{\scriptscriptstyle{\text{\rm g}}}$, then ${\rm P}_{1}={\rm P}_{2}$.

This characterization result has been long known, though: the first proof was
provided in Theorem 2.5 of Koltchinski (1997) and was refined recently in
Konen (2022).

Geometric quantile and distribution functions further enjoy the following
natural equivariance properties; see Girard and Stupfler (2017) and Konen and
Paindaveine (2022).

  1. (vi)

Let $\boldsymbol{\theta}\in{\mathbb{R}}^{d}$ and denote by ${\bf O}$ an
arbitrary $d\times d$ orthogonal matrix. For ${\bf X}\sim{\rm P}$, write ${\bf
F}_{{\bf X}}^{\scriptscriptstyle{\text{\rm g}}},{\bf F}_{{\bf
X}}^{\scriptscriptstyle{\text{\rm g}}}$, etc. instead of ${\bf F}_{\rm
P}^{\scriptscriptstyle{\text{\rm g}}},{\bf F}_{\rm
P}^{\scriptscriptstyle{\text{\rm g}}}$, etc. We have

$\displaystyle{\bf F}_{\boldsymbol{\theta}+{\bf O}{\bf
X}}^{\scriptscriptstyle{\text{\rm g}}}(\boldsymbol{\theta}+{\bf O}{\bf
x})={\bf O}{\bf F}_{\bf X}^{\scriptscriptstyle{\text{\rm g}}}({\bf x})\quad$
$\displaystyle{\bf x}\in\mathbb{R}^{d},$ $\displaystyle{\bf
Q}_{\boldsymbol{\theta}+{\bf O}{\bf X}}^{\scriptscriptstyle{\text{\rm
g}}}({\bf O}{\bf u})=\boldsymbol{\theta}+{\bf O}{\bf Q}_{\bf
X}^{\scriptscriptstyle{\text{\rm g}}}({\bf u})\quad$ $\displaystyle{\bf
u}\in\mathbb{B}_{d},$
$\displaystyle\boldsymbol{\mathscr{R}}_{\boldsymbol{\theta}+{\bf O}{\bf
X}}^{\scriptscriptstyle{\text{\rm g}}}(\tau)=\boldsymbol{\theta}+{\bf
O}\boldsymbol{\mathscr{R}}_{{\bf X}}^{\scriptscriptstyle{\text{\rm
g}}}(\tau)\quad$ $\displaystyle\tau\in[0,1),$ and
$\displaystyle\boldsymbol{\mathscr{C}}_{\boldsymbol{\theta}+{\bf O}{\bf
X}}^{\scriptscriptstyle{\text{\rm g}}}(\tau)=\boldsymbol{\theta}+{\bf
O}\boldsymbol{\mathscr{C}}_{\bf X}^{\scriptscriptstyle{\text{\rm
g}}}(\tau)\quad$ $\displaystyle\tau\in[0,1).$

In particular, if $\rm P$ enjoys symmetry properties (rotational, central, or
with respect to some hyperplanes), ${\bf F}^{\scriptscriptstyle{\text{\rm
g}}}_{\rm P}$, ${\bf Q}^{\scriptscriptstyle{\text{\rm g}}}_{\rm P}$,
$\boldsymbol{\mathscr{R}}_{\rm P}^{\scriptscriptstyle{\text{\rm g}}}(\tau)$,
and $\boldsymbol{\mathscr{C}}_{\rm P}^{\scriptscriptstyle{\text{\rm
g}}}(\tau)$ enjoy the same symmetry properties as $\rm P$ for all
$\tau\in~{}\\![0,1)$.

Recall that an essential property of distribution functions in dimension $d=1$
is their non-decreasing monotonicity. In higher dimensions, the concept of
monotonicity still makes sense; we distinguish two notions of monotonicity,
which both reduce to the usual concept of (non-decreasing) monotonicity when
$d=1$: monotonicity and cyclical monotonicity. We say that a map ${\bf
G}:{\mathbb{R}}^{d}\to{\mathbb{R}}^{d}$ is monotone if

$({\bf G}({\bf y})-{\bf G}({\bf x}))^{\prime}({\bf y}-{\bf x})\geq 0$

for all ${\bf x},{\bf y}\in{\mathbb{R}}^{d}$. We say that $G$ is cyclically
monotone if for any $m\geq 1$ we have

$\sum_{k=1}^{m}\big{(}{\bf G}({\bf x}_{k+1})-{\bf G}({\bf
x}_{k})\big{)}^{\prime}{\bf x}_{k+1}\geq 0$

for every cycle $x_{1},x_{2},\ldots,x_{m+1}:=x_{1}$. For $m=2$, cyclical
monotonicity yields

$\big{(}{\bf G}({\bf x}_{2})-{\bf G}({\bf x}_{1})\big{)}^{\prime}{\bf
x}_{2}+\big{(}{\bf G}({\bf x}_{1})-{\bf G}({\bf x}_{2})\big{)}^{\prime}{\bf
x}_{1}\geq 0,$

and cyclical monotonicity thus implies monotonicity. Theorem 24.8 in
Rockafellar (1970) essentially entails that ${\bf G}$ is cyclically monotone
precisely when it writes as the gradient ${\bf G}=\nabla\psi$ of some convex
(potential) function $\psi:{\mathbb{R}}^{d}\to{\mathbb{R}}$.

The following then holds.

  1. (vii)

The maps ${\bf F}^{\scriptscriptstyle{\text{\rm g}}}_{\rm P}$ and ${\bf
Q}^{\scriptscriptstyle{\text{\rm g}}}_{\rm P}$ are cyclically monotone, hence
also monotone, over ${\mathbb{R}}^{d}$.

Indeed, similarly to the center-outward distribution functions based on
measure transportation introduced in Section 2.1, geometric distribution
functions are gradients of convex functions. This is clear from (7), which
immediately yields ${\bf F}^{\scriptscriptstyle{\text{\rm g}}}_{\rm
P}=\nabla\phi_{\rm P}^{\scriptscriptstyle{\text{\rm g}}}$, where

${\bf z}\mapsto\phi_{\rm P}^{\scriptscriptstyle{\text{\rm g}}}({\bf z}):={\rm
E}\big{[}\|{\bf z}-{\bf X}\|-\|{\bf X}\|\big{]},\quad{\bf
z}\in{\mathbb{R}}^{d}$

is a convex function and ${\bf X}\sim{\rm P}$. It follows that ${\bf
F}^{\scriptscriptstyle{\text{\rm g}}}_{\rm P}$ is cyclically monotone, which
implies that its inverse ${\bf Q}^{\scriptscriptstyle{\text{\rm g}}}_{\rm P}$
is cyclically monotone too (hence the gradient of a convex function).

However, the geometric distribution function ${\bf F}_{\rm
P}^{\scriptscriptstyle{\text{\rm g}}}$, in general, does not push ${\rm P}$
forward to the spherical uniform ${\rm U}_{d}$ over $\mathbb{B}_{d}$. Would
this be case, ${\bf F}^{\scriptscriptstyle{\text{\rm g}}}_{\rm P}$ and ${\bf
F}^{\scriptscriptstyle\pm}_{\rm P}$ both being gradients of convex functions,
the a.s. uniqueness result of McCann (1955) would imply ${\bf
F}^{\scriptscriptstyle{\text{\rm g}}}_{\rm P}={\bf
F}^{\scriptscriptstyle\pm}_{\rm P}$ a.s. for all ${\rm
P}\in\mathcal{P}_{d}^{+}$ (see Section 2.2), which cannot hold in general.
Indeed, ${\bf F}^{\scriptscriptstyle{\text{\rm g}}}_{\rm P}$ can be defined as
the solution of the previously mentioned linear PDE (see Konen (2022)), while
${\bf F}^{\scriptscriptstyle\pm}_{\rm P}$ can be defined as the solution of a
PDE of the Monge-Ampère type which is nonlinear (see Hallin et al. (2021)).
Therefore, the solutions ${\bf F}^{\scriptscriptstyle{\text{\rm g}}}_{\rm P}$
and ${\bf F}^{\scriptscriptstyle\pm}_{\rm P}$, for general $\rm P$, do not
coincide and ${\bf F}^{\scriptscriptstyle{\text{\rm g}}}_{\rm P}$, in general,
does not push ${\rm P}$ forward to ${\rm U}_{d}$.

Nevertheless, recalling that ${\bf F}^{\scriptscriptstyle{\text{\rm g}}}_{\rm
P}$ is a homeomorphism between ${\mathbb{R}}^{d}$ and $\mathbb{B}_{d}$, we see
that each point ${\bf x}\in{\mathbb{R}}^{d}$ belongs to exactly one geometric
contour, namely ${\bf x}\in\boldsymbol{\mathscr{C}}_{\rm
P}^{\scriptscriptstyle{\text{\rm g}}}(\tau)$ where $\tau=\|{\bf
F}^{\scriptscriptstyle{\text{\rm g}}}_{\rm P}({\bf x})\|$. Therefore, we can
relabel quantile regions so that, in the new labeling, the region associated
with an order $\tau\in[0,1)$ has ${\rm P}$-probability content $\tau$. For
this purpose, it is enough to observe that the map

$\tau\mapsto\kappa_{\rm P}^{\scriptscriptstyle{\text{\rm g}}}(\tau):={\rm
P}[\boldsymbol{\mathscr{R}}_{\rm P}^{\scriptscriptstyle{\text{\rm
g}}}(\tau)],\quad\tau\in[0,1),$

from $[0,1)\to[0,1)$ is a continous and monotone increasing bijection. In
fact, $\kappa_{\rm P}^{\scriptscriptstyle{\text{\rm g}}}$ is the univariate
distribution function of the random variable $\|{\bf
F}^{\scriptscriptstyle{\text{\rm g}}}_{\rm P}({\bf X})\|$ with ${\bf
X}\sim~{}\\!{\rm P}$.Consequently, we can define the relabeled geometric
distribution function of $\rm P$ (of ${\bf X}\sim{\rm P}$) as

${\bf x}\mapsto\tilde{\bf F}_{\rm P}^{\scriptscriptstyle{\text{\rm g}}}({\bf
x}):={\bf B_{\rm P}^{\scriptscriptstyle{\text{\rm g}}}}\big{(}{\bf
F}^{\scriptscriptstyle{\text{\rm g}}}_{\rm P}({\bf x})\big{)},$ (11)

where

${\bf B_{\rm P}^{\scriptscriptstyle{\text{\rm g}}}}({\bf
x}):=\begin{cases}\kappa_{\rm P}^{\scriptscriptstyle{\text{\rm g}}}(\|{\bf
x}\|)\frac{{\bf x}}{\|{\bf x}\|}&\textrm{ if }{\bf x}\neq 0\\\ 0&\textrm{ if
}{\bf x}=0\end{cases}$

is a homeomorphism. Consequently, the relabeled geometric distribution
function $\tilde{\bf F}_{\rm P}^{\scriptscriptstyle{\text{\rm g}}}$ is also a
homeomorphism, hence continuously invertible: define the relabeled geometric
quantile function as its inverse $\tilde{\bf Q}_{\rm
P}^{\scriptscriptstyle{\text{\rm g}}}:=(\tilde{\bf F}_{\rm
P}^{\scriptscriptstyle{\text{\rm g}}})^{-1}$. The collections of regions and
contours associated with $\tilde{\bf F}_{\rm P}^{\scriptscriptstyle{\text{\rm
g}}}$ are the same as those associated with ${\bf
F}^{\scriptscriptstyle{\text{\rm g}}}_{\rm P}$—only the labels have changed.
Call $\tilde{\boldsymbol{\mathscr{R}}}_{\rm P}^{\scriptscriptstyle{\text{\rm
g}}}(\tau):=\tilde{\bf Q}_{\rm P}^{\scriptscriptstyle{\text{\rm
g}}}(\tau\mathbb{B}_{d})$ the relabeled geometric quantile region of order
$\tau$. Then, $\tilde{\boldsymbol{\mathscr{R}}}_{\rm
P}^{\scriptscriptstyle{\text{\rm g}}}(\tau)$ satisfies the essential property
of a quantile region of order $\tau$:

  1. (viii)

the relabeled geometric quantile region $\tilde{\boldsymbol{\mathscr{R}}}_{\rm
P}^{\scriptscriptstyle{\text{\rm g}}}(\tau)$ has ${\rm P}$-probability content
$\tau$ for all $\tau\in~{}\\![0,1)$, irrespective of $\rm P$.

We similarly define the relabeled geometric quantile contours

$\tilde{\boldsymbol{\mathscr{C}}}_{\rm P}^{\scriptscriptstyle{\text{\rm
g}}}(\tau):=\tilde{\bf Q}_{\rm P}^{\scriptscriptstyle{\text{\rm
g}}}(\tau\mathbb{S}^{d-1}),\quad\tau\in[0,1).$

One may wonder whether $\tilde{\bf F}_{\rm P}^{\scriptscriptstyle{\text{\rm
g}}}$ now coincides with ${\bf F}^{\scriptscriptstyle\pm}_{\rm P}$ (hence,
$\tilde{\bf Q}_{\rm P}^{\scriptscriptstyle{\text{\rm g}}}$ with ${\bf
Q}^{\scriptscriptstyle\pm}_{\rm P}$). First, let us stress that, even though
$\|\tilde{\bf F}_{\rm P}^{\scriptscriptstyle{\text{\rm g}}}({\bf X})\|$, just
like $\|{\bf F}^{\scriptscriptstyle\pm}_{\rm P}({\bf X})\|$, is uniformly
distributed over $[0,1)$ when ${\bf X}\sim{\rm P}$, the behavior of ${\bf
F}^{\scriptscriptstyle{\text{\rm g}}}_{\rm P}({\bf X})/\|{\bf
F}^{\scriptscriptstyle{\text{\rm g}}}_{\rm P}({\bf X})\|$, to the best of our
knowledge, remains unknown, and has not been investigated so far. Second,
observe that ${\bf B}_{\rm P}^{\scriptscriptstyle{\text{\rm g}}}$ writes as
the gradient $\nabla h_{\rm P}^{\scriptscriptstyle{\text{\rm g}}}$ of the
convex function

${\bf x}\mapsto h_{\rm P}^{\scriptscriptstyle{\text{\rm g}}}({\bf
x}):=\int_{0}^{\|{\bf x}\|}\kappa_{\rm P}^{\scriptscriptstyle{\text{\rm
g}}}(s)\,ds,\quad{\bf x}\in{\mathbb{R}}^{d}.$

It follows that $\tilde{\bf F}_{\rm P}^{\scriptscriptstyle{\text{\rm
g}}}=\nabla h_{\rm P}^{\scriptscriptstyle{\text{\rm g}}}\circ\nabla\phi_{\rm
P}^{\scriptscriptstyle{\text{\rm g}}}$. However, it is well known that the
composition of two gradients of convex functions, in general, is not the
gradient of a convex function. For these reasons, which do not constitute
rigorous proofs, we believe that there is no obvious reason why the equality
$\tilde{\bf F}_{\rm P}^{\scriptscriptstyle{\text{\rm g}}}={\bf
F}^{\scriptscriptstyle\pm}_{\rm P}$ should hold for general $\rm P$. The
simulations in Section 5) below—which do not constitute mathematical
proofs—actually suggest that $\tilde{\boldsymbol{\mathscr{C}}}_{\rm
P}^{\scriptscriptstyle{\text{\rm g}}}(\tau)$ and
$\boldsymbol{\mathscr{C}}_{\rm P}^{\scriptscriptstyle\pm}(\tau)$, for some
distributions $\rm P$, very strongly differ.

###  4.2 Center-outward quantile regions and contours

Let ${\rm P}\in\mathcal{P}_{d}^{+}$. It is shown in Hallin et al. (2021) that

  1. (i)

the regions $\boldsymbol{\mathscr{R}}_{\rm P}^{\scriptscriptstyle\pm}(\tau)$
are compact and arc-connected for all $\tau\in(0,1)$;

  2. (ii)

the regions $\boldsymbol{\mathscr{R}}_{\rm P}^{\scriptscriptstyle\pm}(\tau)$
are stricly nested as $\tau$ increases: $\boldsymbol{\mathscr{R}}_{\rm
P}^{\scriptscriptstyle\pm}(\tau_{1})\subsetneq\boldsymbol{\mathscr{R}}_{\rm
P}^{\scriptscriptstyle\pm}(\tau_{2})$ for any $0<\tau_{1}<\tau_{2}<1$;

  3. (iii)

the median region $\boldsymbol{\mathscr{R}}_{\rm
P}^{\scriptscriptstyle\pm}(0)$ is compact, convex, and has Hausdorff dimension
at most $(d-1)$;

  4. (iv)

if, for some $k\in\mathbb{N}$ and $\alpha\in(0,1)$, the density of ${\rm P}$
is locally of class $\mathcal{C}^{k,\alpha}$, then the contour
$\boldsymbol{\mathscr{C}}_{\rm P}^{\scriptscriptstyle\pm}(\tau)$ is a
$(d-1)$-dimensional manifold of class $\mathcal{C}^{k+1,\alpha}$ for all
$\tau\in(0,1)$;

  5. (v)

for $d=1$, ${\bf F}^{\scriptscriptstyle\pm}_{\rm P}$, ${\bf
Q}^{\scriptscriptstyle\pm}_{\rm P}$, $\boldsymbol{\mathscr{R}}_{\rm
P}^{\scriptscriptstyle\pm}$, and $\boldsymbol{\mathscr{C}}_{\rm
P}^{\scriptscriptstyle\pm}$ reduce to $F^{\scriptscriptstyle\pm}_{\rm P}$,
$Q^{\scriptscriptstyle\pm}_{\rm P}$, ${\mathscr{R}}_{\rm
P}^{\scriptscriptstyle\pm}$, and ${\mathscr{C}}_{\rm
P}^{\scriptscriptstyle\pm}$.

Note that (iv) was first established in del Barrio et al. (2020) when the
support of ${\rm P}$ is convex, then refined in del Barrio and Gonzáles-Sanz
(2023). By construction, we also have the following:

  1. (vi)

the $\rm P$-probability content of $\boldsymbol{\mathscr{R}}_{\rm
P}^{\scriptscriptstyle\pm}(\tau)$ equals the ${\rm U}_{d}$-probability content
of $\tau\mathbb{B}_{d}$, which is $\tau$: $\boldsymbol{\mathscr{R}}_{\rm
P}^{\scriptscriptstyle\pm}(\tau)$ and $\boldsymbol{\mathscr{C}}_{\rm
P}^{\scriptscriptstyle\pm}(\tau)$, thus, qualify as the center-outward
quantile region and quantile contour of order $\tau$ of $\rm P$;

  2. (vii)

if ${\bf X}\sim{\rm P}$, then ${\bf F}^{\scriptscriptstyle\pm}_{\rm P}({\bf
X})\sim{\rm U}_{d}$; if ${\bf U}\sim{\rm U}_{d}$, then ${\bf
Q}^{\scriptscriptstyle\pm}_{\rm P}({\bf U})\sim{\rm P}$: both ${\bf
F}^{\scriptscriptstyle\pm}_{\rm P}$ and ${\bf Q}^{\scriptscriptstyle\pm}_{\rm
P}$, thus, fully characterize $\rm P$.333Note, however, that the collections
of center-outward regions or center-outward contours alone (that is, without
the mappings ${\bf F}^{\scriptscriptstyle\pm}_{\rm P}$ or ${\bf
Q}^{\scriptscriptstyle\pm}_{\rm P}$ between ${\mathbb{R}}^{d}$ and
$\mathbb{B}_{d}$) do not characterize $\rm P$.

  3. (viii)

the maps ${\bf F}^{\scriptscriptstyle\pm}_{\rm P}$ and ${\bf
Q}^{\scriptscriptstyle\pm}_{\rm P}$ are cyclically monotone, hence monotone,
over ${\mathbb{R}}^{d}$.

Finally, ${\bf F}^{\scriptscriptstyle\pm}_{\rm P}$ and ${\bf
Q}^{\scriptscriptstyle\pm}_{\rm P}$ enjoy the following
invariance/equivariance properties (see Hallin et al. (2022)).

  1. (ix)

Let $\boldsymbol{\theta}\in{\mathbb{R}}^{d}$ and denote by ${\bf O}$ an
arbitrary $d\times d$ orthogonal matrix. Letting ${\bf X}\sim{\rm P}$, write
${\bf F}_{{\bf X}}^{\scriptscriptstyle\pm}$, ${\bf Q}_{{\bf
X}}^{\scriptscriptstyle\pm}$, etc. for ${\bf F}_{{\rm
P}}^{\scriptscriptstyle\pm}$, ${\bf Q}_{\rm P}^{\scriptscriptstyle\pm}$, etc.
Then,

$\displaystyle{\bf F}_{\boldsymbol{\theta}+{\bf O}{\bf
X}}^{\scriptscriptstyle\pm}(\boldsymbol{\theta}+{\bf O}{\bf x})={\bf O}{\bf
F}_{\bf X}^{\scriptscriptstyle\pm}({\bf x})\quad$ $\displaystyle{\bf
x}\in\mathbb{R}^{d},$ $\displaystyle{\bf Q}_{\boldsymbol{\theta}+{\bf O}{\bf
X}}^{\scriptscriptstyle\pm}({\bf O}{\bf u})=\boldsymbol{\theta}+{\bf O}{\bf
Q}_{\bf X}^{\scriptscriptstyle\pm}({\bf u})\quad$ $\displaystyle{\bf
u}\in\mathbb{B}_{d},$
$\displaystyle\boldsymbol{\mathscr{R}}_{\boldsymbol{\theta}+{\bf O}{\bf
X}}^{\scriptscriptstyle\pm}(\tau)=\boldsymbol{\theta}+{\bf
O}\boldsymbol{\mathscr{R}}_{{\bf X}}^{\scriptscriptstyle\pm}(\tau)\quad$
$\displaystyle\tau\in[0,1),$ and
$\displaystyle\boldsymbol{\mathscr{C}}_{\boldsymbol{\theta}+{\bf O}{\bf
X}}^{\scriptscriptstyle\pm}(\tau)=\boldsymbol{\theta}+{\bf
O}\boldsymbol{\mathscr{C}}_{\bf X}^{\scriptscriptstyle\pm}(\tau)\quad$
$\displaystyle\tau\in[0,1).$

In particular, if $\rm P$ enjoys symmetry properties (rotational, central, or
with respect to some hyperplanes), ${\bf F}^{\scriptscriptstyle\pm}_{\rm P}$,
${\bf Q}^{\scriptscriptstyle\pm}_{\rm P}$, $\boldsymbol{\mathscr{R}}_{\rm
P}^{\scriptscriptstyle\pm}(\tau)$, and $\boldsymbol{\mathscr{C}}_{\rm
P}^{\scriptscriptstyle\pm}(\tau)$ enjoy the same symmetry properties as $\rm
P$ for all $\tau\in~{}\\![0,1)$. This is not the case for the distribution and
quantile functions resulting from a transport to the Lebesgue uniform over the
unit cube $[0,1]^{d}$ which are favored by some authors.

These properties, which for dimension $d=1$ reduce to those of
$F^{\scriptscriptstyle\pm}$ and $Q^{\scriptscriptstyle\pm}$ and constitute the
essence of the concepts of distribution and quantile functions, fully justify
the terminology center-outward distribution and quantile functions for ${\bf
F}^{\scriptscriptstyle\pm}_{\rm P}$ and ${\bf Q}^{\scriptscriptstyle\pm}_{\rm
P}$. The only difference, here, between the properties of ${\bf
F}^{\scriptscriptstyle\pm}_{\rm P}$ and ${\bf Q}^{\scriptscriptstyle\pm}_{\rm
P}$ and those of ${\bf F}^{\scriptscriptstyle{\text{\rm g}}}_{\rm P}$ and
${\bf Q}^{\scriptscriptstyle{\text{\rm g}}}_{\rm P}$ is (vi) above, and this
difference disappears if the relabeled functions $\tilde{\bf F}_{\rm
P}^{\scriptscriptstyle{\text{\rm g}}}$ and $\tilde{\bf Q}_{\rm
P}^{\scriptscriptstyle{\text{\rm g}}}$ are substituted for ${\bf F}_{\rm
P}^{\scriptscriptstyle{\text{\rm g}}}$ and ${\bf Q}_{\rm
P}^{\scriptscriptstyle{\text{\rm g}}}$, respectively.

In view of this, the relabeled geometric and the center-outward concepts,
which, it seems, do not coincide, apparently both qualify as equally
satisfactory multivariate concepts of distribution and quantile functions. Is
there a way to break that tie and select one of them rather than the other?
Since analytical expressions are impossible, one has to recur to numerical
approximations—in statistical terms, simulations—in order to perform
comparisons. This is what we are doing in the next section, where comparisons
are based on simulations of the quantile contours
$\boldsymbol{\mathscr{C}}_{\rm P}^{\scriptscriptstyle\pm}(\tau)$ and
$\boldsymbol{\mathscr{C}}_{\rm P}^{g}(\tau)$.

##  5 Geometric and center-outward quantile contours: an empirical comparative
study

###  5.1 Numerical computation of contours

Except for some very special distributions, such as the spherical ones,
geometric and center-outward quantiles in dimension $d\geq 2$, unfortunately,
are sharing the inconvenient absence of analytical expressions. Numerical
solutions are possible, though, based on simulations with high numbers of
replications. For all $N\in\mathbb{N}$, denote by ${\rm P}_{N}$ the empirical
distribution of a random sample of size $N$ drawn from ${\rm P}$. The accuracy
of these numerical solutions relies on consistency results for the empirical
counterparts ${\bf F}_{{\rm P}_{N}}^{\scriptscriptstyle{\text{\rm g}}}$ and
${\bf Q}_{{\rm P}_{N}}^{\scriptscriptstyle{\text{\rm g}}}$ in the geometric
case (obtained by substituting ${\rm P}_{N}$ for ${\rm P}$ in (7)), for ${\bf
F}_{{\rm P}_{N}}^{\scriptscriptstyle\pm}$ and ${\bf Q}_{{\rm
P}_{N}}^{\scriptscriptstyle\pm}$ in the measure-transportation setting (see
Hallin et al. (2021) for precise definitions and details).

In the geometric case, Mottonen et al. (1997) established a Glivenko-Cantelli
result of the form

$\sup_{{\bf z}\in{\mathbb{R}}^{d}}\|{\bf F}_{{\rm
P}_{N}}^{\scriptscriptstyle{\text{\rm g}}}({\bf z})-{\bf
F}^{\scriptscriptstyle{\text{\rm g}}}_{\rm P}({\bf z})\|\to 0$ (12)

with ${\rm P}$-probability one as $N\to\infty$, provided that ${\rm P}$ admits
a bounded density. In fact, one can show that the bounded density assumption
is superfluous: similarly to distribution functions in dimension $d=1$, the
Glivenko-Cantelli result (12) is valid for an arbitrary probability measure.
The same convergence holds for the sequence of relabeled counterparts
$\tilde{{\bf F}}_{\rm P}^{\scriptscriptstyle{\text{\rm g}}}$ of ${{\bf
F}}_{\rm P}^{\scriptscriptstyle{\text{\rm g}}}$ and $\tilde{{\bf F}}_{{\rm
P}_{N}}^{\scriptscriptstyle{\text{\rm g}}}$ of ${{\bf F}}_{{\rm
P}_{N}}^{\scriptscriptstyle{\text{\rm g}}}$; see Konen (2022).

In the measure transportation approach, Hallin et al. (2021) proved a similar
Glivenko-Cantelli result for an interpolated version $\overline{{\bf F}}_{{\rm
P}_{N}}^{\,\scriptscriptstyle\pm}$ of their empirical center-outward
distribution functions ${\bf F}_{{\rm P}_{N}}^{\scriptscriptstyle\pm}$:

$\sup_{{\bf z}\in{\mathbb{R}}^{d}}\|\overline{{\bf F}}_{{\rm
P}_{N}}^{\,\scriptscriptstyle\pm}({\bf z})-{\bf
F}^{\scriptscriptstyle\pm}_{\rm P}({\bf z})\|\to 0$

with ${\rm P}$-probability one as $N\to\infty$, when ${\rm
P}\in\mathcal{P}_{d}^{+}$ (see Section 2.2 for the definition of
$\mathcal{P}_{d}^{+}$).

###  5.2 Numerical comparisons of quantile contours

Figure 1 shows, for orders $\tau\in\\{.25,.50,.75\\}$, plots of the empirical
contours $\tilde{\boldsymbol{\mathscr{C}}}_{{\rm
P}_{N}}^{\scriptscriptstyle{\text{\rm g}}}(\tau)$and
$\boldsymbol{\mathscr{C}}_{{\rm P}_{N}}^{\scriptscriptstyle\pm}(\tau)$ based
on random samples $\mathbb{X}_{N}=\\{{\bf X}_{1},\ldots,{\bf X}_{N}\\}$ of
size $N=2,400$ from four different probability distributions on
${\mathbb{R}}^{2}$. Given the regular grid

$\Big{\\{}{\bf
u}_{k}=\big{(}\cos(\theta_{k}),\sin(\theta_{k})\big{)}:\theta_{k}=\frac{2k\pi}{70}\textrm{
and }1\leq k\leq K\Big{\\}}$

of directions on $\mathbb{S}^{1}$, with $K=70$, each contour
$\tilde{\boldsymbol{\mathscr{C}}}_{{\rm P}_{N}}^{\scriptscriptstyle{\text{\rm
g}}}(\tau)$ is estimated by

$\Big{\\{}\tilde{{\bf Q}}_{{\rm P}_{N}}^{\scriptscriptstyle{\text{\rm
g}}}(\tau{\bf u}_{k}):1\leq k\leq K=70\Big{\\}},$

where $\tilde{{\bf Q}}_{{\rm P}_{N}}^{\scriptscriptstyle{\text{\rm g}}}$
stands for the relabeled geometric quantile map defined in Section 3.1. In
other words, to each level $\tau$ we first associate the level $\tau^{\prime}$
such that

$\tilde{\boldsymbol{\mathscr{C}}}_{{\rm P}_{N}}^{\scriptscriptstyle{\text{\rm
g}}}(\tau)=\boldsymbol{\mathscr{C}}_{{\rm
P}_{N}}^{\scriptscriptstyle{\text{\rm g}}}(\tau^{\prime}).$

Then, the geometric contour $\boldsymbol{\mathscr{C}}_{{\rm
P}_{N}}^{\scriptscriptstyle{\text{\rm g}}}(\tau^{\prime})$ is obtained as the
collection $\\{{\bf z}_{1},\ldots,{\bf z}_{K}\\}$ such that ${\bf z_{k}}$
minimizes the empirical objective function $O_{\tau^{\prime},{\bf
u_{k}}}^{{\rm P}_{N}}$.

Each measure-transportation-based contour is estimated by

$\Big{\\{}{\bf Q}_{{\rm P}_{N}}^{\scriptscriptstyle\pm}(\tau{\bf u}_{k}):1\leq
k\leq K\Big{\\}}.$

To construct the map ${\bf Q}_{{\rm P}_{N}}^{\scriptscriptstyle\pm}$, we first
fix a regular $n_{R}\times n_{S}$ grid ${\mathfrak{G}}_{N}$ over the punctured
unit ball $(0,1)\times\mathbb{S}^{1}$ such that $n_{R}n_{S}=N$;444In the case
of a genuine sample, $N$ is factorized into $N=n_{R}n_{S}+n_{0}$ with
$n_{0}<\min(n_{R},n_{S})$. We refer to Hallin and Mordant (2023) for details
on the choice of $n_{R}$ and $n_{S}$. Here, however, we are dealing with
simulations, and can choose $N$ such that $n_{0}=0$. our simulation relies on
$N=2,400$ replications, with $n_{R}=40$ and $n_{S}=60$. Then, we seek an
optimal coupling between $\mathbb{X}_{N}$ and ${\mathfrak{G}}_{N}$, i.e., we
associate to each ${\bf X}_{k}\in\mathbb{X}_{N}$ a unique ${\bf F}_{{\rm
P}_{N}}^{\scriptscriptstyle\pm}({\bf X}_{k})\in\mathfrak{G}$ in such a way
that

$\sum_{k=1}^{N}\|{\bf X}_{k}-{\bf F}_{{\rm
P}_{N}}^{\scriptscriptstyle\pm}({\bf
X}_{k})\|^{2}=\min_{\pi\in\Pi}\sum_{k=1}^{N}\|{\bf X}_{\pi(k)}-{\bf F}_{{\rm
P}_{N}}^{\scriptscriptstyle\pm}({\bf X}_{k})\|^{2},$

where $\pi$ ranges over the set $\Pi$ of permutations of $\\{1,\ldots,N\\}$. A
continuous map $\overline{\bf F}_{{\rm P}_{N}}^{\,\scriptscriptstyle\pm}$ then
is obtained from the values ${\bf F}_{{\rm
P}_{N}}^{\,\scriptscriptstyle\pm}({\bf X}_{k})$ of ${\bf F}_{{\rm
P}_{N}}^{\,\scriptscriptstyle\pm}$ at ${\bf X}_{k}$ by interpolation
techniques; see Hallin et al. (2021) for furter details.

These samples were generated from

  1. (i)

the centered bivariate normal distribution with covariance matrix
$\Sigma=\begin{pmatrix}2&1\\\ 1&1\end{pmatrix}$;

  2. (ii)

the bivariate distribution with independent exponential marginals with mean
one;

  3. (iii)

the standard skew-$t$ distribution with four degrees of freedom and slant
vector $\alpha=(10,10)$— see Azzalini and Capitanio (2014)—and

  4. (iv)

the banana-shaped distribution considered in Hallin et al. (2021), viz. the
mixture

$\frac{3}{8}\mathcal{N}(\mu_{1},\Sigma_{1})+\frac{3}{8}\mathcal{N}(\mu_{2},\Sigma_{2})+\frac{1}{4}\mathcal{N}(\mu_{3},\Sigma_{3}),$

of three bivariate normal distributions where

$\mu_{1}=\begin{pmatrix}-3\\\ 0\end{pmatrix},\quad\mu_{2}=\begin{pmatrix}3\\\
0\end{pmatrix},\quad\mu_{3}=\begin{pmatrix}0\\\ -\frac{5}{2}\end{pmatrix},$

and

$\Sigma_{1}=\begin{pmatrix}5&-4\\\
-4&5\end{pmatrix},\quad\Sigma_{2}=\begin{pmatrix}5&4\\\
4&5\end{pmatrix},\quad\Sigma_{3}=\begin{pmatrix}4&0\\\ 0&1\end{pmatrix}.$

Figure 1: Geometric (red) and measure-transportation-based (blue) contours of
levels $\tau=.25,.5,.75$ for samples of size $N=2400$ drawn from four
different distributions. Geometric contours are relabeled according to their
probability content, and measure-transportation contours are obtained through
the interpolation scheme (see Hallin et al. (2021)) of the optimal coupling
between the observations and a $40\times 60$ regular grid ${\mathfrak{G}}_{N}$
of the unit ball in ${\mathbb{R}}^{2}$; see Section 4.2 for details. Empirical
quantiles, both geometric and center-outward, are connected by lines in the
plots.

Inspection of Figure 1 reveals the following facts.

  1. (a)

First, the contours obtained by measure transportation seem to better reflect
the geometry of the underlying distributions. This is particularly obvious for
the banana-shaped distribution (iv), where the concave shape of the mixture is
well picked up by the center-outward contours but completely missed by the
relabeled geometric ones. The same phenomenon also takes place, to a less
spectacular extent, under distributions (i)-(iii): the relabeled geometric
quantile region of order .75, under the exponential distribution (ii), largely
extends beyond the exponential support, something the center-outward contours
do not; the same geometric contours, for the skew-$t$ distribution (iii),
contrary to the center-outward ones, hardly exhibit any skewness.

  2. (b)

Second, while relabeled geometric quantiles (as a direct consequence of their
re-labelization rather than because of their L1 characterization) do provide
the desired control over the probability content of their regions, hence
yield, for each ${\bf X}_{i}\in\mathbb{X}_{N}$, a sound concept
$R^{(N)}_{i}\\!\\!\\!\coloneqq\text{rank}(\|{\bf F}_{{\rm
P}_{N}}^{\scriptscriptstyle{\text{\rm g}}}({\bf X}_{i})\|)$ of center-outward
multivariate rank, they do not provide any clear notion of signs. Such signs,
quite on the contrary, are well defined, in the measure-transportation-based
approach, for the sample points ${\bf X}_{i}\in~{}\\!\mathbb{X}_{n}$, as the
unit vectors ${\bf F}_{{\rm P}_{N}}^{\scriptscriptstyle\pm}({\bf
X}_{i})/\|{\bf F}_{{\rm P}_{N}}^{\scriptscriptstyle\pm}({\bf X}_{i})\|$,
$i=1,\ldots,N$, which are i.i.d. (uniform over the unit vectors associated
with the $N$ gridpoints in ${\mathfrak{G}}_{N}$ and independent of the ranks
$R^{(N)}_{1},\ldots,R^{(N)}_{N}$).

  3. (c)

On the other hand, geometric contours are more regular than their measure-
transportation-based counterparts. This is a consequence of (iii) and (iv) in
Sections 3.1 and 3.2, respectively: geometric contours are continuous as soon
as ${\rm P}$ is non-atomic and smooth (at least $\mathcal{C}^{d-1}$) while the
oscillations, hence also the variability, of the interpolated measure-
transportation-based contours are further increased by the interpolation
scheme.

While points (a) and (b) above are pleading for the measure-transportation-
based concepts, the decisive argument in favor of the latter follows from a
theoretical result by Girard and Stupfler (2017) on the “intriguing behavior”
(a clear understatement) of high-order geometric quantile contours.

###  5.3 Numerical comparison of extreme contours

Girard and Stupfler (2017) establish the following disturbing result for
extreme geometric quantiles: if ${\rm P}$ is a probability measure on
${\mathbb{R}}^{d}$ with finite second-order moments, i.e. ${\rm E}_{\rm
P}[\|{\bf X}\|^{2}]<\infty$ for ${\bf X}\sim{\rm P}$ with covariance matrix
$\boldsymbol{\Sigma}$, then, denoting by Tr$(\bf A)$ the trace of a square
matrix $\bf A$,

$\|{\bf Q}^{\scriptscriptstyle{\text{\rm g}}}_{\rm P}(\tau{\bf
u})\|^{2}(1-\tau)\to\frac{1}{2}\Big{(}\operatorname{Tr}{\boldsymbol{\Sigma}}-{\bf
u}^{\prime}{\boldsymbol{\Sigma}}{\bf u}\Big{)}$

as $\tau\to 1$ for any ${\bf u}\in\mathbb{S}^{d-1}$. First, let us stress that
the quantity $\operatorname{Tr}{\boldsymbol{\Sigma}}-{\bf
u}^{\prime}{\boldsymbol{\Sigma}}{\bf u}$ is always non-negative, and is
strictly positive as soon as $\boldsymbol{\Sigma}$ is non-degenerate. In
addition to providing the rate at which extreme geometric quantiles are
“escaping to infinity,” this result also provides valuable insights on the
shape of the corresponding extreme contours. For instance, taking for ${\bf
u}$ a unit eigenvector ${\bf u}_{\text{\rm max}}$ associated with
$\boldsymbol{\Sigma}$’s largest eigenvalue yields a norm $\|{\bf
Q}^{\scriptscriptstyle{\text{\rm g}}}_{\rm P}(\tau{\bf u}_{\text{\rm max}})\|$
for the contour $\boldsymbol{\mathscr{C}}_{\rm
P}^{\scriptscriptstyle{\text{\rm g}}}(\tau)$ of order $\tau$ in direction $\bf
u$ which is minimal for $\tau\to 1$; the same contour, in the direction ${\bf
u}_{\text{\rm min}}$ associated with the smallest eigenvalue of
$\boldsymbol{\Sigma}$, asymptotically yields the largest norm $\|{\bf
Q}^{\scriptscriptstyle{\text{\rm g}}}_{\rm P}(\tau{\bf u}_{\text{\rm
min}})\|$. More generally, due to the fact that the shape factor
$\operatorname{Tr}{\boldsymbol{\Sigma}}-{\bf
u}^{\prime}{\boldsymbol{\Sigma}}{\bf u}$ is a decreasing function of ${\bf
u}^{\prime}{\boldsymbol{\Sigma}}{\bf u}$, the asymptotic behavior of extreme
quantiles—moving away fast from the center along the eigendirections of
$\boldsymbol{\Sigma}$ associated with small eigenvalues, and slowly along the
eigendirections of $\boldsymbol{\Sigma}$ associated with large eigenvalues—is
exactly the opposite of what it is expected to be. This is highly
counterintuitive and, actually, totally unacceptable.

Figure 2 provides a finite-sample empirical illustration of this theoretical
fact. We generated a sample of size $N=1,000$ from a centered bivariate
Gaussian distribution $\mathcal{N}({\boldsymbol{0}},{\boldsymbol{\Sigma}})$
with covariance matrix $\boldsymbol{\Sigma}=\begin{pmatrix}1/8&0\\\
0&3/4\end{pmatrix}.$ Then, for $\tau\in\\{.90,.95,.99\\}$, we computed the
empirical relabeled geometric and transportation-based quantile contours of
order $\tau$, as described in Sections 5.1 and 5.2, with $K=50$, $n_{R}=20$,
and $n_{S}=50$ for the latter (see Section 5.2 for notation).

Figure 2: Geometric (red) and measure-transportation-based (blue) contours of
order $\tau=.90,.95,.99$ for samples of size $N=1,000$ drawn from a
(nonspherical) bivariate Gaussian distribution on ${\mathbb{R}}^{2}$.
Empirical quantiles, both geometric and center-outward, are connected by lines
in the plots.

Figure 2 perfectly illustrates the pathological phenomenon described by Girard
and Stupfler (2017): the geometric contours are much wider in the (horizontal)
directions associated with $\boldsymbol{\Sigma}$’s smaller eigenvalue, and
narrower in the (vertical) directions associated with $\boldsymbol{\Sigma}$’s
larger eigenvalue. This pathological behaviour is increasingly severe as
$\tau$ approaches one: for $\tau=.90$, the contour is approximately circular.
In sharp contrast, the center-outward contours, irrespective of $\tau$, are
correctly accounting for the shape of the distribution and the relative
magnitudes of $\boldsymbol{\Sigma}$’s eigenvalues.

##  6 Conclusions

Among all concepts of multivariate quantiles proposed in the literature,
geometric quantiles (after due re-labelization) and measure-transportation-
based center-outward quantiles are the most convincing ones: they both reduce
to classical concepts in dimension one, and both satisfy most properties
expected from quantiles, quantile regions, and quantile contours. Empirical
comparisons and a highly pathological property of extreme geometric quantiles
established by Girard and Stupfler (2017), however, are tilting the balance in
favour of the measure-transportation-based concept.

## Acknowledgments

Marc Hallin acknowledges the support of the Czech Science Foundation grant
GAČR22036365. Dimitri Konen is supported by a research fellowship from the
Centre for Research in Statistical Methodology (CRiSM) of the University of
Warwick. We thank Eustasio del Barrio and Alberto González-Sanz for kindly
providing their codes for the computation of center-outward quantile functions
and contours.

## References

  * Azzalini and Capitanio (2014) Azzalini, A. and A. Capitanio (2014).  The Skew-Normal and Related Families.  IMS Monograph series. Cambridge University Press.
  * Chaudhuri (1996) Chaudhuri, P. (1996).  On a geometric notion of quantiles for multivariate data.  Journal of the American Statistical Association 91, 862–872.
  * del Barrio and Gonzáles-Sanz (2023) del Barrio, E. and A. Gonzáles-Sanz (2023).  Regularity of center-outward distribution functions in non-convex domains.  arXiv preprint arXiv:2303.16862.
  * del Barrio et al. (2020) del Barrio, E., A. González-Sanz, and M. Hallin (2020).  A note on the regularity of optimal-transport-based center-outward distribution and quantile functions.  Journal of Multivariate Analysis 180, 104671.
  * Figalli (2018) Figalli, A. (2018).  On the continuity of center-outward distribution and quantile functions.  Nonlinear Analysis 177, 413–21.
  * Girard and Stupfler (2017) Girard, S. and G. Stupfler (2017).  Intriguing properties of extreme geometric quantiles.  REVSTAT 15, 107–139.
  * Hallin et al. (2021) Hallin, M., E. del Barrio, J. Cuesta-Albertos, and C. Matrán (2021).  Distribution and quantile functions, ranks and signs in dimension $d$: a measure transportation approach.  Annals of Statistics 49, 1139–1165.
  * Hallin and Mordant (2023) Hallin, M. and G. Mordant (2023).  On the finite-sample performance of measure-transportation-based multivariate rank tests.  In M. Yi and K. Nordhausen (Eds.), Robust and Multivariate Statistical Methods: Festschrift in Honor of David E. Tyler, pp. 87–119. Springer.
  * Koltchinski (1997) Koltchinski, V. (1997).  M-estimation, convexity and quantiles.  Annals of Statistics 25, 435–477.
  * Konen (2022) Konen, D. (2022).  Recovering a probability measure from its geometric cdf.  arxiv preprint arXiv:2208.11551.
  * Konen and Paindaveine (2022) Konen, D. and D. Paindaveine (2022).  Multivariate $\rho$-quantiles: a spatial approach.  Bernoulli, 28, 1912–1934.
  * McCann (1995) McCann, R. J. (1995).  Existence and uniqueness of monotone measure-preserving maps.  Duke Mathematical Journal 80, 309–323.
  * Mottonen et al. (1997) Mottonen, J., H. Oja, and J. Tienari (1997).  On the efficiency of multivariate spatial sign and rank tests.  Annals of Statistics 25, 542–552.
  * Paindaveine and Virta (2020) Paindaveine, D. and J. Virta (2020).  On the behavior of extreme $d$-dimensional spatial quantiles under minimal assumptions.  In A. Daouia and A. Ruiz-Gazen (Eds.), Advances in Contemporary Statistics and Econometrics, pp. 243–259. Springer.
  * Rockafellar (1970) Rockafellar, R. T. (1970).  Convex Analysis.  Princeton Mathematical Series. Princeton University Press.