# Massive parallelization of projection-based depths

Leonardo Leone, Pavlo Mozharovskyi, David Bounie

► To cite this version:

Leonardo Leone, Pavlo Mozharovskyi, David Bounie. Massive parallelization of projection-based depths. 2025. hal-05109296

HAL Id: hal-05109296

https://telecom-paris.hal.science/hal-05109296v1

Preprint submitted on 12 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Massive parallelization of projection-based depths

Leonardo Leone

CREST, CNRS, LTCI, Télécom Paris, Institut Polytechnique de Paris

Pavlo Mozharovskyi

LTCI, Télécom Paris, Institut Polytechnique de Paris

David Bounie

CREST, CNRS, Télécom Paris, Institut Polytechnique de Paris

June 09, 2025

## Abstract

This article introduces a novel methodology for the massive parallelization of projection-based depths, addressing the computational challenges of data depth in high-dimensional spaces. We propose an algorithmic framework based on Refined Random Search (RRS) and demonstrate significant speedup (up to 7,000 times faster) on GPUs. Empirical results on synthetic data show improved precision and reduced runtime, making the method suitable for large-scale applications. The RRS algorithm (and other depth functions) are available in the **Python-library data-depth** with ready-to-use tools to implement and to build upon this work.

**Keywords:** Data Depth, Parallel Computing, Projection-Based Depths, Refined Random Search, GPU Acceleration.
# 1 Introduction

The concept of data depth has become a cornerstone in multivariate statistics, offering a non-parametric measure of centrality and outlyingness in high-dimensional spaces. Introduced by Tukey [1975], data depth functions quantify the centrality of a point within a data set, providing a robust alternative to traditional distance-based metrics. However, the computational complexity of calculating data depth, especially in high-dimensional settings, has hindered its practical application. Current work addresses these challenges by introducing a novel methodology for the massive parallelization of projection-based depths, significantly enhancing computational efficiency and scalability.

Projection-based depths, such as halfspace depth, projection depth, and asymmetric projection depth, are particularly well-suited for high-dimensional data analysis due to their affine invariance and robustness. These depths satisfy the weak projection property, allowing multivariate depths to be computed as the minimum of univariate depths across projections in all directions. However, exploring all possible directions in high-dimensional spaces is computationally unfeasible, necessitating optimization routines like the Refined Random Search (RRS). While RRS has shown promise, it is often limited by sequential computation, failing to fully leverage modern parallel computing architectures.

In this article, we propose an algorithmic framework that harnesses GPU acceleration to achieve massive parallelization of projection-based depths. Our approach builds upon the RRS algorithm and extends it by decomposing the computation into highly parallelizable tasks suited for GPU architecture. First, we formalize characteristic time equations for both our parallel approach ($T_G$) and the standard approach ($T_C$), and derive speedup equations ($S_p$) for RRS inspired by Amdahl's law [Amdahl, 2007]. The theoretical speedup $S_p$ is expressed as a quotient involving $\lambda \cdot d$, where $\lambda$ is the ratio of CPU to GPU frequency ($\lambda = \frac{f_{CPU}}{f_{GPU}}$), $d$ is the dimension, and $g$ is the number of available GPU cores. This reevaluation of Amdahl's law highlights speedup limitations due to core speed, core quantity, and sequential operations.

Second, we decompose the computation of projection-based depths into parallelizable tasks, achieving significant speedups of up to 7,000 times faster than state-of-the-art methods, particularly when implemented on GPUs like the NVIDIA GeForce RTX 4080 (see [Dyckerhoff and Mozharovskyi, 2016, Fojtík et al., 2023], and for projection depth [Zuo and Lai, 2011]). The key innovation lies in the parallelization of the three main operations involved in depth computation: direction generation, data projection, and univariate depth calculation. By breaking these operations into simple, *embarrassingly parallel* tasks, we achieve efficient use of GPU resources, resulting in dramatic reductions in runtime without sacrificing precision.
Third, we provide insights into hyperparameter fine-tuning, such as the selection of spherical cap shrinkage, the number of refinements, and directions, to optimize convergence and runtime. These findings ensure the robustness and adaptability of our method across different data sets and dimensionalities. The center-outward ordering precision, evaluated through ranking correlation techniques, confirms the reliability of our approach in high-dimensional data analysis.

Our contributions are multifold. First, we introduce a novel algorithmic framework for the massive parallelization of projection-based depths, achieving significant speedups on GPUs. Second, we provide a comprehensive analysis of the computational challenges associated with data depth and demonstrate the effectiveness of our approach in overcoming these challenges. Third, we explore hyperparameter fine-tuning for the RRS algorithm, providing insights into optimal parameter selection for large-scale scenarios. Finally, we conduct a thorough runtime analysis, demonstrating the empirical speedup achieved by our method and validating our theoretical model. These contributions not only advance the state-of-the-art in data depth computation but also pave the way for broader adoption of depth-based methods in practical applications. By overcoming the computational bottlenecks, we enable researchers and practitioners to leverage the full potential of data depth in large-scale, high-dimensional settings, opening new avenues for robust statistical analysis.

The remainder of the paper is organized as follows. Section 2 provides a detailed background on data depth, focusing on projection-based depths and their computational challenges. Section 3 introduces our methodology for massive parallelization, detailing the parallelization of each operation and the theoretical underpinnings of our approach. Section 4 presents the complete massively-parallel algorithm, including the integration of parallel directions generation, projection, and depth computation. In Section 5, we explore hyperparameter fine-tuning for the RRS algorithm, providing insights into optimal parameter selection for large-scale scenarios. Section 6 conducts a thorough runtime analysis, demonstrating the empirical speedup achieved by our method and validating our theoretical models. Finally, Section 7 evaluates the precision of our approach through center-outward ordering and depth benchmarking, confirming its effectiveness in high-dimensional data analysis. Section 8 concludes.

# 2 Background on data depth

In this section, we review the most relevant aspects of data depth from an optimization perspective. For simplicity, we focus on the empirical version of data depth throughout this article. Following Tukey's seminal idea [Tukey, 1975], data depth is a statistical function that measures the degree to which an arbitrary observation belongs to a data set. Consider a data set $X = \{x_1, ..., x_n\} \subset \mathbb{R}^d$
(using set notation for convenience, though ties are generally not excluded) and a point $z \in \mathbb{R}^d$.

**Definition 1.** Data depth is a function

$$D : \mathbb{R}^d \times \mathbb{R}^{n \times d} \to [0, 1], (z, \mathbf{X}) \mapsto D(z|\mathbf{X}),$$

which satisfies the properties of affine invariance, vanishing at infinity, monotonicity on rays originating from the deepest point, and upper-semicontinuity.

For a comprehensive and formal discussion on data depth and its axiomatics, we refer the reader to Zuo and Serfling [2000] and Mosler and Mozharovskyi [2022]. In this article, we focus on a property that is particularly relevant to our current analysis: the weak projection property. Let $\mathbb{S}^{d-1} = \{\mathbf{x} \in \mathbb{R}^d : \|\mathbf{x}\| = 1\}$ denote the unit sphere and $\mathbb{T}$ be the transposition operator.

**Definition 2 (Dyckerhoff [2004]).** The function $D$ from Definition 1 satisfies the weak projection property if for any $z \in \mathbb{R}^d$, $\mathbf{X} \subset \mathbb{R}^d$, it holds:

$$D(z|\mathbf{X}) = \inf_{\mathbf{u} \in \mathbb{S}^{d-1}} D(z^\top \mathbf{u} | \mathbf{X}^\top \mathbf{u}).$$

In the following, we will omit the term “weak”, and simply refer to depths that satisfy the projection property. These include Mahalanobis [Mahalanobis, 2018], halfspace [Tukey, 1975], zonoid [Koshevoy and Mosler, 1997], projection [Zuo and Serfling, 2000], asymmetric projection [Dyckerhoff, 2004], (continuous) expected convex hull [Cascos, 2007], geometrical [Dyckerhoff and Mosler, 2011] depths. Definition 2 leads to the following computationally important corollary:

**Corollary 1.** Let $\mathbf{u}_1, \dots, \mathbf{u}_k \subset \mathbb{S}^{d-1}$, then:

$$D(z|\mathbf{X}) \le \min_{j=1,\dots,k} D(z^\top \mathbf{u}_k | \mathbf{X}^\top \mathbf{u}_k).$$

In other words, multivariate data depths can be computed as the minimum of univariate depths across projections in all directions. Since exploring all possible directions (an infinite set) is computationally infeasible, a finite subset must be examined using an optimization routine. Within this routine, only univariate depths need to be computed.

For any $\mathbf{u} \in \mathbb{S}^{d-1}$, let $\mathbb{R} \ni y = z^\top \mathbf{u}$ and $\mathbb{R}^n \ni \mathbf{y} = \mathbf{X}^\top \mathbf{u} = \{\mathbf{x}_1^\top \mathbf{u}, \dots, \mathbf{x}_n^\top \mathbf{u}\}$, such that a depth satisfying the projection property can be written as $\inf D(y|\mathbf{y})$. A depth satisfying the projection property can then be expressed as $\inf D(y|\mathbf{y})$.

Among the various depths that satisfy the projection property, we will focus on three robust and geometrically affine-invariant notions throughout this article:
• Halfspace depth:

$$D_H(y|\mathbf{y}) = \frac{1}{n} \min\{\#(\mathbf{y} \le y), \#(\mathbf{y} \ge y)\},$$

where # denotes the cardinal of a set, with the vector-scalar operation applicable element-wise.

• Projection depth:

$$D_P(y|\mathbf{y}) = \left(1 + \frac{|y - \text{med}(\mathbf{y})|}{\text{MAD}(\mathbf{y})}\right)^{-1},$$

where med(\mathbf{y}) denotes the (usual) univariate median and MAD(\mathbf{y}) = med(|\mathbf{y} - med(\mathbf{y})|) stands for the median absolute deviation for the median.

• Asymmetric projection depth:

$$D_{AP}(y|\mathbf{y}) = \left(1 + \frac{(y - \text{med}(\mathbf{y}))_+}{\text{MAD}_+(\mathbf{y})}\right)^{-1},$$

where $(a)_+ = \max\{a, 0\}$ and MAD$_+(\mathbf{y})$ is the median of the positive deviations from the median.

It is noteworthy that all three of the aforementioned univariate notions are computable with the best achievable complexity of $O(n)$. However, implementation details may reveal important differences, as we will discuss in Section 4.

## 3 Massive Parallel Computation of Data Depth

We propose a viable methodology to compute the depth of a point considering a large number of directions and samples in the data set.

A projection-based notion follows three (main) operations to compute the final depth:

1. *Generating k directions in $\mathbb{S}^{d-1}$.*

2. *Projecting the data onto these directions.*

3. *Computing the univariate depth within each projection.*

The time importance of each operation is illustrated in Figure 1. Here, directions are randomly generated following a uniform distribution in $\mathbb{S}^{d-1}$ (i.e., without refinements). The univariate depth is computed using three depth notions: $D_H$ in blue lines, $D_P$ in red lines, and $D_{AP}$ in black lines.

The left graph shows the time percentage of each operation with respect to the number of directions $(k)$, with the space dimension $d$ fixed at 150 and the data set $\mathbf{X}$ containing 1,000 samples.
<figure><img src="image_2.png" /><figcaption>Figure 1: The time breakdown for computing projection depth on CPU without parallelization is illustrated as follows: blue lines represent the time spent generating directions in the unit sphere, following a normal distribution; red lines represent the time spent projecting onto these directions; black lines represent the time spent computing the univariate (projection) depth. Left graph: Displays the percentage of time spent on each operation, with *n* (number of points) fixed at 1,000 and *d* (space dimension) to 150. Right graph: Displays the percentage of time spent on each operation with respect to the space dimension *d*, with *n* and *k* (number of directions) fixed at 1,000.</figcaption></figure>

We note that the time percentage remains nearly constant for all *k*, as the time for both main operations (direction generation and projection) grows proportionally with the number of directions. Overall, projection operations dominate the runtime, being the most time-consuming step.

The right graph, using 1,000 directions and a data set with 1,000 samples, examines the effect of the space dimension on time percentage. Here, $D_P$ and $D_{AP}$ show a greater emphasis on univariate depth computation, as these notions rely on (practically) more time-consuming operations like the median. In contrast, $D_H$ relies solely on value comparisons, making it less computationally intensive in this context.

The time complexity of these three operations depends on the number of points *n*, directions *k*, and space dimensions *d*, reflecting their computational significance.

1. Generating directions: This operation has a time complexity of $O(kd)$ and is the least computationally intensive among the three, as it involves generating *k* random directions in a *d*-dimensional space.

2. Projecting data: This step has a time complexity of $O(knd)$, which grows with the space dimension *d* due to the tensor multiplications required to project the data onto the generated directions. This operation is more computationally demanding than generating directions.

3. Computing univariate depths: The time complexity of this operation depends on the chosen depth notion. For instance, in the case of projection depth ($D_P$), the complexity is $O(kn)$. For
a generic depth notion, the complexity is denoted as $O(D)$. This step can be computationally significant, especially for depth notions that involve more complex calculations, such as the median.

With a resulting time complexity of $O(kd) \cup O(knd) \cup O(D)$, the objective of this article is to develop a method for massive parallelization of the rate-determining operation, $O(knd)$, using an optimization routine called Refined Random Search (RRS). This approach aims to significantly reduce computational time while maintaining the precision of the algorithm.

Several optimization methods have been used to identify the direction that minimizes the depth of a point, including (best performing) RRS, coordinate descent, and the Nelder-Mead method [Dyckerhoff et al., 2021]. Among these methods, we selected RRS due to its balance of precision and parallelization potential, making it the most suitable for constructing our framework methodology. This algorithm is further explored in Section 4.

Our approach leverages the complete parallel computation of depth, including both projection and univariate depth calculation, for a point with respect to different directions within the same refinement. By focusing on data projection, this approach allows our framework to decompose $O(knd)$ into synchronous, parallel calculations.

Let us introduce our approach with the following model. All mentioned operations result in a characteristic time $T_C$, where each operation has an associated cost and respective time complexity. We denote the total number of directions per refinement as $m$, where $m = \lfloor \frac{k}{r} \rfloor$, $T_C$ can be written as:

$$T_C = C_c + r \cdot (C_{rv} \cdot m \cdot d + C_p \cdot m \cdot n \cdot d + C_{d(1)} \cdot D), \quad (1)$$

where $C_c$ is a constant computational cost, $C_{rv}$ is the cost of generating random directions, $C_p$ is associated with projecting data, $C_{d(1)}$ represents the additional multiplicative factor of the univariate depth computation. These operations have distinct computational characteristics and are parallelized in different ways. Our goal is to break down each of these steps into multiple simple tasks suitable for massive parallel (e.g., with GPU) computation.

Most of these operations can be reduced to tensor multiplication(s), floating-point operations (FLOPS), and conditional statements. In parallel computing, these operations are referred to as *embarrassingly parallel*, as most of their steps can be performed concurrently [Navarro et al., 2014]. In *g*-massive parallel computation (e.g., in a GPU with *g* computing units):

* generating random directions can be performed in $\lfloor \frac{m \cdot d}{g} \rfloor$ steps,

* depth computation can be completed in $\lfloor \frac{D}{g} \rfloor$ steps,
• projecting data can be executed in $\lfloor \frac{m \cdot n \cdot d}{g} \rfloor$ steps.

Let us denote with the term $\lfloor \frac{d}{d(m,n)_{max}} \rfloor$ the breakdown of the inner product $(xu^\top)$. Depending on the size of the vectors, vector multiplication might be further divided into smaller kernels to accommodate memory usage:

$$xu^\top = \sum_{i=1}^{d} x_i \cdot u_i = \sum_{i=1}^{p} x_i \cdot u_i + \sum_{j=p+1}^{d} x_j \cdot u_j.$$

Further, let $\lambda$ denote the ratio between the massive parallel frequency (e.g., a GPU core frequency) and the basic (non-parallel) frequency (e.g., a CPU core frequency) to be considered in the parallel time equation. With the frequency ratio $\lambda = \frac{f_{CPU}}{f_{GPU}}$, characteristic time for our parallel approach $T_G$ can be written as:

$$T_G = C_c + r \cdot \lambda \cdot \left( C_{rv} \cdot \left\lfloor \frac{m \cdot d}{g} \right\rfloor + C_p \cdot \left\lfloor \frac{d}{d_{max(m,n)}} \right\rfloor \cdot \left\lfloor \frac{m \cdot n}{g} \right\rfloor + C_{d(1)} \cdot \left\lfloor \frac{D}{g} \right\rfloor \right). \quad (2)$$

Following Amdahl's law for parallel speedup, we can construct a theoretical time ratio to model the behavior of our approach in a system with $g$ massively-parallel (e.g., GPU) cores [Amdahl, 2007, Gustafson, 1988, Marowka, 2012]. Speedup ($S_p$) is formalized as the ratio between characteristic times, expressed as:

$$S_p = \frac{T_C}{T_G}. \quad (3)$$

After formalizing these equations, we now propose an optimized algorithmic framework for massively-parallel depth parallelization.

# 4 Algorithmic framework for massive parallelization

Refined Random Search algorithm combines purely random approximation, using independent directions uniformly distributed over the unit sphere, with an optimization element. Specifically, this optimization routine consists of three steps, repeated within each refinement. A brief explanation of the main operation is provided in Algorithm 1. For more details on the RRS algorithm, see Dyckerhoff et al. [2021] and Section A in the Supplementary Material.

Step 3 of Algorithm 1 presents the pseudo-code for the three main operations performed within each refinement: generation, projection, and computing of univariate depth. At the end of each refinement, a new pole $p$ is selected—the direction with the minimal computed depth in the current refinement step. This selection reduces the search region by an angle $\alpha$, a process known as spherical cap shrinkage.
**Algorithm 1** Classical refined random search

**Input:** $\mathbf{X} \in \mathbb{R}^{n \times d}$, $\mathbf{z} \in \mathbb{R}^d$, $k, r, \alpha$.

**Step 1.** Set $\mathbf{p} = \mathbf{e}_1$, $D_{min} = 1$, $m = \lfloor \frac{k}{r} \rfloor$, $\epsilon = \frac{\pi}{2}$.

**Step 2.** Loop for each refinement: For $l = 1$ to $r$:

**Step 3.** Minimize depth in the $\epsilon$-region: For $j = 1$ to $m$:

(a) **Generating** direction: $\mathbb{R}^d \ni \mathbf{u}_j \sim \mathcal{U}(\mathbb{S}^{d-1})$ inside an $\epsilon$-region around $\mathbf{p}$.

(b) **Projecting** the data: $\mathbf{y}_j = \mathbf{X}^\top \mathbf{u}_j$ and $\mathbf{y}_j = \mathbf{z}^\top \mathbf{u}_j$.

(c) **Computing** univariate depth: $D_j = D(\mathbf{y}_j|\mathbf{y}_j)$.

**Step 4.** Minimal depth $D_l = \min_{j \in \{1,...,m\}} D_j$ and its respective direction $\mathbf{u}_l$.

**Step 5.** Reevaluate depth: if $D_l < D_{min}$ update values:

(a) Update minimal depth $D_{min} = D_l$.

(b) Update pole $\mathbf{p} = \mathbf{u}_l$.

**Step 6.** Reduce the size of the region $\epsilon = \epsilon \cdot \alpha$.

**Output:** $D_{min}$.

An illustration of the RRS algorithm for dimension $d = 3$ is shown in Figure 2, where four refinements are necessary to find the minimal depth. In this figure, each point represents a unique direction $\in \mathbb{S}^{d-1}$, with colors' gamma indicating the depth values. Points marked with a cross correspond to the poles, and solid lines represent the borders of the respective spherical caps. Poles and caps can be matched by color, while the cross of the “new” color represents the new pole.

A crucial observation enables massive parallelism: while each refinement depends on the previous one, the three operations within the same refinement – generation, projection, and computation – are entirely independent for each direction and can thus be run in parallel. Given that the number of random directions within the same refinement step can be substantial, this allows for dramatic acceleration (see the above Section 3 for more details, cf. $T_C$ in (1) and $T_G$ in (2)).

## 4.1 A parallel perspective on the Refined Random Search algorithm

The goal is not merely to identify which parts of Algorithm 1 should be parallelized, but to develop a viable algorithm with massive parallelism. To achieve this, we classify the operations into sequential and parallel components. As illustrated in Figure 1, during purely sequential execution, the most time-consuming operation (with time complexity $O(knd)$, namely projection) is highly parallelizable within each refinement. In this part, we focus on parallelizing the projection operation. The comprehensive algorithm, including the parallelization of further operations, is discussed in the following Subsection 4.2.

Let $\mathbf{U} = (\mathbf{u}_1, ..., \mathbf{u}_m)^\top$ be a $m \times d$ matrix consisting of directions $\mathbb{R}^d \ni \mathbf{u}_1, ..., \mathbf{u}_m \in \mathbb{S}^{d-1}$ (as rows). With the data set denoted by $\mathbf{X}$, the operation to parallelize is the matrix multiplication
<figure><img src="image_3.png" /><figcaption>Figure 2: RRS illustration. Points marked with "X" indicate poles, round points represent random directions, and solid lines indicate the limits of the search region. The color of the round points shows the local $D_P$. The pole with the same color as the boundary is the center of the current region while the pole with a different color is the center of the following refinement.</figcaption></figure>

$P = U X^\top$, where $P$ is a $m \times n$ matrix where each row corresponds to the projection of all observations in $X$ onto a direction $u_j, j = 1, ..., m$.

Algorithm 2 provides a simplified view of data set projection $XU^\top$ as a matrix multiplication, which can be performed in parallel according to hardware constraints [Navarro et al., 2014]. Note that the only non-parallel(izable) loop (over $l$) requires values from previous iteration(s) for the final sum, an operation known as **reduce**, which is also potentially parallelizable. In what follows, $(U_{j,\cdot})^\top$ (and similar) naturally denote the (column) vector corresponding to pointwise projection of $X$ onto $u_j$, i.e., $j^{th}$ row of matrix $U$.

<caption>**Algorithm 2 Parallel projection**</caption>

function PROJPARALLEL($X$, $U$)
    **Initialisation**
    $P \leftarrow 0_{m \times n}$
    **Step 1:** Parallel loop through directions
    for $j \leftarrow 1$ to $m$ do (parallel)
        **Step 2:** Loop through all points
        for $i \leftarrow 1$ to $n$ do (parallel)
            for $l \leftarrow 1$ to $d$ do
                $P_{j,i} = P_{j,i} + X_{i,l} \cdot U_{j,l}$
        return $P$
<caption>**Algorithm 3** Massively-parallel projection</caption>

**function** PROJMASSIVE($X$, $U$, $\kappa$)
**Initialisation**
$\mathcal{B} \leftarrow \emptyset$
$\mathbf{P} \leftarrow \mathbf{0}_{m \times n}$
Load $X$, $U$ into the (GPU's) shared memory
$t \leftarrow 0$
$b \leftarrow 0$
**Step 1:** Create threads
**for** $i \leftarrow 1$ **to** $n$ **do**
    **for** $j \leftarrow 1$ **to** $m$ **do**
        Create thread $\mathcal{t} : P_{j,i} \leftarrow X_{i,\cdot}(U_{j,\cdot})^{\top}$
        if $t \bmod \kappa = 0$ **then**
            $t \leftarrow 0$
            $b \leftarrow b + 1$
            Add $\mathcal{t}$ to $\mathcal{B}[b]$
            $t \leftarrow t + 1$
    **end for**
    Add $\mathcal{t}$ to $\mathcal{B}[b]$
**Step 2:** Run threads in parallel (on GPU)
**for** $k \leftarrow 1$ **to** $b$ **do**
    Load $\mathcal{B}[k]$ into available streaming multiprocessor
    Start executing (in parallel) threads of $\mathcal{B}[k]$
    **while** no streaming multiprocessor is available **do**
        wait for a free streaming multiprocessor
**return** $\mathbf{P}$

▷ Queue of threadblocks
▷ Thread count
▷ Threadblock count
▷ Loop over the data set
▷ Loop over the directions
▷ Reset counting for another threadblock

While the main focus is on data projection, the operations of direction generation and computation of univariate depth can also be parallelized. With a primary interest in massive parallelization, preferably using GPU architectures, algorithmic tasks should be broken down into simple algebraic operations translatable into GPU-routine instructions. The breakdown of the projection task for further processing on a GPU is illustrated in Figure 3, and its implementation is suggested in Algorithm 3.

For relevant details on GPU architecture, the reader is referred to Section C of the Supplementary Material. Here, we mention several crucial details. Even with powerful infrastructure, resources have hierarchies and limitations. The hierarchy is expressed by threadblocks (assembling multiple threads), and the computational limitation (excluding memory limitations for simplicity) is represented by the maximum number of threadblocks that can run simultaneously within a streaming multiprocessor (see Figure 3).

## 4.2 The complete massively-parallel algorithm

Algorithm 4 integrates the previously presented logic into a single algorithm for the massive parallelization of depth computations that satisfy the projection property. The function ProjMassive can be easily applied to a single point, for example, if it is represented as a $1 \times d$ matrix. Not
<figure><img src="image_4.png" /><figcaption>Figure 3: Thread creation and parallel projection of $XU^T$ on the GPU.</figcaption></figure>

only can the projection operation be parallelized, but the generation of directions (see Section C in the Supplementary Material for more details and a precise algorithm) and the computation of univariate depths can also be parallelized in a similar manner. Algorithm 4 is generic with respect to the depth notion, as long as it satisfies the projection property.

# 5 Hyperparameter fine-tuning

The second part of our methodology focuses on exploring hyperparameters for RRS in large-scale scenarios, addressing a gap in the literature. Fine-tuning is necessary for two main reasons: to ensure the best convergence of the algorithm to the lowest depth and to optimize the convergence for the shortest runtime. For RRS, we study the optimization based on three hyperparameters: the spherical shrinkage coefficient $\alpha$, the number of refinements $r$, and the total number of used directions $k$. Our goal is to understand how these hyperparameters impact the convergence to the minimum value.

In our exploration, a multivariate normal distribution is used for the base and the points analyzed, $X \sim N(\mu, \Sigma)$, centered at the origin ($\mu = 0^d$). The covariance matrix $\Sigma$ is a Toeplitz matrix defined as $\Sigma = (\sigma_{ij})_{i=0,...,d-1, j=0,...,d-1}$, with $\sigma_{ij} = 2^{-|i-j|}$.
**Algorithm 4** Massively-parallel refined random search

**function** PROJDEPTHMASSIVE(z, X, k, r, α, κ)

**Initialisation**

$p \leftarrow e_1$

$m \leftarrow \lceil \frac{k}{r} \rceil$

$\epsilon \leftarrow \frac{\pi}{2}$

$D_{min} \leftarrow 1$

**for** $l \leftarrow 1$ **to** $r$ **do**

**Generation:** Generate $m$ random directions ($\epsilon$-around $p$)

$U \leftarrow GenMassive(p, m, \epsilon)$

▷ Parallel directions' generation

**Projection:** Project $X$ and $z$ onto all directions $U$

$P^{(X)} \leftarrow ProjMassive(X, U, \kappa)$

$P^{(z)} \leftarrow ProjMassive(z, U, \kappa)$

**Computation:** Compute $m$ univariate depths within projections

**for** $j \leftarrow 1$ **to** $m$ **do** **(parallel)**

$D \leftarrow D((P^{(z)})_j|(P^{(X)})_{j,\cdot})$

▷ Compute a single univariate depth

**if** $D < D_{min}$ **then**

$D_{min} \leftarrow D$

▷ Update the depth and its corresponding direction

$p \leftarrow U_{j,\cdot}$

$\epsilon \leftarrow \epsilon \cdot \alpha$

**return** $D_{min}$

## 5.1 Convergence and spherical cap shrinkage

First, we study how $\alpha$ affects the final depth value. A data set of 50,000 samples is constructed using a multivariate normal distribution, with the spacial dimension fixed at $d = 50$. From this data set, 50 points were randomly selected and analyzed. For each point, a reference value for the (true) minimal depth was computed using 100 refinements, $\alpha$ set to 0.9, and 3,000,000 directions. The minimal value was selected over three such different runs, each with distinct random directions.

A grid of hyperparameters was used to study the convergence to the minimum value found in the described above reference. Spherical shrinkage values ranged from 0.6 to 0.9, the number of directions ranged from 200 to 90,000, and the number of refinements ranged from 1 to 100. Results for refinements $r \in \{25, 30, 35\}$ are shown in Figure 4, where the mean squared error (MSE) between the minimal value and computed value is measured (see Section ?? of the Supplementary Material for further results.)

Overall, higher $\alpha$ values indicate faster and better convergence to minimal depth. This pattern is observed for all refinement choices and across all analyzed points. Additionally, the volatility between the mean MSE and individual values decreases with higher $\alpha$ values. A slower decrease in the search surface (corresponding to higher values of $\alpha$) may help to avoid getting trapped into local minima and explore a greater part of $\mathbb{S}^{d-1}$, resulting in better RRS convergence. For the subsequent analysis, $\alpha$ is fixed at 0.9, allowing us to focus on the effects of refinements and directions.
<figure><img src="image_5.png" /><figcaption>Figure 4: Convergence study using MSE for $D_P$ with respect to a grid of refinements, directions, and spherical shrinkage values, with $d = 50$. The number of refinements takes value in {25, 30, 35}; the number of directions range from 200 to 90,000 and $\alpha \in \{0.6, 0.7, 0.8, 0.9\}$. The data set size is fixed at 50,000. MSE is computed for 50 random points from $\mathbf{X}$ with respect to a reference value. Light and continuous lines represent the MSE of each point, while bold dashed lines stand for the mean value of all MSEs.</figcaption></figure>

## 5.2 Interdependence of refinements and directions

In this subsection, we explore the impact of refinements and directions on the convergence of the optimization routine. These variables are studied together because they have a more direct correlation with each other. Specifically, $m$, the number of directions used in each refinement, is computed as $\lceil \frac{k}{r} \rceil$.

This investigation is conducted using a grid of values, where the number of directions ranges from 200 to 100,000, the number of refinements ranges from 1 to 175, and the number of dimensions ranges from 5 to 175. The data set $\mathbf{X}$ was constructed from the same multivariate normal distribution as previously, containing 10,000 points. From this data set, 50 random points were selected for analysis. For each point in each dimension, the lowest computed depth was taken as its reference value. All computed depths are compared with this reference to calculate the MSE.

Figure 5 shows RRS convergence with two distinct views. On the left, curves represent the mean number of refinements with a fixed number of directions for points to converge. Only points with $MSE \le 10^{-4}$ (considered as converged) are taken into account. On the right, curves represent the convergence of all analyzed points. On the left side of each curve, not all points have converged yet, and we take the mean number of minimal refinements for converged points. On the right side, all points have $MSE \le 10^{-4}$, and we analyze the combination of refinements and directions needed for convergence when all points have already converged.

The left graph indicates that not all points converge at the same rate; some require fewer refinements and directions, while others may need a combination of higher refinements and directions. The minimal refinement required reaches a plateau, indicating the convergence of all points, which can be seen in both graphs. The right plot shows that a minimal quantity of directions is needed
<figure><img src="image_6.png" /><figcaption>Figure 5: Convergence of the algorithm regarding the interdependence of directions and refinements, with MSE$\le 10^{-4}$ for a multivariate normal distribution. Depth is computed using projection depth for a data set of 10,000 samples. The number of refinements ranges from 10 to 175, the number of dimensions ranges from 1 to 175, and the number of directions ranges from 200 to 100,000. Projection depth is computed for 50 random points from the distribution. Left: Minimal refinement for point-wise convergence. Right: Upper and lower boundaries for the depth convergence of all points with respect to space dimension.</figcaption></figure>

for all points to converge in every dimension. For the number of refinements with fixed directions, the right plot suggests the existence of upper and lower bounds based on the space dimension. A lower refinement boundary is necessary for the searched region to converge, but a high refinement boundary may indicate an insufficient quantity of directions in each refinement for convergence.

We also explore convergence regions for a multivariate exponential distribution (Figure 6). For this distribution, it seems that fewer directions and refinements are needed to reach the lowest depth. Convergence frontiers are closer together on the right side of the figure, while the left part shows earlier convergence even for higher dimensions.

This methodology for hyperparameter selection helps us to understand how the algorithm behaves in higher dimensions and determine the optimal parameter combinations for depth convergence. Our next goal is to verify time gains by comparing both approaches and further developing (3) to explain the observed time behaviors.
<figure><img src="image_7.png" /><figcaption>Figure 6: Convergence of the algorithm regarding the interdependence of directions and refinements with MSE ≤ 10⁻⁴ for an exponential distribution. Depth is computed using projection depth for a data set of 10,000 samples. The number of refinements ranges from 10 to 175; the number of dimensions ranges from 1 to 175, and the number of directions ranges from 200 to 100,000. Projection depth is computed for 50 random points from the distribution. Left: Minimal refinement required for point-wise convergence. Right: Upper and lower boundaries for the depth convergence of all points.</figcaption></figure>

# 6 Runtime analysis in large-scale scenarios

In this section, we investigate the runtime performance in large-scale scenarios and evaluate how the observed temporal behavior aligns with our proposed equations.

## 6.1 Breakdown time analysis

We analyze the significance of each operation highlighted in (2) by measuring their individual runtime contributions. This analysis employs three different approaches to examine the relative impact of each operation on the total runtime.

Figure 7 illustrates the first two approaches, where the time breakdown for the three depth notions ($D_P$, $D_H$, and $D_{AP}$) is presented. The results are averaged over 100 points. This setup allows us to study the impact of varying the number of approximate directions. Overall, univariate depth computation emerges as the rate-determining step. This contrasts with the standard RRS results shown in Figure 1, where projecting the data set consumed a higher percentage of the total time.

Projecting the data set is a significant factor in the overall runtime, accounting for approximately
<figure><img src="image_8.png" /><figcaption>Figure 7: Breakdown time percentage for projection-based depth notions with refinements fixed at 1 and dimensions at 150. The red line represents $D_P$, blue lines $D_H$, and black lines $D_{AP}$. Dashed lines indicate direction generation time, solid lines represent data set projection **X**, and dot-dashed lines denote univariate depth computation. Left: $n$ is fixed at 10,000 while the number of directions is varied. Right: The number of directions $k$ is fixed at 10,000 while the data set size is varied.</figcaption></figure>

one-third of the total time. This demonstrates how our algorithm has decreased the projection time to the point where it is no longer the rate-determining operation. This change is evident in the inversion of the impacts between univariate depth computation and data set projection.

Interestingly, even with $d = 150$, the number of dimensions does not significantly affect the final runtime, as the univariate computation dominates the depth calculation time. Parallelizing these operations results in a final time decomposition that differs from traditional time complexity analysis, where we observe that $Time(O(knd)) < Time(O(kn))$. These findings suggest that the cost of projecting the data, $C_P$, is lower than the cost of univariate depth, $C_{d(1)}$ (since matrix multiplication, which consists of FLOPS, can be more efficiently parallelized on the GPU). Returning to our model, when $d$ is fixed, both data projection and univariate depth computation scale with $n \cdot k$, exhibiting similar performance trends as either $k$ or $n$ increases.

Finally, space dimension is also crucial for time computation and complexity. Figure 8 provides a breakdown of the main steps for $D_P$ calculation. This computation is performed with one refinement, where the number of approximate directions and the data set size are fixed at equal values. As $d$ increases, measures based on projected values become less important, while data set projection regains its role as the most significant operation for depth approximation in higher dimensions. The importance of generating directions also increases with higher dimensions, surpassing the runtime impact of univariate depth computation. However, it remains overshadowed by the data set projection operation.

These results are key for understanding the importance of each factor in Equation 2, and how the algorithm performs in extreme environments. Generating directions becomes more significant
in high-dimensional cases but is overshadowed by other operations. Univariate depth computation emerges as the primary time-consuming operation overall, losing its dominance to data set projection only in very high dimensions (around 10³ variables). Given a consistent computer configuration, the algorithm behavior aligns with the equations proposed in Section 3.

<figure><img src="image_9.png" /><figcaption>Figure 8: Breakdown of operation times for projection depth with respect to space dimension. The number of directions and data set size are fixed according to line color: red for n, k = 10,000, blue for n, k = 15,000, black for n, k = 20,000. Operations are indicated by line style: dashed lines represent direction generation, solid lines represent data set projection, and dot-dashed lines indicate univariate depth computation.</figcaption></figure>

Beyond the breakdown time analysis, the proposed equations describing the speedup factor (3) can also be validated with empirical results.

## 6.2 Empirical speedup analysis

In this part, we compare the runtime of both approaches using the same hyperparameters, and we verify how well the empirical results match the theoretical speedup factor described in (3).

Different notions of depth will exhibit varying speedup factors, as they are influenced by $\mathcal{D}$, as formalized in (1) and (2). The variation in speedup for the three notions is illustrated in Figure 9. Overall, the speedup factor increases with the number of directions until it reaches a plateau, which depends on the dimension $d$ and the depth notion. This result is aligned with (3). Note that the speedup scale in the graphs is adjusted for better visualization of the phenomenon.

Using previous results, we can develop $S_p$. Given that $C_{rv}$ and $C_c$ are negligible compared to other operations and considering that $\mathcal{D} \sim m \cdot n$, the speedup factor becomes:

$$S_p = \frac{m \cdot n \cdot (C_p \cdot d + C_{d(1)})}{\lceil \frac{m \cdot n}{g} \rceil \cdot \lambda \cdot \left( C_p \cdot \lceil \frac{d}{d_{max(m,n)}} \rceil + C_{d(1)} \right)}$$
<figure><img src="image_10.png" /><figcaption>Figure 9: Speedup factor for projection-based depth notions as a function of approximating directions and space dimension. From left to right: halfspace, projection, and asymmetric projection depth. Ordinate scale are set to better accommodate computed values providing better visualization. Colors represent the space dimension, orange $d = 10$, purple $d = 25$, red $d = 50$, blue $d = 100$, and black $d = 150$.</figcaption></figure>

From this equation, we observe two key behaviors. First, the increase in $S_p$ occurs for smaller values of $m \cdot n$. As this product increases, the function shows an amortized convergence to the maximum speedup value, with small fluctuations perceived as a plateau value. The latter depends only on the dimension $d$. To better illustrate this plateau, we consider $m \cdot n$ as a multiple of $g$, i.e., $m \cdot n = w \cdot g$, where $w$ is a factor that cancels out, resulting in (4):

$$S_p = \frac{g}{\lambda} \cdot \frac{C_p \cdot d + C_{d(1)}}{C_p \cdot \lceil \frac{d}{d_{max(m,n)}} \rceil + C_{d(1)}} \quad (4)$$

As stated, this behavior depends on the product $m \cdot n$ and the dimension, which is supported by the results shown in the left graph of Figure 10, exploring $D_p$ computation. With $d$ fixed at 150, approximately the same $S_p$ plateau value is reached by increasing either $n$ or $m$. The rate at which the plateau is reached depends on $n \cdot m$, so the curves do not overlap with lower numbers of directions. Time was measured as the mean time over 100 points.

Continuing the exploration of projection notions, the impact of $d$ on $S_p$ is examined in the right graph of Figure 10. Time is computed with a fixed data set of 10.000 samples, 10.000 directions, and one refinement. The space dimension is increased from 5 to 50.000 in smaller increments than in previous explorations. A similar pattern is observed, showing an increase in speedup up to a point, followed by a fluctuating plateau where the speedup value oscillates around this value with an amortized convergence. This behavior is similar to how $\lceil \frac{d}{d_{max(m,n)}} \rceil$ behaves like $\lceil \frac{m \cdot n}{g} \rceil$ in (4).

Compared to previous state-of-the-art approaches, we achieved a time reduction of three orders of magnitude, with the importance of data set projection and univariate depth inverted. In the next section, we shift our focus from comparing $T_C$ with $T_G$ to analyzing and visualizing the absolute runtime values observed with different parameters.
<figure><img src="image_11.png" /><figcaption>Figure 10: Speedup factor analysis. Left: Projection depth with increasing approximate directions and data set size, with dimension fixed at 150; each color represents a different data set size. Right: Speedup factor for different projection-based notions as the dimension increases. Both directions and data set samples are fixed at 10,000. In both cases, computations are performed with one refinement.</figcaption></figure>

## 6.3 Numerical time analysis

In these experiments, computational time is measured for a grid of parameters including the number of directions, dimensions, data set sizes, and refinements. These measurements are conducted using $D_P$, since all depth notions exhibit similar time behavior. The results (in seconds) are presented in Table 1, where $D_P$ is evaluated with a data set of 10,000 samples and a single refinement.

The time appears to grow proportionally with the number of directions after a threshold of approximately 5,000 directions, aligning with the observations in Figure 9. Interestingly, the dimension seems to have little to no impact on the final runtime, except when the number of approximate directions is high. The last row deviates from this proportional pattern, exhibiting a higher runtime than expected.

<table><caption>Table 1: Computational time (in seconds) to calculate depth for a data set of 10,000 points using one refinement.</caption><thead><tr><th rowspan="2">Directions</th><th colspan="8">Space dimension</th></tr><tr><th>5</th><th>10</th><th>25</th><th>50</th><th>75</th><th>100</th><th>125</th><th>150</th></tr></thead><tbody><tr><td>1,000</td><td>0.01</td><td>0.01</td><td>0.01</td><td>0.01</td><td>0.01</td><td>0.01</td><td>0.01</td><td>0.01</td></tr><tr><td>2,500</td><td>0.01</td><td>0.01</td><td>0.01</td><td>0.01</td><td>0.01</td><td>0.01</td><td>0.01</td><td>0.01</td></tr><tr><td>5,000</td><td>0.02</td><td>0.02</td><td>0.02</td><td>0.02</td><td>0.02</td><td>0.02</td><td>0.02</td><td>0.02</td></tr><tr><td>10,000</td><td>0.03</td><td>0.03</td><td>0.03</td><td>0.03</td><td>0.03</td><td>0.03</td><td>0.03</td><td>0.03</td></tr><tr><td>25,000</td><td>0.06</td><td>0.07</td><td>0.07</td><td>0.07</td><td>0.07</td><td>0.07</td><td>0.07</td><td>0.08</td></tr><tr><td>50,000</td><td>0.13</td><td>0.14</td><td>0.14</td><td>0.14</td><td>0.14</td><td>0.15</td><td>0.15</td><td>0.16</td></tr><tr><td>100,000</td><td>0.26</td><td>0.28</td><td>0.28</td><td>0.30</td><td>0.28</td><td>0.28</td><td>0.30</td><td>0.30</td></tr><tr><td>250,000</td><td>2.49</td><td>2.49</td><td>2.52</td><td>2.59</td><td>2.60</td><td>2.60</td><td>2.64</td><td>2.66</td></tr></tbody></table>
This time increase can be partially explained by referring back to (4). After reaching the maximum speedup, the final change in value $S_p$ is $\lceil \frac{d}{d_{max(m,n)}} \rceil$, indicating memory saturation. However, the most significant increment is in the univariate depth computation due to its memory-bound operations, flow divergence, and low-level conditional steps [Fung et al., 2007]. Data set projection decreases in importance for the final runtime, dropping from approximately 30% to around 16%, while univariate depth computation increases from 70% to 83% of the final time.

Refinements can play a crucial role in computational time, as they not only introduce the optimization character into the algorithm, but can also reduce the memory used in each refinement, mitigating memory overload and potential efficiency loss. Table 2 shows the computational time for $D_P$ with a 10,000-point data set using 250,000 directions, with the number of refinements ranging from 1 to 100. Significant time reductions are observed with the inclusion of more refinements, which helps alleviate memory overload and brings the runtime closer to the expected values from Table 1. Beyond the initial decrease, further increases in refinements do not significantly affect the runtime.

<table><caption>Table 2: Computational time (in seconds) to calculate $D_P$ with a data set of 10,000 samples and 250,000 approximate directions.</caption><thead><tr><th rowspan="2">Refinements</th><th colspan="8">Space dimension</th></tr><tr><th>5</th><th>10</th><th>25</th><th>50</th><th>75</th><th>100</th><th>125</th><th>150</th></tr></thead><tbody><tr><td>1</td><td>2.49</td><td>2.49</td><td>2.52</td><td>2.59</td><td>2.60</td><td>2.60</td><td>2.64</td><td>2.66</td></tr><tr><td>10</td><td>0.62</td><td>0.67</td><td>0.67</td><td>0.69</td><td>0.69</td><td>0.69</td><td>0.74</td><td>0.76</td></tr><tr><td>25</td><td>0.61</td><td>0.65</td><td>0.65</td><td>0.68</td><td>0.69</td><td>0.69</td><td>0.70</td><td>0.71</td></tr><tr><td>50</td><td>0.61</td><td>0.65</td><td>0.65</td><td>0.67</td><td>0.68</td><td>0.68</td><td>0.70</td><td>0.73</td></tr><tr><td>100</td><td>0.66</td><td>0.67</td><td>0.67</td><td>0.68</td><td>0.70</td><td>0.69</td><td>0.71</td><td>0.72</td></tr></tbody></table>

The size of the data set is also a key factor in computational time. Table 3 focuses on varying the data set size while keeping refinements fixed at 1 and the number of directions fixed at 10,000. Our results indicate a similar behavior between increasing the number of directions and increasing the data set size, as expected from (2).

Overall, the numerical time analysis reveals that the number of directions and data set size have a more significant impact on computational time than the space dimension. Time does not increase significantly with higher dimensions.

Beyond computational time, we also investigate data ranking and ordering with respect to the distribution center to verify the precision of data ordering.
<table><caption>Table 3: Computational time (in seconds) to calculate $D_P$ for a data set with 20,000 directions and one refinement.</caption><thead><tr><th></th><th colspan="8">Space dimension</th></tr><tr><th>Data set size</th><th>5</th><th>10</th><th>25</th><th>50</th><th>75</th><th>100</th><th>125</th><th>150</th></tr></thead><tbody><tr><th>1,000</th><td>0.01</td><td>0.01</td><td>0.01</td><td>0.01</td><td>0.01</td><td>0.01</td><td>0.01</td><td>0.01</td></tr><tr><th>2,500</th><td>0.01</td><td>0.01</td><td>0.01</td><td>0.01</td><td>0.01</td><td>0.01</td><td>0.01</td><td>0.01</td></tr><tr><th>5,000</th><td>0.01</td><td>0.02</td><td>0.01</td><td>0.01</td><td>0.02</td><td>0.02</td><td>0.02</td><td>0.02</td></tr><tr><th>10,000</th><td>0.03</td><td>0.03</td><td>0.03</td><td>0.03</td><td>0.03</td><td>0.03</td><td>0.03</td><td>0.03</td></tr><tr><th>25,000</th><td>0.06</td><td>0.06</td><td>0.07</td><td>0.07</td><td>0.07</td><td>0.07</td><td>0.08</td><td>0.08</td></tr><tr><th>50,000</th><td>0.14</td><td>0.14</td><td>0.15</td><td>0.14</td><td>0.15</td><td>0.15</td><td>0.17</td><td>0.16</td></tr><tr><th>100,000</th><td>0.29</td><td>0.28</td><td>0.30</td><td>0.29</td><td>0.29</td><td>0.33</td><td>0.31</td><td>0.33</td></tr></tbody></table>

# 7 Center-outward ordering precision and depth benchmark

The natural ordering provided by depth notions can be evaluated in higher dimensions using ranking correlation techniques. A data set with 100,000 samples is constructed using a normal distribution and a Toeplitz matrix, as previously formalized. Since the distribution is known, its probability density function (pdf) can serve as a reference to ensure the correct ordering of the data.

Three depth notions are used for this ranking assessment: $D_P$, $D_{AP}$, and Mahalanobis depth ($D_M$). The Mahalanobis depth is based on the Mahalanobis distance [Mahalanobis, 2018] and is computed as:

$$D_M(\mathbf{z}|\mathbf{X}) = (1 + (\mathbf{z} - \boldsymbol{\mu})'(\boldsymbol{\Sigma}^{-1})(\mathbf{z} - \boldsymbol{\mu}))^{-1},$$

where $\boldsymbol{\mu}$ is the mean of $\mathbf{X}$ and $\boldsymbol{\Sigma}$ is its variance-covariance matrix.

Table 4 presents the ranking correlation of 5,000 randomly chosen points from the data set. Each point was analyzed using $D_P$ and $D_{AP}$, both using 40 refinements, a total of 100,000 directions, and $\alpha$ set to 0.9 based on the previous analysis. $D_M$ was computed using two different estimators for $\mu$ and $\Sigma$:

* $D_{M(L)}$: Uses the maximum likelihood estimator (MLE) for $\Sigma$ and an unbiased estimator for the location, $\mu = \frac{1}{n} \sum_{i=1}^{n} X_i$.

* $D_{M(R)}$: Employs a robust estimation method, the minimum covariance determinant (MCD) with a factor of 0.5, for both the location $\mu$ and the scatter $\Sigma$ [Rousseeuw and Leroy, 2003, Rousseeuw and Driessen, 1999].

The study uses two robust ranking correlation techniques, namely Kendall's $\tau$ and Spearman's $\rho$ correlation [Spearman, 1987, Kendall, 1938, Croux and Dehon, 2010]. Both $D_P$ and $D_{AP}$ show very strong correlations with the real pdf ranking across all studied dimensions, although the correlation
values (as expected) decrease with higher space dimension. Overall, $D_M$ exhibits high correlation values for both estimation approaches, since both techniques are well-suited for normal distributions.

<table><caption>Table 4: Ranking correlation for depth ordering with multivariate Gaussian distribution.</caption><thead><tr><th rowspan="2">Correlation</th><th colspan="6">Dimensions</th></tr><tr><th>5</th><th>10</th><th>25</th><th>50</th><th>100</th><th>150</th></tr></thead><tbody><tr><th colspan="7">Spearman ρ</th></tr><tr><td>PDF x D<sub>P</sub></td><td>.9999</td><td>.9998</td><td>.9995</td><td>.9987</td><td>.9963</td><td>.9912</td></tr><tr><td>PDF x D<sub>AP</sub></td><td>.9987</td><td>.9969</td><td>.9903</td><td>.9812</td><td>.9549</td><td>.9140</td></tr><tr><td>PDF x D<sub>M(L)</sub></td><td>.9999</td><td>.9998</td><td>.9998</td><td>.9996</td><td>.9993</td><td>.9990</td></tr><tr><td>PDF x D<sub>M(R)</sub></td><td>.9999</td><td>.9998</td><td>.9998</td><td>.9996</td><td>.9994</td><td>.9990</td></tr><tr><td>D<sub>P</sub> x D<sub>AP</sub></td><td>.9987</td><td>.9971</td><td>.9907</td><td>.9810</td><td>.9536</td><td>.9115</td></tr><tr><th colspan="7">Kendall's τ</th></tr><tr><td>PDF x D<sub>P</sub></td><td>.9929</td><td>.9891</td><td>.9821</td><td>.9687</td><td>.9483</td><td>.9216</td></tr><tr><td>PDF x D<sub>AP</sub></td><td>.9697</td><td>.9530</td><td>.9172</td><td>.8820</td><td>.8157</td><td>.7468</td></tr><tr><td>PDF x D<sub>M(L)</sub></td><td>.9945</td><td>.9921</td><td>.9892</td><td>.9848</td><td>.9793</td><td>.9741</td></tr><tr><td>PDF x D<sub>M(R)</sub></td><td>.9962</td><td>.9932</td><td>.9902</td><td>.9851</td><td>.9802</td><td>.9745</td></tr><tr><td>D<sub>P</sub> x D<sub>AP</sub></td><td>.9704</td><td>.9545</td><td>.9187</td><td>.8813</td><td>.8123</td><td>.7431</td></tr></tbody></table>

Using the same number of points, directions, $\Sigma$, $\mu$, and $d = 50$, we also ranked four Student's $t$ distributions with degrees of freedom $\nu = 1, 2, 5, 10$. The results are displayed in Table 5. For smaller degrees of freedom, $D_P$ and $D_{AP}$ exhibit higher $\rho$ and $\tau$ values than $D_{M(L)}$. Overall, the correlation results are higher with a Student's $t$ than with an underlying Gaussian distribution.

<table><caption>Table 5: Ranking correlation for depth ordering with Student's $t$ distribution.</caption><thead><tr><th rowspan="2">Correlation</th><th colspan="4">ν</th><th rowspan="2">Gaussian</th></tr><tr><th>1</th><th>2</th><th>5</th><th>10</th></tr></thead><tbody><tr><th colspan="6">Spearman ρ</th></tr><tr><td>PDF x D<sub>P</sub></td><td>.9999</td><td>.9998</td><td>.9995</td><td>.9991</td><td>.9987</td></tr><tr><td>PDF x D<sub>AP</sub></td><td>.9996</td><td>.9992</td><td>.9982</td><td>.9967</td><td>.9812</td></tr><tr><td>PDF x D<sub>M(R)</sub></td><td>.9999</td><td>.9999</td><td>.9999</td><td>.9999</td><td>.9996</td></tr><tr><td>PDF x D<sub>M(L)</sub></td><td>.9936</td><td>.9978</td><td>.9999</td><td>.9999</td><td>.9996</td></tr><tr><td>D<sub>P</sub> x D<sub>AP</sub></td><td>.9995</td><td>.9990</td><td>.9978</td><td>.9960</td><td>.9810</td></tr><tr><th colspan="6">Kendall's τ</th></tr><tr><td>PDF x D<sub>P</sub></td><td>.9907</td><td>.9871</td><td>.9804</td><td>.9741</td><td>.9699</td></tr><tr><td>PDF x D<sub>AP</sub></td><td>.9836</td><td>.9771</td><td>.9638</td><td>.9514</td><td>.8820</td></tr><tr><td>PDF x D<sub>M(R)</sub></td><td>.9973</td><td>.9964</td><td>.9944</td><td>.9928</td><td>.9848</td></tr><tr><td>PDF x D<sub>M(L)</sub></td><td>.9332</td><td>.9608</td><td>.9936</td><td>.9934</td><td>.9851</td></tr><tr><td>D<sub>P</sub> x D<sub>AP</sub></td><td>.9817</td><td>.9742</td><td>.9597</td><td>.9461</td><td>.8813</td></tr></tbody></table>
Ranking correlation is performed to assess how well the collective centrality of points is measured with respect to their underlying empirical distribution and how these values can be properly ranked. A key factor for optimized results is the choice of hyperparameters, which is based on our exploration method. $D_P$ and $D_{AP}$ can be further fine-tuned to achieve more precise results by reducing the MSE value of the depth computation.

# 8 Conclusion

This article introduced a novel methodology for the massive parallelization of projection-based depths, addressing the longstanding computational challenges associated with data depth in high-dimensional spaces. By leveraging the Refined Random Search algorithm and modern GPU architectures, our approach achieves unprecedented speedup (up to 7,000 times faster than in the state-of-the-art work by Dyckerhoff et al. [2021], and this on a laptop graphics card) while maintaining the precision of depth calculations. This advancement enables the practical application of data depth in large-scale and high-dimensional data sets, opening new avenues for robust statistical analysis in fields such as anomaly detection, classification, and dimension reduction.

Our framework decomposes the computation of projection-based depths into parallelizable tasks, focusing on the three main operations: direction generation, data projection, and univariate depth computation. By optimizing these operations for GPU execution, we significantly reduce runtime without compromising accuracy. The empirical results, validated through synthetic data, demonstrate the effectiveness of our approach in handling extreme data scenarios, achieving both efficiency and precision.

Furthermore, our study provides insights into hyperparameter fine-tuning, such as the selection of spherical cap shrinkage, the number of refinements, and directions, to optimize convergence and runtime. These findings ensure the robustness and adaptability of our method across different data sets and space dimensions. The center-outward ordering precision, evaluated through ranking correlation techniques, confirms the reliability of our approach in high-dimensional data analysis.

The implications of this work are far-reaching. By overcoming the computational bottlenecks of data depth, our methodology paves the way for broader adoption (e.g., also in multivariate functional spaces Claeskens et al. [2014] or the space of curves Lafaye De Micheaux et al. [2021]) of this powerful statistical tool in large-scale applications. Future research could explore the extension of this framework to other depth notions, further optimization of GPU utilization, and the application of our method to real-world data sets. Additionally, investigating the integration of our approach with other statistical techniques, such as robust estimation and machine learning algorithms, could yield even more powerful analytical tools.
Finally, we have made the RRS algorithm (and other depth functions) available in the **Python-library data-depth**, providing researchers and practitioners with the ready-to-use tools to implement and to build upon our work.

## Acknowledgments

This research has been conducted within the Research Chair "Digital Finance" under the aegis of the Risk Foundation, a joint initiative by Cartes Bancaires CB, Telecom Paris and University of Paris 2 Panthéon-Assas. The authors further gratefully acknowledge the support of the Young Researcher Grant of the French National Agency for Research (ANR JCJC 2021) in category Artificial Intelligence registered under the number ANR-21-CE23-0029-01 and the support of the CIFRE grant number 2021/1739. We would like to thank seminar participants at the 2024 International Conference on Statistics and Data Science (ICSDS). We also thank Samuel Willy and Hafid Chakir for their support.

## Supplementary material

This article contains supplementary materials, which provide:

* Additional information about simple random search and refined random search.

* Further results and illustrations regarding depth convergence and depth outcomes.

* A more detailed explanation of the parallel construction of threads and their key constituents.

## References

G. M. Amdahl. Validity of the single processor approach to achieving large scale computing capabilities, reprinted from the afips conference proceedings, vol. 30 (atlantic city, nj, apr. 18–20), afips press, reston, va., 1967, pp. 483–485, when dr. amdahl was at international business machines corporation, sunnyvale, california. *IEEE Solid-State Circuits Society Newsletter*, 12(3):19–20, 2007.

I. Cascos. The expected convex hull trimmed regions of a sample. *Computational Statistics*, 22: 557–569, 2007.

G. Claeskens, M. Hubert, L. Slaets, and K. Vakili. Multivariate functional halfspace depth. *Journal of the American Statistical Association*, 109(505):411–423, 2014.
C. Croux and C. Dehon. Influence functions of the spearman and kendall correlation measures. *Statistical methods & applications*, 19:497-515, 2010.

R. Dyckerhoff. Data depths satisfying the projection property. *Allgemeines Statistisches Archiv*, 88: 163-190, 2004.

R. Dyckerhoff and K. Mosler. Weighted-mean trimming of multivariate data. *Journal of Multivariate Analysis*, 102(3):405-421, 2011.

R. Dyckerhoff and P. Mozharovskyi. Exact computation of the halfspace depth. *Computational Statistics & Data Analysis*, 98:19-30, 2016.

R. Dyckerhoff, P. Mozharovskyi, and S. Nagy. Approximate computation of projection depths. *Computational Statistics & Data Analysis*, 157:107166, 2021.

V. Fojtík, P. Laketa, P. Mozharovskyi, and S. Nagy. On exact computation of tukey depth central regions. *Journal of Computational and Graphical Statistics*, pages 1-26, 2023.

W. W. Fung, I. Sham, G. Yuan, and T. M. Aamodt. Dynamic warp formation and scheduling for efficient gpu control flow. In *40th Annual IEEE/ACM International Symposium on Microarchitecture (MICRO 2007)*, pages 407-420. IEEE, 2007.

J. L. Gustafson. Reevaluating amdahl's law. *Communications of the ACM*, 31(5):532-533, 1988.

M. G. Kendall. A new measure of rank correlation. *Biometrika*, 30(1-2):81-93, 06 1938. ISSN 0006-3444. doi: 10.1093/biomet/30.1-2.81.

G. Koshevoy and K. Mosler. Zonoid trimming for multivariate distributions. *The Annals of Statistics*, 25(5):1998-2017, 1997.

P. Lafaye De Micheaux, P. Mozharovskyi, and M. Vimond. Depth for curve data and applications. *Journal of the American Statistical Association*, 116(536):1881-1897, 2021.

P. C. Mahalanobis. On the generalized distance in statistics. *Sankhyā: The Indian Journal of Statistics, Series A (2008-)*, 80:S1-S7, 2018.

A. Marowka. Extending amdahl's law for heterogeneous computing. In *2012 IEEE 10th International Symposium on Parallel and Distributed Processing with Applications*, pages 309-316. IEEE, 2012.

K. Mosler and P. Mozharovskyi. Choosing among notions of multivariate depth statistics. *Statistical Science*, 37(3):348-368, 2022.
C. A. Navarro, N. Hitschfeld-Kahler, and L. Mateu. A survey on parallel computing and its applications in data-parallel problems using gpu architectures. *Communications in Computational Physics*, 15(2):285-329, 2014.

P. J. Rousseeuw and K. V. Driessen. A fast algorithm for the minimum covariance determinant estimator. *Technometrics*, 41(3):212-223, 1999.

P. J. Rousseeuw and A. M. Leroy. *Robust regression and outlier detection*. John wiley & sons, 2003.

C. Spearman. The proof and measurement of association between two things. *The American Journal of Psychology*, 100(3/4):441-471, 1987. ISSN 00029556.

J. W. Tukey. Mathematics and the picturing of data. In *Proceedings of the International Congress of Mathematicians, Vancouver, 1975*, volume 2, pages 523-531, 1975.

Y. Zuo and S. Lai. Exact computation of bivariate projection depth and the stahel-donoho estimator. *Computational Statistics & Data Analysis*, 55(3):1173-1179, 2011.

Y. Zuo and R. Serfling. General notions of statistical depth function. *Annals of statistics*, pages 461-482, 2000.

J. Nickolls and W. J. Dally. The gpu computing era. *IEEE micro*, 30(2):56-69, 2010.

S. Ryoo, C. I. Rodrigues, S. S. Baghsorkhi, S. S. Stone, D. B. Kirk, and W.-m. W. Hwu. Optimization principles and application performance evaluation of a multithreaded gpu using cuda. In *Proceedings of the 13th ACM SIGPLAN Symposium on Principles and practice of parallel programming*, pages 73-82, 2008.
# Supplementary material to the article
## "Massive parallelization of projection-based depths"

Leonardo Leone, Pavlo Mozharovskyi, David Bounie

This supplementary material provides additional information about simple random search and refined random search in Section A. Section B offers further results and illustrations regarding depth convergence and depth outcomes. Finally, a more detailed explanation of the parallel construction of threads and their key constituents is presented in Section C.

## A Simple random search and refined random search

**Simple random search.** The first structure for the proposed method is developed using RS. Depth is computed with respect to a set of randomly generated directions using an algorithm called (Simple) Random Search (RS). In RS, directions are randomly generated on the unit sphere ($\mathbb{S}^{d-1} = \{x \in \mathbb{R}^d : \|x\| = 1\}$) following a normal distribution, as shown in Algorithm S1. The data set $\mathbf{X}$ is then projected, and the respective depth is computed. The smallest depth value for $k$ directions is kept as the approximated depth.

<table><caption>**Algorithm S1** Create a random direction on the unit sphere from a normal distribution</caption><tr><td><b>function</b> RANDOMSPHERE(d)</td></tr><tr><td><b>u</b> &larr; <b>0</b><sub>d</sub></td></tr><tr><td><b>s</b> &larr; 0</td></tr><tr><td><b>for</b> i &larr; 1 to d do</td></tr><tr><td><b>u</b><sub>i</sub> &larr; N(0,1)</td></tr><tr><td><b>s</b> &larr; <b>s</b> + <b>u</b><sub>i</sub><sup>2</sup></td></tr><tr><td><b>for</b> i &larr; 1 to d do</td></tr><tr><td><b>u</b><sub>i</sub> &larr; <sup><b>u</b><sub>i</sub></sup>&frasl;<sub>&radic;<b>s</b></sub></td></tr><tr><td><b>return</b> <b>u</b><sub>i</sub></td></tr></table>

Algorithm S2 presents the respective pseudocode for a generic depth notion computation using RS.

**Refined random search.** Refined random search (RRS) is an enhanced variation of RS that leverages refinements to concentrate the search direction within specific regions. This method involves partitioning the total number of search directions $k$ into $r$ distinct refinements, where each refinement is performed using $m = \lfloor \frac{k}{r} \rfloor$ directions. The underlying concept of RRS is to perform
<caption>**Algorithm S2** Computing depth for simple random search</caption>

```python
function SIMPLERANDOMSEARCH(z, X, d, k)
    D_min ← 1
    for i ← 1 to k do
        u ← randomSphere(d)
        D ← D(zu$^{T}$|Xu$^{T}$)
        D_min ← min(D, D_min)
    return D_min
```

depth computation through a series of iterative steps. This process employs a set of procedures to identify and select a region of interest, subsequently generating search directions within the relevant subspaces. By focusing the search effort in targeted regions, RRS aims to improve the efficiency and effectiveness of the optimization process.

The first step is performed using RS, with a focus on a hemispherical cap. It is important to note that for this purpose, the depth function satisfies the property $D(zu^T|Xu^T) = D(-zu^T| -Xu^T)$.

At each refinement, the depth is computed with respect to a set of directions $U = (u_1, ..., u_m)$, which are generated around a central point referred to as the pole. Among the $m$ mapped depth values, the lowest value is selected as the current minimal depth. The corresponding direction $u$ is associated with this minimal depth is then designated as the new pole for subsequent iterations.

This selected direction serves as the center of the region in which the subsequent mapping is conducted. The size of the explored region is decreased by a factor $\alpha$, a process referred to as spherical cap shrinkage. A new set of directions is then generated within this smaller spherical cap, and the entire process is repeated $r$ times. During each refinement, a new pole is identified based on the direction corresponding to the lowest depth value.

This approach can be decomposed into granular logical functions, each serving specific purposes, in alignment with our simple function approach. The primary steps for RRS are described in Algorithm S3. It is important to note that certain arguments have been modified in comparison to the random search (RS) presented in Algorithm S2. In this context, the number of directions $k$ denotes the total number of directions considered across all refinements combined, where $m$ is the number of directions in each of the $r$ refinements. Spherical shrinkage refers to the factor by which the size of the spherical cap is reduced ($\epsilon$), initialized as $\frac{\pi}{2}$.

Generating random directions within the search region is not as straightforward as in RS direction generation. To account for region constraints, Algorithm S1 is modified with a series of reflections and transformations, as outlined in Algorithm S4, to map the direction within the desired region.
<caption>**Algorithm S3** Computing depth with Refined Random Search</caption>

**function** REFINEDRANDOMSEARCH($z, X, k, r, \alpha, d$)
    $p \leftarrow e_1$
    $m \leftarrow \lceil \frac{k}{r} \rceil$
    $\epsilon \leftarrow \frac{\pi}{2}$
    $D_{min} \leftarrow 1$
    **for** $l \leftarrow 1$ **to** $r$ **do**
        **for** $i \leftarrow 1$ **to** $m$ **do**
            $u \leftarrow$ randomSpherePole($p, \epsilon$)
            $D \leftarrow D(zu^\top | Xu^\top)$
            **if** $D < D_{min}$ **then**
                $D_{min} \leftarrow D$
                $u_j \leftarrow u$
            $p \leftarrow u_j$
            $\epsilon \leftarrow \epsilon \cdot \alpha$
        **return** $D_{min}$

<caption>**Algorithm S4** Creating a random direction on the unit sphere within the search region</caption>

**function** RANDOMSPHEREPOLE($p, \epsilon, d$)
    $u \leftarrow 0_d$
    $u_1 \leftarrow \cos(U([0, 1]) \cdot \epsilon)$
    $u_{2:d} \leftarrow \sqrt{1 - u_1^2} \cdot \text{randomSphere}(d-1)$
    $\lambda \leftarrow (pu^\top - u_1) / (1 - p_1)$
    $u_1 = u_1 + \lambda$
    **for** $i \leftarrow 0$ **to** $d$ **do**
        $u_i = u_i - \lambda$
    **return** $u$

## B Depth convergence and illustration

The convergence of depth to its minimum value is analyzed using a multivariate normal distribution, as detailed in Section 5 of the paper. Figure S1 illustrates this convergence for 50 sample points, examining the relationship between spherical cap shrinkage ($\alpha$) - represented by different colored lines - and the number of directions. The space dimension is fixed at $d=50$. Each subplot corresponds to a different number of refinements, ranging from 20 to 100. The results consistently demonstrate that higher values of ($\alpha$) lead to faster convergence, irrespective of the number of refinements employed.

The relative squared error of computed projection depth values is presented in Figure S2, comparing them against minimum computed values. The analysis examines various combinations of space dimensions, refinement quantities, and direction counts. The findings reveal that the computation exhibits greater sensitivity to refinement levels when fewer approximate directions are used.
<figure><img src="image_12.png" /><figcaption>Figure S1: Convergence study using MSE for $D_P$ with respect to a grid of refinements, directions, and spherical shrinkage values, with $d = 50$. The number of refinements takes value $r \in \{20, 40, 45, 50, 55, 60, 65, 75, 100\}$; the number of directions range from 200 to 90,000 and $\alpha \in \{0.6, 0.7, 0.8, 0.9\}$. The data set size is fixed at 50,000. MSE is computed for 50 random points from $\mathbf{X}$ with respect to a reference value.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure S2: Relative squared error depth box-plot for projection depth. Computation takes into consideration a data set with 10,000 points generated from a multivariate normal distribution centered in $0_d$ and a Toeplitz matrix for the covariance, as described in Section 5 of the main text. From this data set, 50 points are selected for depth analysis.</figcaption></figure>
C Graphical processing unit parallel programming

The objective of parallel computation is to expedite independent calculations by harnessing the power of multiple CPU processors or, as utilized in this article, multi-core graphical processing units (GPUs). Transitioning from a CPU-based to a GPU-based program is not a trivial task; several critical steps are required for an application to operate effectively in such a multi-core and multi-threaded environment. Key considerations include thread and device synchronization, memory allocation, cached memory access, and flow divergence. Each of these aspects must be meticulously addressed to ensure optimal performance and efficiency in a parallel computing framework.

Calculations performed on the graphics card utilize a multithreaded approach, wherein computations and functions are decomposed into smaller instructions known as threads. Each thread is constructed based on a specific task assigned to the computer, referred to as the (kernel code). Threads are organized and grouped together, typically in sets of 32, forming a unit known as a warp. A collection of warps constitutes a threadblock. A collection of threadblocks is called a grid, and threadblocks within the same grid are allocated and executed in a Streaming Multiprocessor (SM).

This architecture is constructed and specified using an execution configuration that delineates the number of blocks per grid and threads per block. At a lower level, the kernel code incorporates predefined variables, namely `blockId` and `threadId`. These variables are instrumental in efficiently monitoring the computations being executed in various sections of the GPU, thereby mitigating potential asynchronicity issues [Ryoo et al., 2008, Nickolls and Dally, 2010].

In addition to the extra considerations that require computational time, it is important to note that GPU cores operate at a lower frequency, measured by the number of cycles executed in a given time period, compared to CPU cores. Consequently, individual computations tend to be slower on graphics units due to factors such as memory allocation and lower core performance. However, graphics units compensate for this by leveraging a significantly higher number of cores, typically numbering in the thousands. In contrast, CPUs usually feature dozens of cores. This architectural difference allows GPUs to excel in parallel processing tasks, despite the slower performance of individual cores.

Moreover, the necessity to break down operations into simpler tasks mandates the creation of threads with straightforward commands. These parallel computations are particularly well-suited for floating-point operations, such as matrix multiplication. Given these advantages and limitations, the creation and execution of threads must be meticulously coordinated, especially for large-scale applications. In extensive computational routines, leveraging parallel multithreading can result
in substantial time savings, effectively addressing one of the most significant challenges in depth computation.

Depth computation can significantly benefit from parallel approaches not only in data set projection but also in other operations. Generating random directions and computing univariate depth, although on a smaller scale, can also leverage these methods effectively. A simplified overview of the complete Refined Random Search (RRS) process is presented in the pseudocode S5. More intricate steps are involved in random generation and univariate depth operations. These include median computation for univariate $D_P$ calculation and conditional comparisons for $D_H$. These more extensive computations are particularly susceptible to flow divergence, necessitating special attention to ensure optimal performance.

<caption>**Algorithm S5** Simple view of parallel computation of refined random search</caption>

function GENERALPARALLELRRS($X, z, r, k, \alpha, d$)
    **Initialisation**
    $D_m \leftarrow 1$
    $\epsilon \leftarrow \frac{\pi}{2}$
    $m \leftarrow \lceil \frac{k}{m} \rceil$
    $U \leftarrow 0_{m \times d}$
    $p \leftarrow e_1$
    **for** $l \leftarrow 1$ **to** $r$ **do**
        **Step 1: Generation**
        $U \leftarrow DirMassive(p, m, \epsilon)$
        **Step 2: Projection**
        $P^{(X)} \leftarrow ProjMassive(X, U, \kappa)$
        $P^{(z)} \leftarrow ProjMassive(z, U, \kappa)$
        **Step 3: Computation**
        $D \leftarrow UniDepthMassive(P^{(z)}, P^{(z)}, \text{notion})$
        $D_c \leftarrow \min(D)$
        $\epsilon \leftarrow \epsilon \cdot \alpha$
        **if** $D_c < D_m$ **then**
            $D_m \leftarrow D_c$
            $p \leftarrow U_c$
            $\triangleright U_c$ is the direction corresponding to $D_c$
    **return** $D_m$