Amir Abboud
email: amir.abboud@gmail.com

Vincent Cohen-Addad

Hussein Houdrouge
email: hussein.houdrouge@polytechnique.edu

Subquadratic High-Dimensional Hierarchical Clustering

We consider the widely-used average-linkage, single-linkage, and Ward's methods for computing hierarchical clusterings of high-dimensional Euclidean inputs. It is easy to show that there is no efficient implementation of these algorithms in high dimensional Euclidean space since it implicitly requires to solve the closest pair problem, a notoriously difficult problem. However, how fast can these algorithms be implemented if we allow approximation? More precisely: these algorithms successively merge the clusters that are at closest average (for average-linkage), minimum distance (for single-linkage), or inducing the least sum-of-square error (for Ward's). We ask whether one could obtain a significant running-time improvement if the algorithm can merge γ-approximate closest clusters (namely, clusters that are at distance (average, minimum, or sum-of-square error) at most γ times the distance of the closest clusters). We show that one can indeed take advantage of the relaxation and compute the approximate hierarchical clustering tree using r Opnq γ-approximate nearest neighbor queries. This leads to an algorithm running in time r

Opndq `n1`Op1{γq for d-dimensional Euclidean space. We then provide experiments showing that these algorithms perform as well as the non-approximate version for classic classification tasks while achieving a significant speed-up.

Introduction

Hierarchical Clustering (HC) is a ubiquitous task in data science. Given a data set of n points with some similarity or distance function over them, the goal is to group similar points together into clusters, and then recursively group similar clusters into larger clusters. The clusters produced throughout the procedure can be thought of as a hierarchy or a tree with the data points at the leaves and each internal node corresponds to a cluster containing the points in its subtree. This tree is often referred to as a "dendrogram" and is an important illustrative aid in many settings. By inspecting the tree at different levels we get partitions of the data points to varying degrees of granularity. Famous applications are in image and text classification [START_REF] Steinbach | A comparison of document clustering techniques[END_REF], community detection [START_REF] Leskovec | Mining of massive datasets[END_REF], finance [START_REF] Tumminello | Correlation, hierarchies, and networks in financial markets[END_REF], and in biology [START_REF] Breyne | Genome-wide expression analysis of plant cell cycle modulated genes[END_REF][START_REF] Diez | A novel brain partition highlights the modular skeleton shared by structure and function[END_REF].

Perhaps the most popular procedures for HC are Single-Linkage, Average-Linkage, and Ward's method. These are so-called agglomerative HC algorithms (as opposed to divisive) since they proceed in a bottom-up fashion: In the beginning, each data point is in its own cluster, and then the 33rd Conference on Neural Information Processing Systems (NeurIPS 2019), Vancouver, Canada. most similar clusters are iteratively merged -creating a larger cluster that contains the union of the points from the two smaller clusters -until all points are in the same, final cluster.

The difference between the different procedures is in their notion of similarity between clusters, which determines the choice of clusters to be merged. In Single-Linkage the distance (or dissimilarity) is defined as the minimum distance between any two points, one from each cluster. In Average-Linkage we take the average instead of the minimum, and in Ward's method we take the error sum-of-squares (ESS). It is widely accepted that Single-Linkage enjoys implementations that are somewhat simpler and faster than Average-Linkage and Ward's, but the results of the latter two are often more meaningful. This is because its notion of distance is too sensitive and a meaningless "chain" in the data can sabotage the resulting clustering. Extensive discussions of these procedures can be found in many books (e.g. [START_REF] Friedman | The elements of statistical learning[END_REF][START_REF] Leskovec | Mining of massive datasets[END_REF][START_REF] Schütze | Introduction to information retrieval[END_REF][START_REF] Abello | Handbook of massive data sets[END_REF]), surveys (e.g. [START_REF] Murtagh | A survey of recent advances in hierarchical clustering algorithms[END_REF][START_REF] Murtagh | Comments on 'parallel algorithms for hierarchical clustering and cluster validity[END_REF][START_REF] Carlsson | Characterization, stability and convergence of hierarchical clustering methods[END_REF]), and experimental studies (e.g. [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF]).

All of these procedures can be performed in nearly quadratic time, and the main question studied by this paper is whether we can reduce the time complexity to subquadratic. The standard quadratic algorithm for Single-Linkage is quite simple and can be described as follows. After computing the n ˆn distance matrix of the points, we find a minimum spanning tree (MST). This first stage takes Opn 2 dq time if the points are in d-dimensional Euclidean space. In the second stage we perform merging iterations, in which the clusters correspond to connected subgraphs of the MST (initially, each point is its own subgraph). We merge the two subgraph whose in-between edge in the MST is the smallest. By the properties of MST, the edge between two subgraphs (clusters) is exactly the minimum distance between them. This second stage can be done with Opnq insertions, deletions, and minimum queries to a data structure, which can be done in near-linear time. The algorithms for Average-Linkage and Ward's are more complicated since the MST edges between two clusters can be arbitrarily smaller than the average distance or the ESS between them, and we must consider all pairwise distances in clusters that quickly become very large. Nonetheless, an Opn 2 log nq algorithm (following a first stage of computing the distance matrix) has been known for many decades [START_REF] Murtagh | A survey of recent advances in hierarchical clustering algorithms[END_REF].

Can we possibly beat quadratic time? It is often claimed (informally) that Ωpn 2 q is a lower bound because of the fist stage: it seems necessary to compute the distance matrix of the points whose size is already quadratic. More formally, we observe that these procedures are at least as hard as finding the closest pair among the set of points, since the very first pair to be merged is the closest pair. And indeed, under plausible complexity theoretic assumptions 1 , there is an almost-quadratic n 2´op1q lower bound for the closest pair problem in Euclidean space with dimension d " ωplog nq [START_REF] Alman | Probabilistic polynomials and hamming nearest neighbors[END_REF][START_REF] Karthik | On closest pair in euclidean metric: Monochromatic is as hard as bichromatic[END_REF]. This gives a quadratic conditional lower bound for all three Single-Linkage, Average-Linkage, and Ward's method.

Achieving subquadratic runtime has been of interest for many decades (as can be deduced from the survey of Murtagh [START_REF] Murtagh | A survey of recent advances in hierarchical clustering algorithms[END_REF]) and it is increasingly desirable in the era of big data. (See also the recent work on quadratic vs. subquadratic complexity of Empirical Risk Minimization problems [START_REF] Backurs | On the fine-grained complexity of empirical risk minimization: Kernel methods and neural networks[END_REF].)

In this work, we focus on worst-case guarantees while allowing for small approximation in the answers: how fast can we perform these procedures if each iteration is allowed to pick an approximately best pair to merge? More precisely, when merging two clusters the algorithm is allowed to do the following. If the best pair of (available) clusters has (minimum, average, or ESS) distance d then the algorithm can choose any pair of clusters whose distance is between d and γ ¨d, where γ ě 1 is a small constant. When approximations are allowed the time complexity of closest pair drops, and so does the conditional lower bound. Even in high dimensions, Locality Sensitive Hashing techniques can find the γ-approximate nearest neighbors (ANN) in L 1 -distance with n Op1{γq time per query [START_REF] Andoni | Beyond localitysensitive hashing[END_REF][START_REF] Andoni | Optimal data-dependent hashing for approximate near neighbors[END_REF]. This gives a subqadratic n 1`Op1{γq algorithm for closest pair 2 , but can we achieve the same speed-up for γ-approximate Average-Linkage? Namely, can we do Average-Linkage as fast as performing r Opnq (approximate) nearest-neighbor queries?

For the simpler γ-approximate Single-Linkage it is rather easy to see that the answer is yes. This essentially follows from the classical Nearest Neighbor Chain algorithm for HC [START_REF] Murtagh | A survey of recent advances in hierarchical clustering algorithms[END_REF]. Here is a simple way to see why subquadratic is possible in this case: The idea is to replace the expensive first stage of the Single-Linkage algorithm (described above) with an approximate MST computation which can be done in subquadratic time [START_REF] Borodin | Subquadratic approximation algorithms for clustering problems in high dimensional spaces[END_REF][START_REF] Har-Peled | Approximate nearest neighbor: Towards removing the curse of dimensionality[END_REF] using ANN queries. Then we continue to perform the second stage of the algorithm with this tree.

Still it is of great interest to speed up the Average-Linkage and Ward's algorithms since they typically give more meaningful results. This is much harder and before this work, no subquadratic time algorithm for Average-Linkage or Ward's with provable guarantees were known. Various algorithms and heuristics have been proposed, see e.g. [START_REF] Schütze | Projections for efficient document clustering[END_REF][START_REF] Franti | Fast agglomerative clustering using a knearest neighbor graph[END_REF][START_REF] Kull | Fast approximate hierarchical clustering using similarity heuristics[END_REF][START_REF] Otair | Approximate k-nearest neighbour based spatial clustering using kd tree[END_REF][START_REF] Jeon | Nc-link: A new linkage method for efficient hierarchical clustering of large-scale data[END_REF][START_REF] Yildirim | K-linkage: A new agglomerative approach for hierarchical clustering[END_REF], that beat quadratic time by either making assumptions on the data or by changing the merging criteria altogether. Intuitively, while in Single-Linkage only Opnq distances are sufficient for the entire computation (the distances in the MST), it is far from clear why this would be true for Average-Linkage and Ward's.

Our Contribution

Our main result is a γ-approximate Ward's algorithm that runs in subquadratic Õpn 1`Op1{γ 2 q ε ´2 ǹdq time, for any γ ą 1 and ε ą 0, when the points are in d-dimensional Euclidean space. For instance, using the precise bounds of [START_REF] Andoni | Optimal data-dependent hashing for approximate near neighbors[END_REF] on the constants in the Op1{γq term, we can get a 2.35approximation in time Õpn 1.1 ε ´2q. We also prove a similar but slightly slower result for Average-Linkage: our running time for a γ-approximation of Average-Linkage is Õpn 1`Op1{γq `ndq. Moreover, our algorithms are reductions to Õpnq approximate nearest neighbor queries in dimension Õpdq with L 2 distance (for Ward's) or L 1 distance (for Average-Linkage). Thus, further improvements in ANN algorithms imply faster approximate HC, and more importantly, one can use the optimized ANN libraries to speed up our algorithm in a black-box way. In fact, this is what we do to produce our experimental results. Our theorems are as follows.

Theorem 1.1. Given a set of n points in R d and a γ-Approximate Nearest Neighbor data structure which supports insertion, deletion and query time in time T n , there exists a γp1 `εq-approximation of Ward's Method running in time Opn ¨T ¨ε´2 logp∆nqq, where ∆ is the aspect ratio of the point set.

Theorem 1.2. Given a set of n points in R d and a data structure for γ-Approximate Nearest Neighbor under the L 1 -norm which supports insertion, deletion and query time in time T n , there exists a γp1 `εq-approximation of Average Linkage running in time n ¨T ¨ε´2 log Op1q p∆nq, where ∆ is the aspect ratio of the point set.

Our algorithm for approximating Ward's method is very simple: We follow Ward's algorithm and iteratively merge clusters. To do so efficiently, we maintain the list of centroids of the current clusters and perform approximate nearest neighbor queries on the centroids to find the closest clusters. Of course, this may not be enough since some clusters may be of very large size compared to others and this has to be taken into account in order to obtain a γ-approximation. We thus partition the centroids of the clusters into buckets that represents the approximate sizes of the corresponding clusters and have approximate nearest neighbor data structure for each bucket. Then, given a cluster C, we identify its closest neighbor (in terms of Ward's objective) by performing an approximate nearest neighbor query on the centroid of C for each bucket and return the best one.

Our algorithm for Average-Linkage is slightly more involved. Our algorithm adapts the standard Average-Linkage algorithm, with a careful sampling scheme that picks out representatives for each large cluster, and a strategic policy for when to recompute nearest neighbor information. The other sections of this paper are dedicated to explaining the algorithm. Implementation-wise it is on the same order of complexity as the standard Average-Linkage algorithm (assuming a nearest neighbor data structure is used as a black-box), while efficiency-wise it is significantly better as it goes below quadratic time. The gains increase as we increase the tolerance for error, in a controlled way.

We focus our empirical analysis on Ward's method. We show that even for a set of parameters inducing very loose approximation gurantees, the hierarchical clustering tree output by our algorithm is as good as the hierarchical clustering tree produced by Ward's method in terms of classification for most of several classic datasets. On the other hand, we show that even for moderately large datasets, e.g.: sets of 20-dimensional points of size 20000, our algorithm offers a speed-up of 2.5 over the popular implementation of Ward's method of sci-kit learn.

Related Works

A related but orthogonal approach to ours was taken by a recent paper [START_REF] Cochez | Twister tries: Approximate hierarchical agglomerative clustering for average distance in linear time[END_REF]. The authors design an agglomerative hierarchical clustering algorithm, also using LSH techniques, that at each step, with constant probability, performs the merge that average linkage would have done. However, with constant probability, the merge done by their algorithm is arbitrary, and there is no guarantee on the quality of the merge (in terms of average distance between the clusters merged compared to the closest pair of clusters). We believe that our approach may be more robust since we have a guarantee on the quality of every merge, which is the crux of our algorithms. Moreover, they only consider Average-Linkage but not Ward's method.

Strengthening the theoretical foundations for HC has always been of interest. Recently, an influential paper of Dasgupta [START_REF] Dasgupta | A cost function for similarity-based hierarchical clustering[END_REF] pointed to the lack of a well-defined objective function that HC algorithms try to optimize and proposed one such function. Follow up works showed that Average-Linkage achieves a constant factor approximation to (the dual of) this function [START_REF] Cohen-Addad | Hierarchical clustering: Objective functions and algorithms[END_REF][START_REF] Moseley | Approximation bounds for hierarchical clustering: Average linkage, bisecting k-means, and local search[END_REF], and also proposed new polynomial time HC algorithms for both worst-case and beyond-worst-case scenarios that can achieve better approximation factors [START_REF] Roy | Hierarchical clustering via spreading metrics[END_REF][START_REF] Charikar | Approximate hierarchical clustering via sparsest cut and spreading metrics[END_REF][START_REF] Cohen-Addad | Hierarchical clustering beyond the worst-case[END_REF][START_REF] Charikar | Hierarchical clustering better than average-linkage[END_REF][START_REF] Charikar | Hierarchical clustering for euclidean data[END_REF]. Other theoretical works prove that Average-Linkage can reproduce a "correct" clustering, under some stability assumptions on the data [START_REF] Balcan | A discriminative framework for clustering via similarity functions[END_REF]. Our work takes a different approach. Rather than studying the reasons for the widespread empirical findings of the utility of HC algorithms (and mainly Average-Linkage and Ward's), we take it as a given and ask: how fast can we produce results that are as close as possible to the output of Average-Linkage and Ward's. In some sense, the objective function we try to optimize is closeness to whatever Average-Linkage or Ward's produce.

A γ-Approximation of Ward's Method 2.1 Preliminaries

Let P Ă R d be a set of n points. Up to rescaling distances we may assume that the minimum distance between any pair of points is 1. Let ∆ denote the aspect ratio of P , namely ∆ " max u,vPP distpu, vq. Let γ ą 1 be a fixed parameter. Our goal is to build a γ-approximation of Ward's hierarchical clustering.

Let C be a cluster, then define the error sum-of-square as Thus, Ward's algorithm constructs a hierarchy of clusters where each level represents a clustering of the points and where clusters at a given level are subsets of clusters of level `1. Ward's algorithm builds this hierarchy in a bottom-up fashion, starting from n clusters (each point is itself a cluster). Then, given the clustering of a given level , Ward's algorithm obtains the clustering of the next level by merging the two clusters that yield the clustering of minimal ESS. More formally, consider a clustering C " tC 1 , . . . , C u. To find the clustering of minimum ESS obtained by merging a pair of clusters of C, it is enough to minimize the increase in the ESS induced by the merge. Therefore, we want to identify the clusters C i , C j that minimize the following quantity.

∆ESSpC i , C j q " |C i ||C j | |C i | `|C j | ||µpC i q ´µpC j q|| 2 2 . (1)
We will also make use of the following fact. Let γ ą 0 be a parameter, P a set of points in R d . Let D be a data structure that for any set P of n points in R d where d " Oplog nq, supports the following operations. Insertion of a point in P in time Opn f pγq q, for some function f . Deletion of a point in P in time Opn f pγq q; Given a point p P P , outputs a point inserted to the data structure at L 2 ´distance at most γ times the distance from p to the closest point inserted to the data structure, in time Opn f pγq q.

There are data structures based on locality sensitive hashing for f pγq " 1 `Op1{γ2 q, see for example [START_REF] Andoni | Optimal data-dependent hashing for approximate near neighbors[END_REF]. The proof of the following lemma is in the appendix. Lemma 2.1. For any ε ą 0, the above nearest neighbour data structure for clusters with parameters γ, ε, Dpγ, εq has the following properties:

2.
• The insertion time is Opn f pγq ε ´1 log nq;

• On Query(C), it returns a clusters C1 such that ESSpC Y C 1 q ´ESSpCq ´ESSpC 1 q ď p1 `εqγ min BPDpε,γq pESSpC Y Bq ´ESSpCq ´ESSpBqq.

• The query time is Opn f pγq ε ´1 logpn∆qq.

The Main Algorithm

We define the value of merging two clusters A,B as ESSpA Y Bq ´ESSpAq ´ESSpBq. Our algorithm starts by considering each point as its own cluster, together with the Nearest Neighbour Cluster Data Structure described above. Then, the algorithm creates a logarithmic number of rounded merge values that partition the range of possible merge values. Let I be the sequence of all possible merge values in increasing order.

Given a set of n points with minimum pairwise distance 1 and maximum pairwise distance ∆, we have that the total number of merge value β is Oplogpn∆qq.

The algorithm maintains a clustering and at each step decides which two clusters of the current clustering should be merged. The clusters of the current clustering are called unmerged clusters.

The algorithm iterates over all merge values in an increasing order while maintaining the following invariant: Invariant 2.2. When the algorithm reaches merge value δ, for any pair of unmerged cluster C, C 1 we have ESSpC Y C 1 q ´ESSpCq ´ESSpC 1 q ě δ{γ.

We now give a complete description of our algorithm.

B. Remove N N pCq from ToMerge and add C 1 to ToMerge; µpC 1 q follows immediately from µpCq, µpN N pCqq, |C| and |N N pCq| (see Fact 1) C. Remove C, N N pCq from L and add C 1 to L

The running time analysis and proof of correctness of the algorithm are deferred to the appendix.

3 A γ-Approximation of Average-Linkage

Preliminaries

For two sets of points A, B, we let avgpA, Bq "

1 |A||B| ř aPA ř
bPB dpa, bq. The following simple lemma is proved in the appendix. Lemma 3.1. Consider three sets of points A, B, C. We have that avgpA, Cq " avgpC, Aq ď avgpA, Bq `avgpB, Cq

Overview and Main Data Structures

Our goal is to design a γ-approximate Average-Linkage algorithm. The input is a set P of n points in a d-dimensional Euclidean space. The algorithm starts with a clustering where each input point is in its own cluster. The algorithm then successively merges pairs of clusters. When two clusters are merged, a new cluster consisting of the union of the two merged clusters is created. The unmerged clusters at a given time of the execution of the algorithm are the clusters that have not been merged so far. More formally, at the start the set of unmerged clusters is the set of all clusters. Then, whenever two clusters are merged, the newly created cluster is inserted to the set of unmerged clusters while the two merged clusters are removed from the set. The algorithm merges clusters until all the points are in one cluster.

To be a γ-approximation to Average-Linkage, our algorithm must merge clusters according to the following rule: If the minimum average distance between a pair of unmerged clusters is v then the algorithm is not allowed to merge two unmerged clusters with average distance larger than γ ¨v.

Let ε ą 0 and γ ě 1 be parameters. We will show how to use a γ-approximate nearest neighbor data structure (on points) to get a γ 1 -approximate Average-Linkage algorithm where γ 1 " p1 `εq ¨γ.

We make use of the following key ingredients.

• We design a sampling scheme that allows to choose at most poly log n points per cluster while preserving the average distance up to p1`εq-factor with probability at least 1´1{n 5 . • We design a data structure that given a set of clusters, allows to answer approximate nearest neighbor queries (on clusters) according to the average distance. • Finally we provide a careful scheme for the merging steps that allows to bound the number of times the nearest neighbor queries for a given cluster have to be performed.

The Algorithm

We are now ready to describe our algorithm. Our algorithm starts with all input points in their own clusters and performs a nearest neighbor query for each of them. The algorithm maintains a partition of the input into clusters that we call the unmerged clusters, identical to average linkage. The algorithm proceeds in steps. Each step consists of merging several pairs of clusters. For each step we associate a value v, which we refer to as the merge value of the step, which is a power of p1 `εq and we will show the invariant that at the end of the step associated with value v, the unmerged clusters are at distance greater than v{pp1 `εq 2 γq.

For each cluster C, we will maintain a sample of its points by applying the sampling procedure (see supplementary materials for more details). To avoid recomputing a sample too often, we set a variable spCq which corresponds to the size of the cluster the last time the sampling procedure was called.

Lazy sampling. Every time two clusters C 1 , C 2 are merged by the algorithm to create a new cluster, the following operations are performed:

1. If |C 1 Y C 2 | ě p1 `ε2 {p1 `γqq maxpspC 1 q, spC 2 qq, then the sampling procedure is called on C 1 Y C 2 and an approximate nearest cluster query is performed using the nearest cluster data structure (see supplementary materials). Then,

spC 1 Y C 2 q is set to |C 1 Y C 2 |.
The resolution parameter for sampling is the value of the current step divided by n. Namely, if the value of the current step is v, we set α C1YC1 " v for the sampling procedure. 2. Otherwise, spC 1 Y C 2 q is set to maxpspC 1 q, spC 2 qq and the algorithm uses the sample of argmax CPtC1,C2u |C| as the sample for

C 1 Y C 2 .
Once the above has been performed, a γ-approximate nearest cluster query is performed using the sample defined for the cluster resulting of the merge.

Thus, at each step, all the clusters have a γp1 `Opεqq-approximate nearest neighbor among the clusters. We denote ν t pCq the approximate nearest neighbor for cluster C at the tth step. This approximate nearest neighbor is computed using our data structure (see supplementary materials). We let νpCq " ν tpCq pCq, where tpCq is the step at which C was created.

We say that a cluster C is active at a step of value v if v ě avgpC, νpCqq{n 2 . For all steps at which a cluster is not active, the algorithm does not recompute an approximate nearest cluster and the algorithm keeps the approximate nearest cluster computed at the creation of the cluster. This is a crucial point for making our algorithm run in subquadratic time. If this cluster has been merged, a pointer to the resulting cluster is kept instead. During the steps at which the cluster is active, a γ-approximate nearest cluster query is performed.

Algorithm at a given step. A step of value v is performed when the minimum distance between pairs of approximate nearest clusters is at least v{γ and no more than p1 `εqv. At a given step of value v, the algorithm is as follows. The clusters are processed one by one. Each cluster is inserted into our nearest cluster data structure (see supplementary materials) and an approximate nearest cluster query is performed. If the approximate nearest clusters returned by the data structure is at distance at most v, then the two clusters are merged and the following operations are performed: the two clusters are removed from the data structure and the new cluster is inserted into the list of the clusters to be processed. Then the algorithm continues to process the list.

See supplementary materials for the proof of correctness.

Experiments

Our experiments focus on Ward's method and its approximation since it is a simpler algorithm in contrast with average-linkage. We implemented our algorithm using C++11 on 2.5 GHz 8 core CPU with 7.5 GiB under the Linux operating system. Our algorithm takes a dynamic Nearest Neighbour data structure as a block box. In our implementation, we are using the popular FLANN library [START_REF] Muja | Scalable nearest neighbor algorithms for high dimensional data[END_REF] and our own implementation of LSH for performing approximate nearest neighbor queries. We compare our algorithm to the sci-kit learn implementation of Ward's method [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF] which is a Python library that also uses C++ in the background.

Our algorithm has different parameters for controlling the approximation factor. These parameters have a significant effect on the performance and the precision of the algorithm. The main parameter that we have is which determines the number of data structures to be used (recall that we have one approximate nearest neighbor data structure for each p1 `εq i , for representing the potential cluster sizes) and the sequence of merge values. Moreover, we make use of FLANN library procedure for finding approximate nearest neighbors using KD-trees. This procedure takes two parameters the number of trees t and the number of leaves visited f . The algorithm builds t randomized KD-trees over the dataset. The number of leaves parameter controls how many leaves of the KD-trees are visited before stopping the search and returning a solution. These parameters control the speed and precision of the nearest neighbor search. For instance, increasing the number of leaves will lead to a high precision but at the expense of a higher running time. In addition, decreasing the number of KD-Tree increases the performance but it decreases the precision. For LSH, we use the algorithm of Datar et al. [START_REF] Datar | Locality-sensitive hashing scheme based on p-stable distributions[END_REF] which has mainly two parameters, H the number of hash functions used and r controlling the 'collision' rate (see details in [START_REF] Datar | Locality-sensitive hashing scheme based on p-stable distributions[END_REF]).

To study the effects of these parameters, we did different experiments that combine several parameters and we report and discuss the main results in We report the normalized mutual information score of the clustering output by the different algorithms compared to the ground-truth labels for each dataset. We note that 0.05 can obtained on Newsgroup through a random labelling of the vertices (up to ˘0.02). Hence LSH seems a more robust approach for implementing approx-ward.

experiments are classic real-world datasets from the UCI repository and the sci-kit-learn library. Iris contains 150 points in 4 dimensions, Digits 1797 in 64 dimensions, Boston 506 points in 13 dimensions, Cancer 569 points in 3 dimensions, and Newsgroup 11314 points in 2241 dimensions.

To measure the speed-up achieved by our algorithm, we focus our attention on a set of parameters which gives classification error that is similar to Ward's on the real-world datasets, and then run our algorithm (with these parameters) on synthetic dataset of increasing sizes. These parameters are precisely " 8, number of trees T " 2, the number of visited leaves L " 10. The datasets are generated using the blobs procedure of sci-kit learn. The datasets generated are d-dimensional for d " t10, 20u and consists of a number of points ranging from 10 000 to 20 000. In both dimensions, we witness a significant speed-up over the sci-kit learn implementation of Ward's algorithm. Perhaps surprisingly, the speed-up is already significant for moderate size datasets. We observe that the running time is similar for LSH or FLANN.

(a) Running time of our algorithm with parameters (ε = 8, T = 2, L = 10) (in red) and of Ward's method, on datasets of sizes ranging from 10 000 points to 20 000 points in R 10 . We observe that our algorithm is more than 2.5 faster on datasets of size 20 000.

(b) Running time of our algorithm with parameters (ε = 8, T = 2, L = 10) (in red) and of Ward's method, on datasets of sizes ranging from 10 000 points to 20 000 points in R 20 . We observe that our algorithm is more than 2.5 faster on datasets of size 20 000. Interestingly, it seems that the dimension has little influence on both our algorithm and Ward's method.

A.3.2 A Data Structure for Approximate Nearest Cluster

In this section, we introduce a data structure for finding approximate nearest clusters. The following theorem is proved in the appendix. Theorem A.5. Let γ ą 0 be a parameter, P a set . Let D be a data structure that for any set P of n points in R d where d " Ωplog nq, supports the following operations:

1. Insertion of a point in P in time Opn f pγq q, for some function f ;

2. Deletion of a point in P in time Opn f pγq q;

3. Given a point p P P , outputs a point inserted to the data structure at L 1 -distance at most γ times the distance from p to the closest point inserted to the data structure, in time Opn f pγq q.

Then, for any ε ą 0, there exists a data structure for pairs pS, wq where S is a set of points in R d and w is a positive value, that supports the following operations:

1. Insertion of a pair (set, value) in time Opηε ´1 log n ¨nfpγq q;

2. Deletion of a pair (set, value) in time Opηε ´1 log n ¨nfpγq q;

3. Given a set of points C in R d and a value w, outputs a pair pC 1 , w 1 q inserted to the data structure that is such that that avgpC, C 1 q `w `w1 is at most γp1 `εq times min pC ˚,w ˚q in the data structure avgpC, C ˚q `w `w˚i n time Opηε ´1 log n ¨nfpγq q.

Proof of Lemma A.3. We claim that we can simply use a constant factor approximation to the median problem to find c i -there is a vast literature of near-linear algorithms producing an Op1qapproximation to the median.

Consider the median of P , namely the point p ˚P P that minimizes ř pPP distpp, p ˚q. We have that

1 |P |´1
ř pPP distpp, p ˚q is at most avgpP q. Thus, consider any point p that is an Op1q-approximation to the median of P . We have that 1

|P |´1

ř pPP distpp, pq " OpavgpP q. Then, the remaining step of the sampling procedure is to evaluate the distance from each point to p to define R i and G i . This can be done in linear time. Finally, the sampling of points in R i can also be done in linear time.

Proof of Lemma A.4. We have, by Lemma 3.1,

avgpC i , C j q " |R i | |C i | avgpR i , C j q `|G i | |C i | avgpG i , C j q ď avgpR i , C j q `|G i | |C i | avgpR i , G i q ď avgpR i , C j q `ε ¨avgpR i , G i q ď avgpR i , C j q `ε ¨avgpC i q Similarly, we have avgpR i , C j q " |R j | |C j | avgpR i , R j q `|G j | |C j | avgpG j , R i q ď avgpR i , R j q `|G j | |C j | avgpR j , G j q
ď avgpR i , R j q `ε ¨avgpR j , G j q ď avgpR i , R j q `ε ¨avgpC j q

Combining yields avgpC i , C j q ď avgpR i , R j q `ε ¨avgpC j q `ε ¨avgpC i q.

Therefore, by applying Corollary 1 to SpC i q, SpC j q, we have that avgpSpC i q, SpC j q " p1 εqavgpR i , R j q and so avgpC i , C j q ď p1 `εq x avgpC i , C j q since the diameter of the points in R i and R j is at most avgpC i q{ε and avgpC j q{ε respectively.

We now aim at proving that avgpC i , C j q ě p1 ´Opεηqq x avgpC i , C j q. Recall that by assumption, we have that avgpC i q, avgpC j q ď η ¨avgpC i , C j q. Thus, again combining with Corollary 1, we have that x avgpC i , C j q ď p1 `εqavgpR i , R j q `2εη ¨avgpC i , C j q. Moreover, as discussed above, we have that avgpC i , C j q ě p1 ´OpεqqavgpR i , R j q and so, rescalling ε, we have x avgpC i , C j q ď p1 `εqavgpC i , C j q, as claimed.

Proof of Theorem A.5. We start with some preprocessing steps and notations. We consider an isometric embedding of all the input points into L 1 with distortion at most p1 `εq, for some sufficiently small ε ą 0. In the remaining, we thus work with the L 1 norm.

For each point p, for each integer i, let p i " p 1 ¨p1 . . . p 1 loooooomoooooon i Namely, the coordinates of p i are obtained by concatenating the coordinates of p i times. Given a set of j points S " tp 1 , p 2 , . . . , p j u and a value w S , we let q i pSq be the point in a pi ¨j ¨d `2q-dimensional space with coordinates p i 1 , p i 2 , . . . p i j , 0, i ¨j ¨wS . Namely, q i pSq is obtained by concatenating p i of all the j points p P S, adding an extra coordinate of value 0 and adding a final coordinate with value i ¨j ¨wS . We also let d i pSq be the point in a pj ¨i ¨d `2q-dimensional space with coordinates p 1 , p 2 , . . . p j , p 1 , p 2 , . . . p j , . . . p 1 , p 2 , . . . p j loooooooooooooooooooooooooomoooooooooooooooooooooooooon i¨j , i ¨j ¨wS , 0. Namely obtained by concatenating the coordinates of the point p 1 , p 2 , . . . p j , i times, adding an extra coordinate of value i¨j ¨wS and adding a final coordinate with value 0. We have the following claim, whose proof follows immediately from the definition. Claim 1. Given two sets A and B, of size i and j respectively, and two values w A , w B , we have that

1 i ¨j ||q j pAq ´di pBq|| 1 " w A `wB `1 i ¨j ÿ aPA ÿ bPB ||a ´b|| 1 .
We now describe our data structure using an approximate nearest-neighbor data structure D for the L 1 distance between points. We make use of an approximate nearest neighbor data structure D i,j,k , for each integers i, j P t1, 2, . . . , ηu.

Let C be a cluster. The insertion is as follows. Let i " |C|. The algorithm inserts the point d j pCq in the data structure D i,j , for all j P t1, 2, . . . , ηu. Deletion of C consists of removing d j pCq from the D i,j it has been inserted into. The time complexities for insertion and deletion follow immediately.

The approximate nearest neighbor query for cluster C is performed as follows. For all j P t1, 2, . . . , ηu, the algorithm creates the point q j pCq, and makes a nearest neighbor query in the data structure D i,j . Let p j be the point returned by the query q j pCq on data structure D i,j and ν j be the cluster corresponding to p j . Claim 1 implies that 1 |C|¨|ν j | ||q j pCq ´di pν j q|| 1 " p1 ˘εqpw C `wν j `avgpC, ν j qq.

Then, let j ˚" argmin j avgpC, ν j q. We now argue that avgpC, ν j ˚q `wC `wν j ˚ď γp1 ὲq min C 1 ‰C pavgpC, C 1 q `wC `wC 1 q.

Let Ĉ " argmin C 1 ‰C pavgpC, C 1 q `wC `wC 1 q and ĵ " | Ĉ|. Consider the data structure D i, ĵ . By its correctness, D i, ĵ returned a point p ĵ such that ||q ĵ pCq ´pĵ || 1 ď γp||q ĵ pCq ´di p Ĉq|| 1 . Thus, applying Claim 1 yields that avg w pC, ν ĵ q `wC `wν ĵ ď γp1 `εqpavgpC, Ĉq `wC `w Ĉ q. By the choice of j ˚, we thus have that avgpC, ν j ˚q ď γavgpC, Ĉq, as claimed.

Invariant. The correctness of the algorithm is captured by the following invariant. The proof, as well as the running time analysis, are deferred to the appendix. Lemma A.6 (Invariant for correctness). The following holds with probability at least 1 ´1{n 3 . Consider the tth step of the algorithm, let v be the merge value at the tth step.

1. At the end of the step, no cluster at (inner) average distance greater than vp1 `εq has been merged by the algorithm so far. For any unmerged clusters C i , C j , we have that x avgpC i , C j q " p1 `OpεqqavgpC i , C j q.

2. For any unmerged cluster C at the end of the step, ν t pCq is an unmerged p1 `Opεqqγapproximate nearest cluster of C.

3. Finally, at the end of a step of value v, there is no pair of clusters at average distance less than v{pp1 `εq 2 γq.

Proof of Lemma A.6. We prove it by induction on the number of steps of the algorithm. This is clearly true at first.

We start with (1). For simplicity, assume that first that the algorithm does not do lazy sampling and runs the sampling procedure after each merge. Then, (1) follows from the definition of the algorithm and the inductive hypothesis on the correctness of the sampling procedure (Lemma A.4). More formally, the definition of the algorithm ensures that no pair of clusters C i , C j such that x avgpC i , C j q ą vp1 `εq are merged Moreover, by the inductive hypothesis, we have that for any cluster C, avgpCq ď vp1 `εq.

Thus, we can apply Lemma A.4 with η " p1 `εq and we deduce that for any pair of clusters C i , C j ,

x avgpC i , C j q " p1 ˘εqavgpC i , C j q with probability at least 1 ´1{n 5 . Taking a union bound over all n steps and n merges of the algorithm and all Opn 2 q pairs of clusters in total concludes the proof of (1) in the case of non-lazy sampling.

To finish the proof of (1), we need to show that lazy sampling does not degrade the quality of the outcome of the sampling by too much. Hence, consider an unmerged cluster resulting from the merge possibly at a previous step of two clusters C 1 , C 2 . If |C 1 Y C 2 | ě p1 `ε2 {p1 `γqq maxpspC 1 q, spC 2 qq, then the sampling procedure is applied and the average distance between the samples of C 1 Y C 2 and any other cluster C 3 is within a p1 `εq factor from the average distance between C 1 Y C 2 and C 3 with probability at least 1 ´1{n 4 and the above analysis applies. Now, if |C 1 Y C 2 | ă p1 `ε2 {p1 `γqq maxpspC 1 q, spC 2 qq, then assume w.l.o.g. that |C 1 | ě |C 2 |. Hence, we have that by Lemma 3.1 that for any other unmerged cluster C 3 avgpC 2 , C 3 q ď avgpC 2 , C 1 q `avgpC 1 , C 3 q. Now, by the inductive hypothesis, we have that avgpC 2 , C 1 q ď γavgpC 1 , C 3 q and so avgpC 2 , C 3 q ď p1 `γqavgpC 1 , C 3 q. It follows that avgpC 1 Y C 2 , C 3 q ď p1 `εqavgpC 1 , C 3 q. Finally, by the induction hypothesis, we have that the sample of C 1 preserves the distance from C 1 to C 3 with probability at least 1 ´1{n 4 up to a p1 `εq factor. Thus, we indeed have that the average distance between the samples of any pair of unmerged clusters is within a factor p1 `εq of the average distance of the pair.

We then turn to (3), thus consider the end of a step. Observe that if there are two clusters C 1 , C 2 that are at pairwise distance less than v{pp1 `εq 2 γq then by the inductive hypothesis, the sampling procedure guarantees the two samples for C 1 , C 2 are at average distance at most v{γ. Therefore, a γ-approximate nearest cluster query returns a cluster at distance less than v. Thus, consider the cluster, say C 2 , that is inserted into the data structure last. When C 2 is processed, a nearest neighbor query is performed and so, since the cluster C 1 has been inserted first in the data structure, C 2 should have had an approximate nearest neighbor at distance less than v and so should have been merged.

We now move to prove (2): We finish by considering unmerged clusters at step t. We show that for any unmerged cluster C, the nearest cluster is at average distance at least 1 γ avgpC, νpCqq and at most p1 `1{nqavgpC, νpCqq. This will conclude the proof of the invariant.

Let i be the step at which C is created. Let C ˚be the nearest cluster to C at the tth step. By Theorem A.5, Lemma A.4 and the inductive hypothesis of the γ-approximate nearest neighbor procedure we have that avgpνpCq, Cq ď γavgpC ˚, Cq. Since the unmerged clusters at step t ą i are the union of the clusters of C i , we have that avgpC, C 0 q ě avgpC, C ˚q for any cluster C 0 of C t0 . It follows that for any i 1 ě i, the cluster of C i 1 that is the nearest to C is at distance at least γ ´1 ¨avgpC, νpCqq.

  T px ´µpCqq where µpCq " 1 |C| ř xPC x. We let the error sum-of-square of a clustering C " tC 1 , . . . , C u be ESSpCq " ÿ CPC ESSpCq.

Fact 1 .

 1 Given two set of points A, B with corresponding centroids µpAq, µpBq respectively, we have that the centroid of A Y B is on the line joining µpAq to µpBq, at distance |B| |AYB| ||µpAq ´µpBq|| 2

  1.1 Finding The Nearest Neighbour ClusterOur algorithm relies on a Nearest Neighbour Data Structure for clusters, where the distance between two clusters A, B is given by ESSpA Y Bq ´ESSpAq ´ESSpBq. Given a parameter ε ą 0, our Nearest Neighbour Data Structure Dpγ, εq for clusters consists of Opε ´1 log nq Nearest Neighbour Data Structures for points defined as follows. There is a data structure D for each P tp1 ` q i | i P r1, . . . , log 1` nsu. The data structure works as follows. Insertion(C): Inserting a cluster of a set C of points is done by inserting µpCq in the D i such that p1 `εq i´1 ď |C| ă p1 `εq i .

Query(C): For each i P tp1 ` q i | i P r1, . . . , log 1` nsu perform a nearest neighbor data query for µpCq in D i , let N N i pCq be the result. Output N N i pCq that minimizes ∆ESS C,N NipCq .

Table 1 .

 1 The main data that is used in these

		Iris Cancer Digits Boston Newsgroup
	Ward's	0.67	0.46	0.82	0.80	0.146
	Ward-FLANN ( = 0.5, T = 16, L = 5) 0.62	0.53	0.79	0.80	ă 0.05
	Ward-FLANN ( = 4, T = 16, L = 128) 0.76	0.47	0.56	0.78	ă 0.05
	Ward-FLANN ( = 8, T = 2, L = 10)	0.75	0.51	0.47	0.80	ă 0.05
	Ward-LSH ( " 10, r " 3, H " n 1{10 ) 0.69	0.58	0.58	0.82	ă 0.05
	Ward-LSH ( " 10, r " 3, H " n 1{2 ) 0.72	0.48	0.73	0.83	0.104
	Ward-LSH ( " 2, r " 3, H " n 1{2 )	0.72	0.57	0.63	0.83	0.113

Table 1 :

 1

These lower bounds hold under the Strong Exponential Time Hypothesis of Impagliazzo and Paturi[23,

[START_REF] Impagliazzo | Which problems have strongly exponential complexity[END_REF] regarding the complexity of k-SAT.2 On the negative side, we know that a p1 `εq approximation requires quadratic time[START_REF] Rubinstein | Hardness of approximate nearest neighbor search[END_REF].

from µpAq.

Let L be the list of unmerged clusters, initially it contains all the points.

For each ν P I: (a) ToMerge Ð L (b) While ToMerge is not empty: i. Pick a cluster C from ToMerge, and remove it from ToMerge. ii. N N pCq Ð Approximate Nearest Neighbour Cluster of C. iii. If ESSpC Y N N pCqq ´ESSpCq ´ESSpN N pCqq ď ν: A. Merge C and N N pCq. Let C 1 be the resulting cluster.

Acknowledgements. Ce projet a bénéficié d'une aide de l' État gérée par l'Agence Nationale de la Recherche au titre du Programme FOCAL portant la référence suivante : ANR-18-CE40-0004-01.

A Missing Proofs and details

A.1 Proofs for the approximate Ward's algorithm Proof of Lemma 2.1. The running time follows almost immediately from the definition: there are Opε ´1 log nq data structure to query. The correctness results from the following argument. Consider the cluster C ˚that has been inserted to the data structure and that minimizes min C0 inserted ∆ESSpC, C 0 q. Let j be the integer such that p1 `εq j´1 ď |C ˚| ď p1 `εq j . Consider the cluster C j returned by the query on D j . We have that |C j | ď p1 `εq|C ˚| and so by the correctness of the data structure ∆ESSpC j , Cq ď γp1 `εq∆ESSpC, C ˚q and the lemma follows.

A.2 Runtime analysis and correctness for the approximate Ward's algorithm Running Time The outer loop of Algorithm 1 iterates β times. The total number of clusters created by the algorithm is Opnq where n is the total number of input points. Thus, The inner for loop takes Opnq times. By Lemma 2.1, the body of the inner loop will have at most the complexity of the nearest neighbour search Opn f pγq ε ´1 logpn∆qq. Summing up all these complexities results in Opn 1`f pγq ε ´1 logpn∆qq.

Proof of Correctness

Lemma A.1. Invariant 2.2 holds.

Proof. We proceed by induction on the merge ν. When the merge value is 1, the invariant trivially holds. Now assume that the invariant holds up to some merge value ν. We first show that there is no pair of clusters C i , C j with ∆ESSpC i , C j q ă ν{γ at the end of the iteration corresponding to merge value ν. Assume toward contradiction that this wasn't the case and consider the cluster of C i , C j that was created the last, say C i . Then, a nearest neighbor cluster query was made on C i and since C j was already in the data structure, Lemma 2.1 implies that the query returned a cluster of C such that ∆pC , C i q ă ν. Hence C i was merged to C and not an unmerged cluster at the end of the iteration.

A.3 Proofs for the approximate Average-Linkage algorithm

Proof of Lemma 3.1. Let U " |A||C||B|. We note that for each a P A, c P C, the triangle inequality implies that dpa, cq ď min bPB pdpa, bq `dpb, cqq and so dpa, cq ď 1 |B| ř bPB pdpa, bq `dpb, cqq.

Approximating Cluster Distance by Sampling

Let C 1 , . . . , C k be a collection of clusters. Let n 2 α i be an upper bound on the average distance between points within C i . Assume that the minimum average distance between any pair of clusters is at least α i {n 2 for all i. For each cluster C i , we make a slight abuse of notation and let avgpC i q denote the average distance between points in C i (i.e.: avgpC i q " avgpC i , C i q). Let c i be a point such that avgpc i , C i q ď avgpC i q{ε and let R i denote the points of C i whose distance to c i is at most avgpC i q{ε 2 . In other words,

We consider the following sampling scheme. Among the points in R i , pick ηε ´6 log 3 n points uniformly at random. Let κ i " avgpG i , R i q. By an immediate averaging argument we have that

We make use of the following lemma by Chen [START_REF] Chen | On coresets for k-median and k-means clustering in metric and euclidean spaces and their applications[END_REF].

Lemma A.2 ([13], Lemma 3.3). Let V be a set of points in a metric space pX, dq, and let λ 1 , ξ ą 0 be given parameters. Let ∆ be the diameter of V . Let U be a sample of size ξ ´2 lnp2{λ 1 q points of V picked independently and uniformly, where each point of U is assigned weight |V |{|U | such that ř uPU wpuq " |V |. For a fixed point p, where p is not necessarily a an element of V , we have that | ř vPV distpv, pq ´řuPU wpuqdistpu, pq| ď ξ|V |∆, with probability at least 1 ´λ1 .

From this, we deduce the following corollary.

Corollary 1. Let V be a set of points in a metric space pX, dq, and let λ 1 , ξ ą 0 be given parameters.

Let ∆ be the diameter of V . Let U be a sample of size ξ ´2 lnp2{λ 1 q points of V picked independently and uniformly. For a fixed point p, where p is not necessarily a an element of V , we have that |avgpV, pq ´avgpU, pq| ď ξ∆, with probability at least 1 ´λ1 .

The proof of the following lemma is in the appendix.

Lemma A.3. Given a set of point C i of size m, the sampling procedure can be performed in time Opm{ε 5 q.

For any two clusters C i , C j let SpC i q, SpC j q denote the set of points sampled by the above procedure. Furthermore, we define x avgpC i , C j q " avgpSpC i q, SpC j qq `εκ i `εκ j . We then have the following crucial lemma, proved in the appendix.

Lemma A.4. Consider a set of clusters tC 1 , . . . , C u such that for any pair of clusters C i , C j , avgpC i q, avgpC j q ď ηavgpC i , C i q for some constant η.

Then, by taking a sampling of size 10ηε ´6 log 3 n, we have x avgpC i , C j q " p1 ˘εqavgpC i , C j q with probability at least 1 ´1{n 5 .

We now show an upper bound on the distance to the cluster C 1 containing νpCq. This follows from applying Lemma 3.1 as follows. Consider the sequence of merges that involve νpCq. Let νpCq Ă νpCq 1 Ă . . . Ă νpCq k denote the clusters that contain νpCq and that are successively merged after step i and until time t. By Lemma 3.1, we have that avgpC, νpCq 1 q ď avgpC, νpCqq àvgpνpCq, νpCq 1 q ď avgpC, νpCqq `avgpC, νpCqq{n 2 since C is not active. Similarly, by the inductive hypothesis (1), avgpC, νpCq 2 q ď avgpC, νpCq 1 q `avgpνpCq 1 , νpCq 2 q. Here again, C is not active and so avgpνpCq 1 , νpCq 2 q ď avgpC, νpCqq{n 2 . Since the overall number of merges is at most n, we conclude that avgpC, νpCq k q ď `avgpC, νpCqq `avgpC, νpCqq{n as claimed.

Therefore, the invariant also holds and so the inductive hypothesis is satisfied.

A.4 Running Time Analysis for the approximate Average-Linkage algorithm

We need to bound the number of times an approximate nearest cluster query is performed, the total time incurred by the sampling procedure, the running time of a step, and the number of steps. This is the purpose of the following section.

A.4.1 Sampling Time

Lemma A.7 bounds the total running time incurred by the sampling procedure. Lemma A.7. The total running time caused by the sampling procedure over the entire execution of the algorithm is at most Opn 1`ρ ε ´2γ log nq.

Proof. The lemma follows from Lemma A.3 and due to the fact that the procedure is only called on clusters resulting from the merge of two clusters

γqq maxpspC 1 q, spC 2 qq. Thus, the number of clusters in which an input point can contribute to the running time of the sampling procedure is Opε ´2γ log nq.

Running Time of a

Step At a given step associated with a certain merge value v, the goal is to merge all clusters whose nearest neighbor is at distance at most v so that at the end of the step, the distance from each cluster to its approximate nearest neighbor is greater than v. Let n v be the number of active clusters at the beginning of the step. Lemma A.8. The total number of nearest neighbor queries made by the algorithm during a step with merge value v is Opn v q.

Proof. Observe that the total number of merges is at most Opn v q. Moreover the total number of nearest neighbor queries is bounded by the total number of merges plus the number of active clusters and so at most Opn v q.

A cluster can remain active throughout the entire algorithm. Hence, the number of active step is a priori only bounded by Opε ´1 log ∆nq which gives the claimed complexity.

A slightly more involved algorithm allows to remove the dependency in log ∆ at the price of a slightly worse approximation guarantee: we were only able to show a γ 2 -approximation instead of a γ-approximation in this case. We defer this to the full version of the paper.


S Codis

BEYOND THE POWER SPECTRUM WITH LARGE DEVIATION THEORY

Keywords: cosmology, theory -large-scale structure of Universe -methods, analytical, numerical

A large-deviation principle is used to model the time-evolution of the large-scale structure of the Universe. This approach allows for analytical predictions in the mildly non-linear regime, beyond what is commonly achievable via other statistics such as correlation functions. The idea is to measure the mean cosmic densities within concentric spheres and study their joint statistics.The spherical symmetry then leads to surprisingly accurate predictions where standard calculations of perturbation theory usually break down. Results for the one-point statistics of the cosmic density field are shown and implications for future large galaxy surveys are discussed.

Introduction

In the first stages of structure formation, cosmic fields evolve linearly and therefore remain Gaussian. They are thus completely described by their two-point correlation function (or power spectrum in Fourier space). Those primordial Gaussian inhomogeneities then grow under the laws of gravity in our expanding Universe and give rise to the large-scale cosmic web. This subsequent evolution on mildly non-linear scales generates gravitational non-Gaussianities and information then leaks from the power spectrum to higher order statistics inducing nonzero contributions to the whole hierarchy of N-point correlation functions. The above mentioned time evolution of the large-scale structure of the Universe is governed by the so-called Vlasov-Poisson system, which for a selfgravitating collisionless fluid (our hypothesis here) simply expresses that the volume in phase space is conserved (Liouville theorem) and the gravitational potential is sourced by density fluctuations (Poisson equation). This highly non-linear set of equations can be solved using numerical simulations (with N-body methods or resolving the full phase space) or analytically in some specific regimes. Let us emphasize that efforts on the analytical side are not vain as they allow us to understand the details of structure formation and provide important tools to analyse observational datasets, in particular when computing power becomes a limitation. This is notably the case when covariance matrices need to be finely modelled, therefore requiring thousands of numerical realisations of the observable Universe, which is far beyond the reach of our existing facilities. Eventually, hybrid approaches combining simulations and first principles calculations will probably be the optimal way to extract cosmological information from the large-scale structure of the Universe. In the context of covariance matrices, the response function formalism is a nice illustration of how analytical methods can reduce the computing needs.

In this proceeding, we focus on the analytical modelling of the distribution of cosmic fields. Sect. 2 first introduces perturbative approaches while Sect. 3 presents a more successful non-perturbative method based on large-deviation theory. Finally, Sect. 4 wraps up.

Cosmological perturbation theories

N-point correlation functions

On large scales or at early times when the fields are weakly non-Gaussian, perturbation theory (hereafter PT) techniques can be implemented in order to derive predictions for any observable that depends on the density and velocity fields. This is achieved by first taking the first two moments of the Vlasov equation and closing the SF2A 2018 system assuming there is no shell crossing, therefore no pressure term. Then the fields can be expanded with respect to the initial fields and can eventually be computed hierarchically (order by order) using convolutions of the successive PT kernels (for a review, see Bernardeau et al. 2002). This approach is valid only where the density contrast is small δ 1. Once applied to predict the power spectrum, it allows to reach scales k of a few 0.1 h/Mpc at most. Going deeper into the non-linear regime is extremely difficult due to the slow convergence of the perturbative expansion on the one hand (coupling with small non-linear scales) and to the limitation of the no-shell crossing condition on the other hand. The price to pay is often to introduce a myriad of free parameters that encode our ignorance of the small-scale physics (as is consistently done for instance in the EFTofLSS approach). However, in the context of cosmological inference, it is still unclear whether there is any gain in going to smaller scales (hence having more statistics given the number of modes available in this regime) if many nuisance parameters have to be introduced and eventually marginalised over.

In particular, non-linearities make optimal extraction of information a difficult problem and therefore the question of which observables to use becomes crucial. The current most common approach is to start from the power spectrum which contains all the information in the Gaussian case and go beyond with higher order Npoint correlation functions. In practice, the bispectrum is already quite challenging and any analysis is usually restricted to the power spectrum, possibly the bispectrum but no more. An alternative route could be to look for other observables which can be predicted from first principles and can probe the mildly non-linear regime (i.e beyond the weakly non-linear regime but before shell crossing). To achieve this goal, one may think that imposing some level of symmetry could potentially provide us with observables that are less prone to the effect of the unknown small scales. To date, one such configuration has been evidenced: count-in-cells.

Count-in-cells : a perturbative approach

The idea behind count-in-cells is to draw spheres of a given size and measure the mean cosmic matter densities within them. The distribution of these mean densities P(δ) (i.e their one point statistics) is initially Gaussian in the first stages of structure formation and subsequently develops some level of non-Gaussianity being in particular increasingly skewed towards small densities as voids are occupying more and more volume and peaks are getting more and more concentrated. Using a perturbative approach, this probability distribution function (PDF) of the density field can be written as an Edgeworth expansion (at second order here) around a Gaussian kernel

P(δ) = G(δ) 1 + σ S 3 3! H 3 δ σ + σ 2 S 4 4! H 4 δ σ + 1 2 S 3 3! 2 H 6 δ σ + • • • (2.1)
where H i are probabilists' Hermite polynomials. This Edgeworth expansion involves the successive field cumulants which in PT are known to scale as increasing power of the amplitude of fluctuation σ, δ i = S i σ 2(i-1) , where the reduced cumulants S i are independent of σ at leading order. As an example, Fig. 1 displays the true PDF with solide lines while a second order Edgeworth expansion is superimposed with dashed lines. The various colours show different redshift. As expected the PDF is initially close to Gaussian and becomes more and more skewed with cosmic time. As non-linearities develop, the Edgeworth expansion performs increasingly badly as expected. An important difficulty here is that a truncated Edgeworth expansion never provides us with a well defined non-Gaussian PDF (normalised to one and positive) as soon as the non-Gaussianity is non zero and even if tiny (this is due to the fact that a cumulant generating function is a polynomial only in the case of a Gaussian random variable). However, we will show in the next section that large deviations theory can be used to predict from first principles this PDF without suffering from the issues of negativity and normalisation.

In this framework, we will see that all cumulants are accounted for and are exact at tree order.

3 Count-in-cells as modelled by Large deviation theory Large deviation theory is at the core of many developments in physics, in particular in statistical physics, in order to describe the tails of exponentially suppressed distributions with the idea that "an unlikely fluctuation is brought about by the least unlikely among all unlikely paths". As such it encompasses and improves upon the central limit theorem which only describes behaviours close to the maximum of the PDF (likely configurations). This theory has recently been applied in cosmology to study the time evolution of the cosmic density field on large scales (Bernardeau & Reimberg 2016). A precise formulation of large deviation theory is beyond the scope of this proceeding and we therefore chose to only describe the qualitative picture here. Further mathematical details can be found in e.g Uhlemann et al. (2016).

The goal here is to predict the statistics of densities in spheres of a given size. In the Gaussian initial conditions, this statistics is Gaussian and well known. However in the late time Universe and on intermediate scales, the field is not Gaussian anymore. As illustrated on Fig. 2, the Lagrangian patches (i.e where the matter is coming from in the initial conditions) corresponding to spheres of same size R and density ρ, show a large variety of shapes, densities and volumes (only the mass is conserved from the initial to the final conditions). This is expected since at a given macroscopic configuration (mean density ρ in the sphere) correspond many different microscopic states that originate from very different initial configurations. In other words, many different paths can connect the initial conditions to the same late time mean density in a sphere. Amongst all this possible paths, there is however one path that dominates the statistics even more so that the event is rare and therefore the variance σ 2 of the field is small. Given that we consider spheres that are maximally symmetric, one can conjecture that the most likely path should also respect this symmetry and may therefore be the so-called spherical collapse model. The spherical collapse mapping is a well known exact non-linear solution of the gravitational evolution when the initial condition is an isolated spherically symmetric fluctuation. Explicit solutions can be found for specific cosmologies (e.g Einstein de Sitter) while in general numerical solutions are necessary. However, it was shown that once expressed in terms of the linear density, the spherical collapse mapping is very weakly dependent on cosmology and is well reproduced by the simple analytic form ζ SC (δ l ) = (1 -δ l /ν) -ν with ν = 21/13. Once the conjecture on the most likely path is made and in the SF2A 2018 limit of small variance, large deviation theory can be used to predict the cumulant generating function of the field and eventually get the resulting PDF which is nothing but the inverse Laplace transform of the cumulant generating function. Eventually, the explicit prediction for the density PDF can be obtained (see e.g equation of Codis et al. (2016b) and the corresponding public code LSSFast distributed freely at http: //cita.utoronto.ca/ ~codis/LSSFast.html). All cumulants are predicted at once by this formalism and are exact at tree order in PT. Thanks to this property, the predicted PDF matches very well the measurements in simulation for variance of order unity as can be seen in Fig. 3, even in the tails that are governed by rare (non-linear) events which are beyond the reach of PT techniques. In addition, all predictions are analytical with an explicit dependence on cosmology and no free parameters, a unique situation in this field of research. Note that this formalism can be implemented not only for one sphere but for an arbitrary number of concentric spheres. Large-deviation theory has been successfully applied to the density one-point distribution but also to the density slopes or profiles (Bernardeau et al. 2015), cosmic velocities, to two-point statistics that allowed us to model the error budget expected from finite volumes probed by galaxy surveys (Codis et al. 2016a). Primordial non-Gaussianities were also investigated in this framework. It was shown that using the full knowledge of the PDF, tighter constraints on cosmology (in particular the equation of state of dark energy) could be obtained (Codis et al. 2016b) and once tracer bias is accounted for (we observe galaxies not the underlying total matter density field), density PDFs could help breaking the degeneracies between biasing (nuisance) parameters and cosmological ones (Uhlemann et al. 2018).

Conclusions

The distribution of cosmic fields cannot be modelled in the tails using perturbative techniques. A large deviation principle can however be implemented and provides us with a detailed and accurate modelling even in the tails of the distributions and in the mildly non-linear regime beyond what is usually achievable via PT, reaching subpercent precision for scales above 10 Mpc/h at redshift zero. These ideas are promising to extract cosmological information from the non-Gaussian large-scale structure of the Universe as observed by future galaxy surveys (Euclid among others).

Fig. 1 .

 1 Fig. 1. Distribution of density contrast in spheres of 10Mpc/h at various redshifts from blue to red. The dashed lines are the corresponding second order Edgeworth expansion.

Fig. 2 .

 2 Fig. 2. Spheres of the same radius and density at late times (right) can originate from very different patches in the initial conditions (left), with various sizes, shapes and densities.

Fig. 3 .

 3 Fig. 3. Distribution of density measured in the Horizon-Run 4 simulation (with error bars) at z = 0.7 and for various radii as labelled and compared to the large deviations theory predictions (solid).

  SC thanks Fondation Merac for funding. This work and SC's attendance to the 2018 SF2A meeting were also partially supported by the Programme National Cosmologie et Galaxies (PNCG) of CNRS/INSU with INP and IN2P3, co-funded by CEA and CNES.

CNRS and Sorbonne Université, UMR 7095, Institut d'Astrophysique de Paris, 98 bis Boulevard Arago, F-75014 Paris, France c Société Francaise d'Astronomie et d'Astrophysique (SF2A)