Keywords: fast multipole method, fast Fourier transform, random projection, dislocation dynamics, covariance matrix, multidimensional scaling . . . . . . . . . . . . isotropic elastic fields . . . . . . . . . . 2.2.1 Isotropic elastic stress through Mura's formula . . . . . . . . . . 2.2.2 Isotropic elastic nodal forces and energy . . . . . . . . . . . . . . 2.3 Existing Fast Multipole DD formulations . . . . . . . . . . . . . .

Advanced techniques for the low-rank approximation of matrices are crucial dimension reduction tools in many domains of modern scientific computing. Hierarchical approaches like H 2 -matrices, in particular the Fast Multipole Method (FMM), benefit from the block low-rank structure of certain matrices to reduce the cost of computing n-body problems to O(n) operations instead of O(n 2 ). In order to better deal with kernels of various kinds, kernel independent FMM formulations have recently arisen such as polynomial interpolation based FMM. However, they are hardly tractable to high dimensional tensorial kernels, therefore we designed a new highly efficient interpolation based FMM, called the Uniform FMM, and implemented it in the parallel library ScalFMM. The method relies on an equispaced interpolation grid and the Fast Fourier Transform (FFT). Performance and accuracy were compared with the Chebyshev interpolation based FMM. Numerical experiments on artificial benchmarks showed that the loss of accuracy induced by the interpolation scheme was largely compensated by the FFT optimization. First of all, we extended both interpolation based FMM to the computation of the isotropic elastic fields involved in Dislocation Dynamics (DD) simulations. Second of all, we used our new FMM algorithm to accelerate a rank-r Randomized SVD and thus efficiently generate multivariate Gaussian random variables on large heterogeneous grids in O(n) operations. Finally, we designed a new efficient dimensionality reduction algorithm based on dense random projection in order to investigate new ways of characterizing the biodiversity, namely from a geometric point of view.

I would like to acknowledge the HiePACS team for providing such a peaceful and inspiring atmosphere of research during my stay at Inria Bordeaux. A huge thanks to the ScalFMM group, in particular Bérenger Bramas, who has always been understanding and has always made himself available for nice enriching discussions. I am very grateful to my team mates as well, namely Jean-Marie, Arnaud, Maria, Yuval, Grégoire, Abdou, Manu, with whom I shared wonderful times. Thanks Thomas Traub and Martin Schanz for there support before and during the thesis. My experience at TU Graz definitely had a huge impact on this work. Thank you Franz, Michaël, Matthias and Bernhard as well! I wish to thank my supervisor Olivier Coulaud for allowing me to address such a wide variety of scientific issues. Thank you Olivier for being always so supportive and patient. Thank you Eric Darve for always suggesting novel ideas and for raising my knowledge to randomized numerical linear algebra and many other great topics in scientific computing. Thank you for hosting me at Stanford University and allowing to meet great people there.

Thank you Laurent Dupuy for introducing me to the world of dislocations and the beautiful community of multiscale material modeling. Finally, many thanks to Alain Franc for the wonderful exchanges we had around randomized algorithms, geometry and biodiversity. Thank you especially for trusting me in the amazing project we started building.

To everyone thank you for those small attentions that made this research experience so delightful, I wish success to all of you!

Résumé

Les techniques avancées pour l'approximation de rang faible des matrices sont des outils de réduction de dimension fondamentaux pour un grand nombre de domaines du calcul scientifique. Les approches hiérarchiques comme les matrices H 2 , en particulier la méthode multipôle rapide (FMM), bénéficient de la structure de rang faible par bloc de certaines matrices pour réduire le coût de calcul de problèmes d'interactions à n-corps en O(n) opérations au lieu de O(n 2 ). Afin de mieux traiter des noyaux d'interaction complexes de plusieurs natures, des formulations FMM dites "kernel-independent" ont récemment vu le jour, telles que les FMM basées sur l'interpolation polynomiale. Cependant elles deviennent très coûteuses pour les noyaux tensoriels à fortes dimensions, c'est pourquoi nous avons développé une nouvelle formulation FMM efficace basée sur l'interpolation polynomiale, appelée Uniform FMM. Cette méthode a été implémentée dans la bibliothèque parallèle ScalFMM et repose sur une grille d'interpolation régulière et la transformée de Fourier rapide (FFT). Ses performances et sa précision ont été comparées à celles de la FMM par interpolation de Chebyshev. Des simulations numériques sur des cas tests artificiels ont montré que la perte de précision induite par le schéma d'interpolation était largement compensées par le gain de performance apporté par la FFT. Dans un premier temps, nous avons étendu les FMM basées sur grille de Chebyshev et sur grille régulière au calcul des champs élastiques isotropes mis en jeu dans des simulations de Dynamique des Dislocations (DD). Dans un second temps, nous avons utilisé notre nouvelle FMM pour accélérer une factorisation SVD de rang r par projection aléatoire et ainsi permettre de générer efficacement des champs Gaussiens aléatoires sur de grandes grilles hétérogènes. Pour finir, nous avons développé un algorithme de réduction de dimension basé sur la projection aléatoire dense afin d'étudier de nouvelles façons de caractériser la biodiversité, à savoir d'un point de vue géométrique.

Mots clés: méthode multipole rapide, transformé de Fourier rapide, projection aléatoire, dynamique des dislocations, matrices de covariance, positionnement multidimensionnel.

Résumé étendu

Introduction Les techniques avancées pour l'approximation de rang faible des matrices sont des outils d'algèbre linéaire numérique cruciaux pour la réduction de dimension. Elles interviennent dans un nombre important de domaines du calcul scientifique, dont certains sont abordés dans cette thèse. Des variantes hiérarchiques comme les matrices dites H 2 , en particulier la méthode multipôle rapide (FMM), peuvent en plus bénéficier de la structure de rang faible par bloc de certains opérateurs matriciels. Pendant les deux dernières décennies, la FMM a été utilisée de manière intensive dans de nombreuses applications allant de la simulation de problèmes à n-corps à la résolution d'équations aux dérivées partielles (EDP). De manière générale, la FMM est une méthode hiérarchique d'algèbre linéaire numérique permettant d'accélérer les calculs mettant en oeuvre des matrices noyaux (ou matrices d'interactions), comme défini ci-après K = {k(∥x i -x j ∥ 2 )} i,j=1...n ou {x i } i=1...n est une distribution spatiale de points et k est une fonction noyaux (fonction de 2 variables d'espace) suffisamment régulière mais pouvant présenter une singularité à l'origine. La FMM permet de réduire le coût de calcul des problèmes d'interactions à ncorps ou d'application de produits matrice-vecteur à O(n) opérations au lieu de O(n 2 ). D'autre part, de plus en plus d'applications nécessitent de mettre en oeuvre des noyaux d'interaction très complexes, en particulier de natures différentes et de dimensions différentes (potentiellement grandes). Afin de gérer ce vaste ensemble de fonctions, des formulations FMM dites indépendantes du noyaux (ou kernel-independent) ont récemment vu le jour. Par exemple, les formulations multipôles rapides basées sur l'interpolation polynomiale reposent sur l'interpolation du noyau d'interaction en p + 1 points d'une grille notée {x a } a=1.. Cependant, de telles approches s'étendent difficilement à des noyaux d'interactions tensoriels lorsque leur dimension est élevée, en particulier parce que l'approximation de rang faible est sous optimale. Nous nous restreignons dans cette thèse au cas des noyaux non-oscillants, mais mentionnons dans ce manuscrit les récentes avancées qui permettent d'étendre notre approche aux noyaux oscillants. L'avantage de cette approche par rapport aux autres méthodes dites kernel-independent est qu'elle ne repose sur aucune hypothèse concernant la nature du noyau, en particulier son lien éventuel avec une EDP. D'autre part, son implémentation est relativement simple car elle repose sur des outils mathématiques relativement basiques et peut bénéficier des outils optimisés pour l'algèbre linéaire numérique dense tels que blas et lapack.

Fast interpolation based FMM Dans cette thèse, nous avons développé une nouvelle formulation FMM hautement performante basée sur l'interpolation polynomiale, appelée Uniform FMM ou ufmm, qui permet de réduire le coût de calcul de l'étape la plus coûteuse de l'algorithme FMM, à savoir l'étape de transfert M2L. Cette technique repose simplement sur une reformulation de l'opérateur M2L sans utiliser d'hypothèse supplémentaire sur le noyau et sans introduire une nouvelle approximation matricielle. Cette approche a été implémentée dans la bibliothèque ScalFMM * , où elle bénéficie des dernières évolutions en terme de vectorisation, de paradigmes pour le calcul parallèle et de moteur d'exécutions. Notre méthode repose sur une grille d'interpolation régulière, ce qui implique une structure dites Toeplitz à l'opérateur M2L et permet de réduire son coût de stockage. De plus, la conversion de l'étape M2L dans le domaine de Fourier par transformée de Fourier rapide (FFT) permet de réduire significativement la complexité de l'étape de transfert M2L. Les performances et la précision de l'approche ont été comparé aux méthodes existantes telles que la FMM basée sur l'interpolation de Chebyshev et ses variantes optimisées par recompression algébrique des opérateurs. Des simulations numériques sur des cas tests artificiels ont montré que la perte de précision introduite par le schéma d'interpolation était largement compensé par l'accélération FFT. D'autre part, l'empreinte mémoire de l'étape M2L a été réduite de manière significative au prix de stocker des développements multipolaires de taille plus importante.

Dynamique des dislocations Dans la première partie de cette thèse, nous avons étendue les formulations FMM par interpolation à l'évaluation des champs isotropiques élastiques mis en oeuvre dans les simulations de Dynamique des Dislocations (DD). Le modèle DD repose sur un représentation intégrale aux frontières du champs de force élastique créé par un ensemble de défauts linéiques en interactions élastiques, ces défauts sont aussi connus sous le nom de dislocations. La simulation d'ensemble massifs de dislocations, par exemple pour des discrétisations allant jusqu'à O(10 7 ) segments, permet de comprendre de manière plus précise les lois de comportement plastique des matériaux cristallins irradiés. Entreprendre de telles simulations a nécessité d'adapter les routines d'interpolation de la bibliothèque ScalFMM aux éléments intégraux et à la nature tensorielle des interactions propres aux dislocations. En utilisant ScalFMM comme moteur d'exécution de la FMM, nous avons étendu les capacités du code DD parallèle OptiDis † en fournissant une évaluation multipôle rapide performante des champs élastiques isotropiques. Notre méthode a permis de réaliser des simulations DD avec n = O(10 7 ) segments en prenant en compte les contributions de champs lointain de la contrainte, de la force et de l'énergie élastique, tout en maintenant une précision de l'ordre de 10 -4 et minimisant le coût de ce champs lointain. D'autre part l'utilisation d'une approche kernel-independent a permis de facilement mettre en oeuvre des formulations optimisées du champs d'interaction pour son évaluation multipôle rapide, et facilitera la mise en oeuvre de noyaux anisotropiques. D'autre part, l'interpolation polynomiale a permis de diminuer la dimensionnalité du noyaux et d'évaluer efficacement les opérateurs intégraux. Enfin, nos travaux ont permis de mener à bien des simulations clés dans le domaine des dislocations dans un temps raisonnable et avec la précision attendue, en particulier nous avons pu simuler les mécanismes associés à l'apparition de bandes claires dans les alliages de Zirconium (Figure 1).

Matrices de covariance

Dans la seconde partie de cette thèse, nous avons mis en oeuvre notre nouvel algorithme FMM afin de générer de manière efficace des variables aléatoires Gaussiennes multidimensionnelles sur des grilles arbitraires et en particulier très hétérogènes. L'approche standard consiste à appliquer une racine carrée de la matrice de covariance, calculée à partir de la donnée d'une fonction de corrélation et d'une grille de n points, à un bruit blanc généré sur cette grille. Cette méthode est particulièrement coûteuse, en effet O(n 3 ) opérations sont nécessaire pour calculer de manière exacte une racine carrée de matrice par une factorisation de Cholesky ou SVD, de plus il est nécessaire d'assembler la matrice pour en extraire sa racine. L'approche que nous proposons repose sur l'accélération du calcul de la racine en combinant deux techniques avancées d'approximation de rang faible:

• une SVD de rang-r approchée par projection aléatoire, aussi connu sous le nom de randomized SVD,

• accélérée par des multiplications de matrices approchées grâce à la FMM.

Le coût de la randomized SVD étant gouverné par 2r produits matrice-vecteur, la FMM permet de réduire le coût de pré-calcul de la racine carré de rang faible de O(rn 2 ) à O(r 2 n) opérations sans qu'il n'y ait besoin de stocker la matrice de covariance. Premièrement, nous avons optimisé notre implémentation ufmm afin de mieux gérer un nombre élevé de vecteur d'entrée, à savoir r = O( 103 ). Dans un second temps, étant donné la régularité de certaines fonctions de corrélation utilisées dans ce type d'applications, e.g., la corrélation Gaussienne, nous avons développé un algorithme d'approximation de rang faible par blocs basé sur l'interpolation polynomiale, en changeant simplement le critère d'admissibilité de la version standard de la FMM par interpolation. Notre approche hiérarchique a permis de générer des champs aléatoires Gaussiens sur des grilles allant jusqu'à n = O (10 6 ) points en quelques secondes. Pour finir, nous avons illustré que pour des distributions de points très hétérogènes, le caractère hiérarchique de notre approche permettait d'augmenter ces performances de manière significative, en effet le remplissage de l'arbre est d'autant plus creux que la grille est hétérogène.

Positionnement multidimensionel Dans la dernière partie de la thèse, nous avons mis au point un algorithme de réduction de dimension efficace pour étudier de nouvelles façons de décrire la biodiversité, à savoir d'un point de vue géométrique. Notre approche est basée sur une technique de visualisation appelée le positionnement multidimensionnel ou multidimensional scaling (MDS). Cette technique permet de représenter un ensemble de n points dans un espace euclidien à r dimensions avec r << n grâce à la seule donnée des distances deux-à-deux entre ces points. Le coût de calcul de la MDS est en général largement dominé par l'étape de factorisation de la matrice de similarité associée à la matrice de distance, en effet les méthodes standards requièrent O(n 3 ) opérations et nécessitent l'assemblage et le stockage de la matrice. La plupart des méthodes rapides pour le positionnement multidimensionnel sont basées sur des techniques de factorisation par sélection aléatoire de colonnes (ou échantillonnage aléatoire), qui ne s'appliquent qu'à des matrices définie-positives, à savoir la méthode de Nyström. D'autre part, une discussion sur l'influence de la présence de valeurs propres négatives sur la distorsion du nuage a également été menée.

xiii Notre nouvelle approche met en oeuvre un algorithme de SVD basé sur la projection Gaussienne dense connu sous le nom de randomized SVD. La projection aléatoire a de multiples avantages sur la sélection de colonnes aléatoire dans ce domaine bien que cette dernière soit moins coûteuse en terme de calcul. Premièrement, la projection aléatoire est par essence plus proche du concept de réduction de dimension, car une illustration plus directe du le lemme de Johnson-Lindenstrauss. Ensuite, elle permet une factorisation sous la forme SVD de manière rapide, stable et précise, en fournissant des informations cruciales sur la structure de la matrice, à savoir une approximation de l'image et du rang numérique. D'autre part, la randomized SVD présentent des avantages significatifs en terme d'implémentation et de performance. La projection aléatoire a permis de faire tourner l'algorithme de positionnement multidimensionnel sur de grands échantillons fournis par la nouvelle génération de technologies de séquençage d'ADN dites ngs, et ce grâce à une implémentation C++ efficace de la randomized SVD (mais aussi sur la décomposition de Nyström) dans la bibliothèque fmr ‡ .

Des méthodes de visualisations basiques ont été utilisées pour représenter les nuages de points dans un nombre relativement faible de dimensions. Nous avons pu traiter de manière fluide et précise des échantillons jusqu'à O (10 5 ) séquences (ou reads) provenant du Lac Léman, et ainsi illustrer deux phénomènes:

• la réunion en cluster des reads associés à des espèces identifiées • et la forte concentration des reads dans des régions très localisées et bien séparées entre elles.

Une étude plus avancée de ce dernier phénomène devrait permettre la définition de nouveaux Operational Taxonomic Units (OTUs) et ainsi une caractérisation de la biodiversité sur un critère purement géométrique.

Conclusion

Nous avons montré sur deux applications l'importance de développer des algorithmes efficaces qui exploitent la structure de rang faible par bloc des matrices d'interaction et prennent en compte l'hétérogénéité du domaine. Une application en physique des matériaux nous a permis de souligner la nécessité de développer des formulations peu coûteuses et générique pour traiter des interactions très complexes. Une application en géostatistique nous a permis d'illustrer les bénéfices d'une représentation hiérarchiques de la matrice de covariance sur une technique efficace pour la factorisation globale des matrices (ou plus précisément dans ce cas le calcul de racines carrés), reposant sur la projection aléatoire. Enfin, une approche plus exploratoire nous a permis d'étendre la factorisation par projection aléatoire aux matrices de similarités mises en oeuvre pour la classification de grands jeux de données en biodiversité. Cette dernière a fourni des résultats prometteurs qui doivent permettre une réduction supplémentaire du coût de calcul en exploitant plus finement la structure des matrices.

Introduction Introduction Introduction

Challenges in Scientific Computing Applied physicists, biologists, geoscientists, climatologists, economists and many other kind of computational scientists keep developing new models in order to simulate real-life phenomenon more accurately and at more realistic scales. However, because of their complex nature, large dimensions and extreme refinement, these models are more and more demanding in terms of computational resources. Therefore, the challenges faced by the scientific computing community for the next few decades is huge, in particular when it comes to tackling the curse of dimensionality induced by the use of high dimensional models or projection spaces. For instance, new approaches developed by researchers in mechanical engineering and materials sciences in order to better understand the behavior of complex structures or materials involve a growing number of parameters, which leads to very intensive computations and often requires advanced HPC solutions and, what is called in this domain, reduced order models [START_REF] Zahm | Model order reduction methods for parameter-dependent equations-Applications in Uncertainty Quantification[END_REF]. On the other hand, Molecular Dynamics (MD) and Multiscale Materials Modeling (MMM) as a whole have pushed the barriers of scientific computing by addressing arbitrary high order models for the interactions at a given scale and also by coupling computations at various scales [START_REF] Elliott | Novel approaches to multiscale modelling in materials science[END_REF][START_REF] Díaz De La | Materials research by means of multiscale computer simulation[END_REF]. A large part of scientific computing is dedicated to the resolution of problems formulated in term of n-body interactions or relying on the discretization of continuous operators on a mesh. Therefore, they usually implement relatively basic matrix computations such as multiplication, inverse, square root or more evolved matrix functions. Over the last two decades many algorithms have already been designed to reduce the complexity of these fundamentals operations. A significant part of these algorithms belong to the framework of sparse Numerical Linear Algebra (NLA) as it makes use of the sparsity of the matrix operators, but they will not be discussed in this thesis. On the other hand, the framework of dense NLA addresses the computation of more general dense matrices, which rapidly involve prohibitive resources consumptions. In particular, since the cost of applying such operators is quadratic in the problem size and factorization in low rank form as well as inversion are cubic, the challenge in this field has been the design of linear scaling algorithm. Among the notable algorithms that contributed most to the field of scientific computing, we would like to emphasize two kinds of algorithms. First of all, hierarchical algorithms such as the Fast Multipole Method (FMM) and H matrices as they allow for fast matrix multiplications as well as fast inversion. Then, randomized algorithms have recently brought important contributions to the fast factorization of matrices in standard forms that are oblivious of the matrix sparsity.

Hierarchical Algorithms

The computation of n-body interactions problems lead to the design of the very popular FMM [START_REF] Greengard | A fast algorithm for particle simulations[END_REF], in order to apply discrete matrix operators at a linear cost. Although the FMM terminology is highly oriented towards physical applications, it has proven its efficiency and relevance in many other areas of computational sciences. In particular, the FMM became very popular in the mechanical engineering community for the resolution of boundary value problems, namely by accelerating the matrix multiplications involved in the iterative solvers [START_REF] Nishimura | Fast multipole accelerated boundary integral equation methods[END_REF][START_REF] Tausch | The variable order fast multipole method for boundary integral equations of the second kind[END_REF][START_REF] Chaillat | Fast Multipole Method for 3-D elastodynamic boundary integral equations[END_REF][START_REF] Messner | Fast directional multilevel summation for oscillatory kernels based on chebyshev interpolation[END_REF][START_REF] Sylvand | La méthode multipôle rapide en électromagnétisme : Performances, parallélisation, applications[END_REF]. In fact, it probably even motivated the design of its algebraic counterpart, namely the more general framework of hierarchical matrices, a.k.a., H-matrices [START_REF] Bebendorf | Hierarchical matrices: a means to efficiently solve elliptic boundary value problems[END_REF][START_REF] Hackbusch | H2-matrix approximation of integral operators by interpolation[END_REF]. Both methods benefit from the block low rank structure of certain discrete or integral operators and rely on a data sparse representation of the matrices. Furthermore, they involve approximations that need to be set w.r.t. to the target application. They have not yet been integrated into generic packages such as lapack, because their implementation need a very special care in order to achieve optimal performance. However, the formulation and implementation of the FMM and H-matrices have become more and more generic in order to easily adapt to any application [START_REF] Ying | A kernel-independent adaptive fast multipole algorithm in two and three dimensions[END_REF][START_REF] Fong | The black-box fast multipole method[END_REF]. Nonetheless, FMM algorithms still need to be improved in order to cope with the growing demand in high dimensional models and extremely refined simulations.

Randomized Algorithms

In the last few years, a new domain of NLA, called randomized NLA, has gained popularity thanks to its potential benefits to data analysis applications, e.g., machine learning or data mining. This sudden interest for randomized algorithms lead to highly efficient methods for the Low-Rank Approximation (LRA) of matrices with application to a wide range of domains [START_REF] Michael | Randomized algorithms for matrices and data[END_REF]. Randomized LRA are random in the sense that they involve sketching of a matrix using random procedures [START_REF] David | Sketching as a tool for numerical linear algebra[END_REF], but most importantly because they rely on error bounds that hold with probability deriving from the theory of the concentration of measure [START_REF] Achlioptas | Database-friendly random projections: Johnson-lindenstrauss with binary coins[END_REF]. All randomized algorithms more or less derive from the beautiful idea, illustrated by Johnson-Lindenstrauss lemma [START_REF] William | Extensions of lipschitz mappings into a hilbert space[END_REF], that for any given set of point in an Euclidean space there exists an arbitrary low-distortion embedding into a lower dimensional Euclidean space. Despite relying on a relatively involved mathematical background [START_REF] Indyk | Approximate nearest neighbors: towards removing the curse of dimensionality[END_REF], the general idea of randomized LRA algorithms remain very intuitive and their implementation straightforward. Moreover, as they rely on very basic matrix computations [START_REF] Halko | Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions[END_REF], they can further benefit from the structure and the sparsity of the input matrix and they parallelize well. In this thesis we will focus on dense random projection based approaches [START_REF] Vempala | The Random Projection Method[END_REF] and there low-rank variants [START_REF] Liberty | Accelerated dense random projections[END_REF].

Scope of study

First of all, our work aim at proposing new variants of the FMM that are better suited for the requirements of modern scientific applications, such as Multiscale Materials Modeling through Dislocation Dynamics (DD) simulations. Second of all, we want to highlight the interest of using the hierarchical algorithms in a larger variety of domains, in particular on geostatistical applications, that involve large grids and non-oscillatory interaction kernels of various nature. Furthermore, we want to benefit from this framework to illustrate the benefits of combining hierarchical algorithms with randomized LRA in terms of performance. In this work, we focus on developing efficient sequential algorithms that fully benefit from the recent advances in NLA and that can easily be ported on parallel applications. Although no HPC solution is actually proposed in this thesis, all our contributions are developed in commonly used and maintained open-source parallel libraries that implement the state-of-the-art parallel paradigms.

Overview

The present manuscript is organized in 3 main parts, focusing on various aspects of dense NLA for computationally demanding modern scientific applications. In particular, it describes the implementation of new efficient hierarchical algorithms, randomized algorithms or combination of both.

Preliminary Part I introduces the various algorithms studied in this thesis and the current state-of-the-art. In particular, Chapter 1 describes the main concepts and features of the FMM and interpolation based variants. Then, Chapter 2 presents the framework of our first application, namely Dislocation Dynamics. Finally, Chapter 3 introduces the fundamental concepts of randomized NLA for the low-rank approximation of dense matrices.

n-body interactions Part II gathers our various contributions to the computation of nbody interaction problems. First of all, Chapter 5 describes a highly efficient interpolation based FMM that applies to any non-oscillatory kernel. Then, Chapter 6 implements this algorithm along with other interpolation schemes for computing isotropic elastic interactions in large dislocation networks.

Covariance matrix Part III focuses on our contributions to covariance matrix computations. In Chapter 7 we describe a linear time algorithm for performing randomized LRA of matrices using our interpolation based FMM for performing fast matrix multiplications with application to the generation of Gaussian Random Fields (GRFs). Finally, Chapter 8 addresses the design of a geometric view on biodiversity by mean of random projection aided dimensionality reduction.

Part I

Scientific background & Related works 1.1.1 Introduction

The simulation of n-body problems remains a crucial issue in a wide variety of scientific domains involving elastic, acoustic, electric, gravitational interactions,…These problems usually boil down to evaluating sums of the form

ϕ(x i ) = n ∑ j=1 k(x i , x j )w(x j ) (1.1)
where x ∈ R n×d denote the points of a d-dimensional grid, i.e., the set of n bodies/particles represented in d dimensions and w(x i ), respectively ϕ(x i ), denotes the density, respectively the potential, associated to the particle i located at x i . Originally introduced by Greengard and Rokhlin [START_REF] Greengard | A fast algorithm for particle simulations[END_REF], the FMM provides an efficient method for approximating ϕ(x i ), i.e., the potential created by all particles on the particle i, for all i ∈ {1, . . . , n} up to a given accuracy

ε F M M in O(n) operations instead of the O(n 2
) operations required by the naïve direct computation of the sums. From an algebraic point of view, it is equivalent to performing an approximate product between a kernel matrix and a vector, where the input and output vectors are respectively w = {w(x i )} i=1,...,n and ϕ = {ϕ(x i )} i=1,...,n and the kernel matrix reads as

K = {k(x i , x j )} i,j∈{1,...,n} .
The FMM relies on 2 main ingredients:

• the hierachical partitioning of the input grid using a binary tree (e.g., quadtree in 2D, octree in 3D),

• the separation or decoupling of variables x and y in k(x, y), e.g., using analytic expansions of the kernel k.

If the input kernel matrix has a block low-rank structure, then hierachical partitioning of the grid allows for efficiently building those blocks. On the other hand, analytical expansions are used for computing the smooth block-to-block interactions in low-rank form.

The Fast Multipole Algorithm

The Fast Multipole Algorithm (FMA) is a 2-stage algorithm that produces approximate potentials φ given a set of particles x, a kernel k(•, •) and densities w using a cluster tree, an admissibility criterion and a low-rank approximation scheme. Let us describe the algorithm in full details in this subsection.

Octree First of all, the FMA relies on a hierarchical partitioning of the domain using a cluster tree structure. Since only 3D grids (or spatial grids) are considered in this thesis, the name octree will be used to denote the tree structure. The root cluster of the tree is the smallest cube enclosing all particles. At the level L of the octree we subdivide all parent cells in 8 cubes of equal size to get all child cells at level L + 1. This recursive partition is stopped at the leaf level L once a suitable criterion is reached, e.g., minimum or average number of particles per leaf cell, minimum leaf cell size. We denote

C (L) x
the level-L cell containing the particle located at x. Its center and its width are denoted respectively c (L) x and ω (L) . In particular, ω (L)

x = ω (0) /2 L with ω (0) the width of the root cluster. Finally, the number of particles contained in

C (L) x is denoted n x .
Admissibility Once the domain is partitioned into clusters we need to define admissible cluster pairs T

(L) xy = (C (L) x , C (L) y )
. At a given level L and for a fixed target cell C (L)

x , the pair

T (L)
xy is uniquely defined by its transfer vector

t (L) xy = c (L) x -c (L) y ω (L) If T (L)
xy is admissible, then the matrix of interactions

K(T (L) xy ) = { k(x i , y j ) } i<nx,j<ny ∈ R nx×ny
can be approximated at a given precision using a low-rank method, i.e.,

K(T (L) xy ) ≈ U(C (L) x )C(T (L) xy )V T (C (L) y ) with U ∈ R nx×r , C ∈ R r×r , V ∈ R ny×r and r = r xy << n x , n y . A pair of clusters T (L)
xy is said admissible if it satisfies ∥t (L) xy ∥ ≥ γ (admissibility criterion)

where γ is a prescribed parameter. While the original FMA [START_REF] Greengard | A fast algorithm for particle simulations[END_REF] enforces γ = 2, i.e., the well-separation of cells or strong admissibility, general H 2 -matrices can also enforce γ = 1, i.e., the adjacency of cells, a.k.a., weak admissibility, e.g., the HSS structure [START_REF] Xia | Fast algorithms for hierarchically semiseparable matrices[END_REF][START_REF] Chandrasekaran | Some fast algorithms for hierarchically semiseparable matrices[END_REF]. See Bebendorf [START_REF] Bebendorf | Hierarchical matrices: a means to efficiently solve elliptic boundary value problems[END_REF] for a comprehensive introduction to the concept of H-matrix.

Nesting A particular aspect of the FMM that produces a H 2 -format is the nesting of the basis between levels. In fact, the basis U resp. V at a given level can be expressed in terms of the basis at a lower level, resp. upper level, by a relation of the following kind

U iα (C (L-1) x ) = U iα ′ (C (L) x )U ′ α ′ α (C (L) x ) V T βi (C (L-1) x ) = V ′ T ββ ′ (C (L) x )V T β ′ i (C (L) x )
1.1. The Fast Multipole Method for all i = 1, . . . , n x , where U ′ , V ′ ∈ R r×r denote the operators used to transfer the basis, and ultimately the expansions, from one level to another.

Neighbor and Interaction lists Let C (L)

x denote a target cell at level L. Cells that form admissible pairs with C

I(C (L) x ) = { C (L) y ∈ tree : T (L)
xy is admissible and

C (L-1) y ∈ N (C (L) x )

}

The latter condition ensures that all contributions are only computed once during the hierarchical summation scheme. For γ = 2, I(C

(L)
x ) contains at most 6 d -3 d cells for any L, i.e., 189 in 3D. Fig. 1.2 shows the neighbor list (colored in red) and interaction list (colored in green) on a portion of a quadtree for a given target cell C (L)

x . On the other hand, two brother cells do not share the same interaction list, therefore there are 7 d -3 d possible configurations of interacting pairs, i.e., 316 in 3D. Given a fictitious target cell C (L)

x , we denote the set of all possible interactions I f ull (C

(L) x ). ( L -1) ( L) ( L) x y a y b y c C ( L) x C ( L-1) x C ( L) y a C ( L) y b C ( L) y c C ( L-1) y c P2P P2M M2M M2L L2L L2P P2M M2L

Fast Multipole Algorithm

The standard FMA computes approximate particle interactions in a hierarchical fashion, i.e., considering multiple levels of the tree. The multi-level fast multipole summation is presented in algorithm 1. It consists in 2 stages, namely 2 passes through the binary tree: an upward pass followed by a downward pass. Figure 1.1 illustrates the interaction and neighbor lists on the last two levels of a simple 1D tree structure. It also schematically shows how the FMA transfers nearfield and farfield contributions by mean of X2Y operations, where X and Y can be equal to • P, when contributions are expressed on the source or target particles contained in leaf cells,

• M, when contributions are expressed as multipole expansions M in source cells,

• L, when contributions are expressed as local expansions L in target cells.

For a general d-dimensional grid, the FMA consists in successively performing the following operations:

• The Upward Pass consists in performing the following operations at each level from bottom (L = L) to top (L = 2) -P2M: Compute all multipole expansions (M) at the leaf level, i.e., accumulate contributions from source particles y to their containing leaf cell

C ( L) y in M(C ( L) y ) M β (C ( L) y ) = ny ∑ j=1 V T βj (C ( L) y )w(y j ), ∀β ≤ r (1.2) -M2M: If L < L, transfer the multipole expansions from 2 d level-(L + 1) child cells to their level-(L) parent cell M β (C (L) y ) = ∑ β ′ ≤r V T ββ ′ (C (L) y )M β ′ (C (L+1) y ), ∀β ≤ r (1.3)
• The Downward Pass consists in performing the following operations at each level from top (L = 2) to bottom (L = L) -M2L: Transfer multipole expansions (M) from a maximum of

6 d -3 d (189 in 3D) source cells C (L) y ∈ I(C (L) x ) into local expansions (L) of target cell C (L) x L α (C (L) x ) = ∑ C (L) y ∈I(C (L) x ) ∑ β≤r C αβ (T (L) xy )M β (C (L) y ), ∀α ≤ r (1.4) -L2L: If L < L,
transfer all the local expansions from level-(L) parent cells to their

2 d level-(L + 1) child cells L α (C (L+1) x )+ = ∑ α ′ ≤r U αα ′ (C (L) x )L α ′ (C (L) x ), ∀α ≤ r (1.5)
-L2P: Aggregate local expansions (L) at the leaf level on the target particles x for all target cells

C ( L) x φ(x i ) = ∑ α≤r U iα (C ( L) x )L α (C ( L) x ) (1.6)
M β (C ( L) y ) = ny ∑ j=1 V T βj (C ( L) y )w(y j ), ∀β ≤ r else // M2M: from 2 d childs to parent M β (C (L) y ) = ∑ β ′ ≤r V T ββ ′ (C (L) y )M β ′ (C (L+1) y ), ∀β ≤ r // Downward pass for level L = 2, . . . , L do for target cell C (L)
x ∈ tree do // M2L: from Multipole to Local expansions

for source cell C (L) y ∈ I(C (L) x ) do L α (C (L) x )+ = ∑ β≤p C αβ (T (L) xy )M β (C (L) y ), ∀α ≤ r if L < L then // L2L: from parent to 2 d childs L α (C (L+1) x )+ = ∑ α ′ ≤r U αα ′ (C (L) x )L α ′ (C (L) x ), ∀α ≤ r else // L2P: from leaf Local Expansion to target Particles φ(x i ) = ∑ α≤r U iα (C ( L) x )L α (C ( L) x )
// P2P: direct computation of the nearfield

for source cell C ( L) y ∈ N (C ( L) x ) do φ(x i )+ = ny ∑ j=1 k(x i , y j )w(y j )

Computational cost and asymptotic complexity

Although the linear asymptotic complexity of the FMM and H 2 -methods (its algebraic counterpart) is well-established, in this subsection we recall the computational cost of the FMA in details regardless of the expansion technique (series expansion, polynomial interpolation,…) and regardless of the interaction kernel k. Therefore, we denote • t k , the time required to compute k(x, y) for one couple (x, y) multiply by the density w(y) and add to the potential ϕ(x)

• t U/V T = O(r), the time required to compute a multipole expansion for a given particle

• t C = O(r 2
), the time required to compute the basis of the expansion U or V T at a given level

Consequently, the time required to transfer an expansion from child to parent or the other way around reads as r × t U/V T = O(r 2 ). For the sake of simplicity, we consider the case, where all cells at all levels are non-empty and the number of particles per leaf is a constant Nearfield interactions Each cell is in nearfield interaction with at most 27 cells (including itself) in 3D. At the leaf level, nearfield interactions are actually computed in a direct fashion, therefore the time t P2P required to compute the P2P step is equal to the time required to compute interactions between n 0 source and 27n 0 target particles time the number of leaf cells:

n 0 = n/2
t P2P = ∑ C ( L) x ∈tree (n 0 × 27n 0 )t k = 2 d L × n 0 × 27n 0 × t k = n × 27n 0 × t k
Hence, the P2P step has a linear complexity in n, i.e., t P2P = O(n). In order to maintain this linear complexity, n 0 has to remain constant or at least bounded as n grows. For a given geometry, if particles are uniformly distributed (i.e., homogeneous density of particles) then increasing the tree depth by one level every time the problem size is multiplied by 2 d ensures a constant n 0 .

Farfield interactions On the other hand, since the size of the interaction list is bounded for a given level and the number of cells scales like 2 d L = n/n 0 then the cost of the farfield is linear in n. Indeed, the asymptotic cost of each step of the FMM can be estimated as follows:

• P2M/L2P t P2M/L2P = (2 d L × n 0 ) × t U/V T = n × t U/V T = O(n) • M2M/L2L t M2M/L2L = ( L-1 ∑ L=2 2 dL × 2 d ) × (r × t U/V T ) ≤ (n/n 0 × 2 d ) × (r × t transf er ) = O(n) • M2L t M2L = L ∑ L=2 2 dL × 189 × t C ≤ n/n 0 × 189 × t C = O(n)
Balance The fundamental idea behind FMM is to replace direct computations by lowrank approximations while using a binary tree and interaction lists in order to ensure an asymptotic linear complexity. The common procedure to setup an FMM algorithm is to experimentally determine the order of the expansion in order to setup the accuracy. At this point the cost of computing interactions between cells in farfield or nearfield interactions is fixed but not necessarily well-balanced. Increasing the number of level will increase the size of the farfield, i.e., the number of interactions that are computed in a low-rank fashion, and decrease the number of nearfield interactions, i.e., direct computations. Therefore, the optimal tree depth of any FMM algorithm is usually obtained when nearfield and farfield computations are perfectly balanced.

Interpolation based FMM

The original FMA as introduced by Greengard and Rokhlin [START_REF] Greengard | A fast algorithm for particle simulations[END_REF] relies on expansions of the kernel k(r) = 1/r in terms of spherical harmonics, i.e., in the basis of Legendre polynomials expressed in spherical coordinates. Since then many algorithms have been derived from this one using the same terminology but implementing other strategies for the expansion of the kernel. They are usually divided in 2 classes, whether the expansion technique depends on the kernel or not.

• A kernel dependent approach implementing -Spherical Harmonics expansions for k(r) = 1/r was introduced in the original paper [START_REF] Greengard | A fast algorithm for particle simulations[END_REF].

-Taylor series expansions for Dislocation Dynamics (elastostatics) was introduced by Arsenlis et al. [8].

-Hankel series expansions for frequency domain elastodynamics k(r) = e ikr /r was introduced by Chaillat et al. [START_REF] Chaillat | Fast Multipole Method for 3-D elastodynamic boundary integral equations[END_REF].

• Kernel independent approaches like -the kifmm, introduced by Ying et al. [START_REF] Ying | A kernel-independent adaptive fast multipole algorithm in two and three dimensions[END_REF], relies on Boundary integral representation.

-the skeletonization-based FMM proposed by Martinsson and Rokhlin [START_REF] Martinsson | An accelerated kernel-independent fast multipole method in one dimension[END_REF].

-the bbfmm, introduced by Fong and Darve [START_REF] Fong | The black-box fast multipole method[END_REF], relies on Chebyshev polynomial interpolation.

In this thesis we focus on polynomial interpolation based variants like the bbfmm. First of all, we present the general framework of polynomial interpolation based FMM. Then, we discuss some useful features of polynomial interpolation. Finally, we present the bbfmm algorithm using a polynomial basis independent terminology.

Chebyshev polynomial interpolation

As discussed in the thesis of Matthias Messner [START_REF] Messner | Fast Boundary Element Method in Acoustics[END_REF] or Michael Messner [89], the Chebyshev polynomial interpolation is usually a good alternative to the Chebyshev series expansion. For instance, as shown here, the computation of the coefficient c n takes the basic form of a discrete summation over p + 1 Chebyshev nodes instead of an integral over [-1, 1].

Basics

The most simple way of constructing a polynomial approximation of order p of a continuous function f on [-1, 1] is to interpolate f at p + 1 equidistant points {x i } i=0,...,p by

f (x) ≈ f p (x) = p ∑ n=0 c n x n
The interpolation procedure leads to the resolution of a system of p + 1 equations. This method is time consuming and becomes numerically unstable as p grows. Moreover, since the set of interpolation points is chosen equally spaced, the interpolating polynomials does not converge uniformly on [-1, 1] (Runge phenomenon). A better choice for {x i } i=0,...,p is the set of Chebyshev points, i.e., the root of the Chebyshev polynomials of the first kind

T p+1 , namely xm = cos ( 2m + 1 p + 1 π 2
)

, for m = 0, . . . , p

The resulting interpolating polynomials converge exponentially on [-1, 1] for analytic functions. Moreover, there is a quasi-best approximation result in ∞-norm, namely

∥f -f p ∥ ∞ ≈ max x∈[a,b] |f (x) -f p (x)|,
that involves the Lebesgue constants, see Rivlin [START_REF] Tj Rivlin | The lebesgue constants for polynomial interpolation[END_REF].

Lebesgue constants

In order to quantify the quality of an interpolation technique, i.e., the choice of interpolation nodes, we often use the Lebesgue constants. These quantities grow only logarithmically with p when interpolating at Chebyshev nodes, whereas they grow exponentially when interpolating on equispaced nodes.

Chebyshev interpolation formula

Let us write the polynomial f p of order p interpolating function f on Chebyshev nodes as a sum of Chebyshev polynomials of the first kind

f (x) ≈ f p (x) = p ∑ n=0 c ′ n T n (x) = 1 2 c 0 + p ∑ n=1 c n T n (x), ∀x ∈ [-1, 1] (1.8)
then the coefficients c n read as

c n = 2 p + 1 p ∑ m=0 f (x m )T n (x m )
As shown in [START_REF] Fong | The black-box fast multipole method[END_REF] and recalled in [START_REF] Messner | Fast Boundary Element Method in Acoustics[END_REF] exponential convergence is also ensured for functions f (•) and kernels k(•, •) when using Chebyshev polynomial interpolation.

A generic interpolation formula

Let us now re-write the interpolation formula (1.8) in a more practical and generic form

f (x) ≈ p ∑ m=0 S(p; x, xm )f (x m ), ∀x ∈ [-1, 1]
(1.9)

where the interpolation operator S is defined as

S(p; x, xm ) = 1 p + 1 + 2 p + 1 p ∑ n=1 T n (x)T n (x m ) (1.10)
The interpolation formula (1.9) is generic in the sense that changing the interpolation polynomials (or nodes) only affects the definition of S. In the following developments we will always consider (1.9) when doing polynomial interpolation. Furthermore, for the sake of clarity we will omit the interpolation order p in the definition of the order-p interpolation operator, i.e.,

S(x, xm ) ≡ S(p; x, xm )

Useful features of polynomial interpolation

Interpolating on an arbitrary interval In the general case, we may need to interpolate a function f on an arbitrary [a, b]. This can be done provided the mapping function Φ defined as follows

                 Φ : [-1, 1] → [a, b] x → Φ(x) = a + b 2 + b -a 2 x Φ -1 : [a, b] → [-1, 1] x → Φ -1 (x) = 2x -b -a b -a (1.11)
The interpolation formula of order p then reads as

f (x) ≈ p ∑ m=0 S(Φ -1 (x), xm )f (Φ(x m )), ∀x ∈ [a, b] , i.e., f (x) ≈ p ∑ m=0 S( 2x -b -a b -a , xm )f ( a + b 2 + b -a 2 xm ), ∀x ∈ [a, b]
The multi-index notation Let f denote a function mapping from R 3 to C, a similar interpolation formula as (1.10) can be derived for f that reads as

f (x) ≈ ∑ |α|≤p S(x, xα )f (x α ), ∀x ∈ [-1, 1] 3
This formula uses the third-dimensional multi-index notation defined as follows α := (α 1 , α 2 , α 3 ), where ∀i ∈ {1, 2, 3} , α i = 0, . . . , p.

The summation over all combinations of α 1 , α 2 and α 3 can be denoted

∑ |α|≤p ≡ p ∑ α 1 =0 p ∑ α 2 =0 p ∑ α 3 =0
where

|α| := max i≤3 (α i )
Please note that the interpolation points are represented in a tensorial way such that

xα := (x α 1 , xα 2 , xα 3 ) ∈ [-1, 1] 3
Using this notation, the interpolation operators can be expressed in a shorter form

S(x, xα ) := S(x 1 , xα 1 )S(x 2 , xα 2 )S(x 3 , xα 3 )
which makes them easier to read and manipulate.

Shifting derivatives Polynomial interpolation techniques allow for approximating the derivatives of f at almost no extra cost by simply deriving the interpolation formula. This feature was described and mathematically analyzed in [START_REF] Börm | Hybrid cross approximation of integral operators[END_REF][START_REF] Börm | Approximation of integral operators by variable-order interpolation[END_REF] and relies on a stability result, namely a generalization of Markov's inequality. Additionally, since derivation decreases the order of the polynomial basis by one, this technique provides an order p -1 polynomial approximation of the derivative of f from the interpolation of f at p + 1 points. Let f denote a univariate analytic function defined on [-1, 1], then for all x ∈ [-1, 1] we have

f ′ (x) ≈ (f p (x)) ′ (1.12)
This result does not hold in general but holds for analytic functions. Therefore, incorporating the expression of f p given in (1.9) into the (1.12) yields the following approximation of

f ′ f ′ (x) ≈ p ∑ m=0 ∂S ∂x (x, xm )f (x m ), ∀x ∈ [-1, 1]
Hence, approximating f ′ boils down to shifting the derivative to the interpolating polynomial S in (1.9). This trick is particularly convenient, since S ′ is known explicitly for most interpolation polynomials. For a function f defined on an arbitrary interval [a, b] the approximation reads as

f ′ (x) ≈ p ∑ m=0 ∂Φ -1 ∂x (x) ∂S ∂x (Φ -1 (x), xm )f (Φ(x m )), ∀x ∈ [a, b] (1.13)
If we replace Φ by its definition (1.11), the approximation formula (1.13) becomes

f ′ (x) ≈ p ∑ m=0 2 b -a ∂S ∂x ( 2x -b -a b -a , xm )f ( a + b 2 + b -a 2 xm ), ∀x ∈ [a, b]
Let f denote a function defined on a 3-dimensional interval [-1, 1] 3 , we have

∇ x f (x) ≈ ∑ |α|≤p ∇ x S(x, xα )f (x α ), ∀x ∈ [-1, 1] d
where the gradient is defined as

∂S ∂x 1 (x, xα ) := S ′ (x 1 , xα 1 )S(x 2 , xα 2 )S(x 3 , xα 3 ) ∂S ∂x 2 (x, xα ) := S(x 1 , xα 1 )S ′ (x 2 , xα 2 )S(x 3 , xα 3 ) ∂S ∂x 3 (x, xα ) := S(x 1 , xα 1 )S(x 2 , xα 2 )S ′ (x 3 , xα 3 )
Applications of this approach to the FMM can be found in [START_REF] Messner | Fast Boundary Element Method in Acoustics[END_REF][START_REF] Traub | A kernel interpolation based fast multipole method for elastodynamic problems[END_REF]. It will often be used in this work in order to decrease the dimensionality of the kernel during the M2L step at the cost of decreasing the accuracy of the approximation scheme.

The bbfmm, a Chebyshev interpolation based FMM

In this subsection, we present the interpolation formula of an arbitrary non-oscillatory kernel function and the fast multipole summation scheme associated to it. The resulting algorithm is an interpolation based FMM oblivious to the polynomial basis, that coïncides with the bbfmm for the basis of Chebyshev polynomials. Most fundamental mathematical results behind this method were covered 7 years prior the bbfmm in the framework of H 2 matrices, see [START_REF] Hackbusch | H2-matrix approximation of integral operators by interpolation[END_REF].

Interpolating interaction kernels

y ∈ I(C (L)
x ), the interaction k(x, y) can be approximated within well-controlled error bounds (see [START_REF] Messner | Fast Boundary Element Method in Acoustics[END_REF] for further details). As a result, the farfield contributions to the potential created at position x i

ϕ f ar (x i ) = ∑ y j ∈I(C (L) x i ) k(x i , y j )w(y j )
can be approximated at a fixed precision by the following expression

ϕ f ar (x i ) ≈ ∑ |α|≤p S(x i , xα ) ∑ |β|≤p k(x α , ȳβ ) ∑ y j ∈I(C (L) x i )
S(y j , ȳβ )w(y j ) using a low-rank representation of the interaction kernel, where U iα (C

(L) x ) = S(x i , xα ), V T βj (C (L) y ) = S(y j , ȳβ ) and C αβ (T (L) xy ) = k(x α , ȳβ
). An important feature of polynomial interpolation is that it only requires the evaluation of k on a fixed grid, which ensures both kernel independence and modularity.

Compact notations

For the sake of clarity we sometimes simplify notations by omitting the dependence in the grid points, i.e., we denote

S α (x) = S(x, xα ), S β (y) = S(y, ȳβ ), Kαβ = k(x α , ȳβ ) (1.15)
Thus, the compact interpolation formula reads as

k(x, y) ≈ ∑ |α|≤p S α (x) ∑ |β|≤p Kαβ S β (y) (1.16)
Fast summation scheme The Fast Multipole summation scheme associated with polynomial interpolation is presented in Algorithm 2, while Algorithm 3 defines the functions that are specific to the original bbfmm as introduced by Fong and Darve [START_REF] Fong | The black-box fast multipole method[END_REF]. Here we recall the equations associated with the main steps of the algorithm.

• The Upward Pass consists in performing the following operations at each level from bottom (L = L) to top (L = 2) -P2M: Aggregate all source contributions from the particles to the multipole expansion at the leaf level

M β (C ( L) y ) = ny ∑ j=1 S(y j , ȳβ )w(y j ), ∀β/|β| ≤ p (1.17)
-M2M: Transfer multipole expansions child to parent cells

M β (C (L) y ) = ∑ |β ′ |≤p S(ȳ β , ȳβ ′ )M β ′ (C (L+1) y ), ∀β/|β| ≤ p (1.18)
• The Downward Pass consists in performing the following operations at each level from top (L = 2) to bottom (L = L) -M2L: Transfer all multipole expansions between all admissible cluster pairs

L α (C (L) x )+ = ∑ |β|≤p Kαβ M β (C (L) y ), ∀α/|α| ≤ p (1.19)
-L2L: Transfer local expansions from parent to child cells

L α (C (L) x ) = ∑ |α ′ |≤p S(x α , xα ′ )L α ′ (C (L-1) x ), ∀α/|α| ≤ p (1.20) -L2P: Aggregate all leaf-level local expansions L(C ( L)
x ) on the target particles x

φ(x i )+ = ∑ |α|≤p S(x i , xα )L α (C ( L) x ) (1
M β (C (L) y ) = ∑ |β ′ |≤p S(ȳ β , ȳβ ′ )M β ′ (C (L+1) y ), ∀β/|β| ≤ p M(C (L) y ) ←transformM(M(C (L) y )) // Downward pass for level L = 2, . . . , L do for target cell C (L)
x ∈ tree do // M2L: transfer expansions between interacting cells applyM2L( L(C

(L) x ), C (L) x ) L(C (L) x ) ←transformL( L(C (L) x )) if L < L then // L2L: interpolation from parent to 2 d childs L α ′ (C (L+1) x ) = ∑ |α ′ |≤p S(x α , xα ′ )L α ′ (C (L)
x ), ∀α/|α| ≤ p else // L2P: interpolation from leaf cell to target particles

φ(x i )+ = ∑ |α|≤p S(x i , xα )L α (C ( L) x )
// P2P: direct computation of the nearfield

for source cell C (L) y ∈ N (C (L) x ) do φ(x i )+ = ∑ y j ∈C (L) y k(x i , y j )w(y j ) 1
y ∈ I f ull (C (L) x ) do Ki αβ (L) = k(x α , ȳβ ), ∀α, β/|α|, |β| ≤ p Function applyM2L( L(C (L) x ), C (L) x ) for i-th source cell C (L) y ∈ I(C (L) x ) do Lα (C (L) x )+ = ∑ |β|≤p Ki αβ (L) Mβ (C (L) y ), ∀α/|α| ≤ p Function transformM(M) return M Function transformL( L) return L

Computational cost and memory footprint

The computational cost and memory footprint of the M2L step is much larger than the other steps in most FMM variants, this subsection shows that this also applies to the interpolation based FMM. In order to accurately describe the cost of interpolation based FMM such as the bbfmm, Table 1. 6 ) memory footprint of these operators represents the largest part of the memory requirements. As a result, computations involving large tree depth and p > 10 exceed the limits of standard desktop computers in term of storage. On the other hand, the cost of applying K to the multipole expansions scales like O(p 6 ). Moreover, these matrices must be applied at all levels and a potentially large number of times (max. 189 interactions per cell per level), which usually makes the M2L step the most computationally expensive step of any Fast Multipole scheme.

Homogeneity

In the case of a homogeneous kernel k, i.e.,

∃α ∈ R/∀(λ, x, y) ∈ R × R d × R d , k(λx, λy) = |λ| α k(x, y), (1.23
L(C (L) x )+ = λ(α, L) Kref M(C (L) y ) (1.24)
This optimization results in a reduction of the memory footprint by a factor of L -1 = 7 for the largest tree depth considered in this thesis, namely L = 8, but also reduces the precomputation time.

Related works on M2L-optimized FMM

The main contribution of this thesis, described in Section 5.1, consists in a new M2Loptimized FMM powered by Fast Fourier Transform (FFT). While optimization of the interpolation based M2L operators and conversion of the FMM to Fourier domain have already been investigated independently in the past, they have never been combined. Here, we discuss those past developments in order to show how the method positions itself in this context.

M2L-optimization of the bbfmm

As investigated in [START_REF] Fong | The black-box fast multipole method[END_REF], the rank of the interpolation-based M2L operators is not optimal, therefore compression of K should be considered. Messner et al. [START_REF] Messner | Optimized m2l kernels for the chebyshev interpolation based fast multipole method[END_REF] describes several optimizations of the bbfmm based on compressing the M2L operators. These methods are briefly described here and will be used as references in our numerical benchmarks (see Section 5.4).

Compressed variant

The first variant called the compressed-bbfmm and presented in Algo 4 relies on the global compression of the set of M2L operators. The resulting factorization of the M2L allows for shifting part of its application to the P2M/L2P step and store the result in transformed expansions. This method usually requires a significant amount of precomputation time and does not always provide optimal compression and application time as opposed to the second variant.

U L , [ C 1 L , . . . , C 316 L ] , V L ←SVD [ K1 (L), . . . , K316 (L) ] Function applyM2L( L(C (L) x ), C (L) x ) for i-th source cell C (L) y ∈ I(C (L) x ) do L α (C (L) x )+ = ∑ |β|≤p (C i L ) αβ M β (C (L) y ), ∀α/|α| ≤ p Function transformM(M) // Apply V T L . M = V T L M Function transformL( L) // Apply U L . L = U L L(, ∀α/|α| ≤ p
Symmetric variant For kernels exhibiting symmetries, like k(x, y) = k(y, x), the symmetricbbfmm presented in Algo 5 allows for significantly reducing the memory footprint of the bbfmm while reducing both precomputation and running time. In the symmetric-bbfmm, the symmetries of the kernel and the interpolation grid allow for reducing the size of the interaction list from |I f ull (C

(L) x )| = 316 to |I sym (C (L)
x )| = 16 for a given level L. Consequently, we need to store less M2L operators and they can be applied to a set of multipole expansions at once, however expansions have to be permuted properly before the M2L translation. On the other hand, individual compression of the operators becomes affordable, which results in a near-optimal low-rank compression and thus a faster application of each M2L operator. In particular, a fully or partially pivoted Adative Cross Approximation (ACA) is used to compress the M2L operators with or without assembling it. Then, recompression can be performed using a QR Decomposition and a subsequent Singular Value Decomposition (SVD). Although this method is kernel dependent it applies to a wide range of kernels, e.g., RBF kernels and kernel satisfying k(x, y) = k(|x -y|) in general.

Algorithm 5: symmetric-bbfmm functions

Function precomputeM2L() // For all levels in tree, or 1 level if kernel is homogeneous for L = 2, . . . , L do // Assemble all possible interactions, i.e., 16

for i-th source cell C (L) y ∈ I sym (C (L) x ) do U i L , V i L ←ACA [ Ki (L) ] Function applyM2L( L(C (L) x ), C (L) x ) for i-th source cell C (L) y ∈ I(C (L) x ) do L α (C (L) x )+ = ∑ |β|≤p ∑ r (U i L ) αr (V i L ) rβ M β (C (L) y ), ∀α/|α| ≤ p

Conversion of a FMM to Fourier domain

Principle An FFT conversion of the original Fast Multipole Algorithm was originally proposed in a technical report by Greengard et al. [START_REF] Greengard | On the efficient implementation of the fast multipole algorithm[END_REF], and later implemented and improved by Elliott et al. [START_REF] William | Fast fourier transform accelerated fast multipole algorithm[END_REF]. For an efficient parallel implementation of this approach please refer to Pan et al. [START_REF] Pan | The power of combining the techniques of algebraic and numerical computing: Improved approximate multipoint polynomial evaluation and improved multipole algorithms[END_REF]. The general idea of the original FFT conversion proposed in these papers consists in rewriting the M2L operations in the form of a 2D linear convolution, then converting it to a 2D circular convolution (using zero-padding) in order to finally perform the M2L operations in Fourier domain at a lower cost using FFT.

Limits This approach differs from our method since it is based on expanding the kernel 1/r in spherical harmonics while our method relies on kernel-independent polynomial interpolations. On the other hand, both variants require zero-padding for the conversion to circular convolution. The FFT conversion of the FMA algorithm additionally suffers from smearing effects, i.e., some entries are artificially not zeroed out due to the conversion to circular convolution form. While our method suffers from the Runge phenomenon due to the interpolation on an equispaced grid, the original variant suffers from various sources of instability. The first occurs when a coefficient with strong variations in magnitude is warped in order to rewrite the expansion as a linear convolution. The second instability results from a scaling problem that amplifies during the FFT operations. Stabilization techniques were proposed in [START_REF] Greengard | On the efficient implementation of the fast multipole algorithm[END_REF] and further discussed and improved in [START_REF] William | Fast fourier transform accelerated fast multipole algorithm[END_REF]. In Section 5.1, we also suggest methods in order to defeat the Runge phenomenon that may occur in our method.

Introduction to Dislocation Dynamics

Developing a deeper understanding of material's behaviour is crucial for the design of new structures. For instance, providing fine models for the plastic behaviour of complex alloys submitted to high gradients of temperature and high levels of radiations has become a key issue in nuclear safety.

A multiscale approach The general multiscale approach to material modeling is represented on Figure 2.1, it can be summarized in 3 main scales:

• Molecular Dynamics (MD) represents the material at the atomistic scale. MD simulations allows to identify elastic laws by analyzing atomistics interactions on very large sets of atoms.

• Dislocation Dynamics (DD) bridges the gap between the MD scale and the FEM scale by modeling the moving plans of atoms represented in terms of line defects called dislocations. DD can be used to simulate either a single grain or a large number of grains. Ultimately, DD produces empirical plasticity laws by simulating the propagation of massive ensembles of dislocations.

• Finite Element Methods (FEM) are the well-known and broadly used methods that aim at modeling elastoplastic or more complex media at the mesoscale. The internal stress or strain states generated by FEM is usually used for design purposes. The DD model Dislocations are by definition the result of an eigenstrain state in a medium, see Mura's book [START_REF] Mura | Micro-mechanics of Defects in Solids[END_REF]. They are usually studied in the context of lattice defects in crystalline materials. In particular, the simulation of massive ensembles of dislocations allows for a better understanding of the mesoscopic plasticity laws. Several models have emerged in order to make such simulations computationally tractable to large ensembles.

In this thesis we focus on the most widely spread approach, namely Dislocation Dynamics (DD). Here, we will only recall the fundamentals of the nodal approach to DD simulations, please refer to Bulatov and Cai [START_REF] Bulatov | Computer simulations of dislocations[END_REF] for an extensive introduction to the topic or see Kubin [START_REF] Kubin | Dislocations, mesoscale simulations and plastic flow[END_REF] for a more physics-oriented approach.

Topology Dislocations consist in closed line defects that are ubiquitous in crystalline materials. As shown on Figure 2.2, a dislocation line (in red) is characterized by a Burger's vector b representing the intensity and direction of the dislocation. At any point x of the line we also define the line tangent t(x). In the nodal approach to DD, dislocations lines are discretized by segments and nodes. In particular, segments are treated as straight dislocation lines with constant line tangent vectors t, and all segments discretizing a dislocation loop share the same Burger's vector b as the loop. Finally, the segment length denoted L is given in Å(short for ångström, 1Å= 0.1nm= 10 -10 m), typically 20Å≤ L ≤ 100Å. In this thesis, all lengths and distances involved in DD are given in Å.

b

x i t x j

Dynamics and algorithm

The motion of dislocations is usually initiated by an external elastic stress state applied to the material. This stress states acts on dislocation lines along with the internal stress state created by the dislocation network itself. The resulting forces applied on the lines create the motion according to advanced mobility laws. While moving, dislocations can collide, merge, annihilate or even multiply, creating a potentially significant modification of the network's topology. Please refer to Arnaud Etcheverry's thesis [START_REF] Etcheverry | Simulation de la dynamique des dislocations à très grande échelle[END_REF] for an extensive presentation of the algorithm and the efficient implementation in OptiDis code * .

Integral representations of the isotropic elastic fields

In this section we present expressions for the isotropic elastic stress, forces and energy fields in integral form. The elastic stress field is first represented in integral form using Mura's formula. Then, we recall the definition of the elastic nodal forces acting on the extremities of finite dislocation lines. Finally, we address the integral representation of the elastic energy.

Isotropic elastic stress through Mura's formula

Mura's formula The elastic stress field created at point x by a dislocation loop (C ′ ) with Burger's vector b ′ is given in [START_REF] Mura | Micro-mechanics of Defects in Solids[END_REF] by Mura's formula as

σ ij (x) = C ijkl (C ′ ) ε lnh C pqmn G kp,q (|x -x ′ |)b ′ m dx ′ h , (2.1)
where G ij is the Green's function for the infinite elastic media characterized by the Hooke's tensor C. We denote dx ′ (x ′ ) the infinitesimal variation of x ′ and t ′ (x ′ ) the line tangent vector at point

x ′ . Thus, if dx ′ (x ′ ) = |dx ′ | then dx ′ = t ′ (x ′ )dx ′ (x ′ ).
Finally, ε denotes the permutation symbol and we omit the summation over repeated indices using Einstein's summation convention. A demonstration of Mura's formula can be found in Mura's book [START_REF] Mura | Micro-mechanics of Defects in Solids[END_REF].

Isotropic elastic stress field

In isotropic elasticity, the Hooke's tensor reduces to the 2 Lamé constants (µ, ν) and the Green's function reads as

G ij (r) = 1 8πµ (δ ij r ,pp - 1 2(1 -ν) r ,ij ), (2.2)
where r = xx ′ and r = |r|. The notation r ,i denotes the partial derivatives of r w.r.t. x i , i.e., r ,i = ∂r ∂x i .

Replacing the expression (2.2) of G in (2.1) yields

σ ij ((C ′ ), x) = µb ′ k 8π (C ′ ) r ,mpp (ε jmk dx ′ i + ε imk dx ′ j ) + 2 1 -ν ε nmk (r ,ijm -δ ij r ,ppm )dx ′ n (2.3)
Let us write (2.3) in a slightly different form and distinguish 2 main contributions

σ ij ((C ′ ), x) = µ 8π (σ A ij + σ A ji )((C ′ ), x) + 2 1 -ν (σ B ij -δ ij σ B pp )((C ′ ), x)) (2.4)
where σ A and σ B are defined as follows

σ A ij ((C ′ ), x) = (C ′ ) r ,mpp (x -x ′ )ε jmk b ′ k t ′ i dx ′ σ B ij ((C ′ ), x) = (C ′ ) r ,ijm (x -x ′ )ε nmk b ′ k t ′ n dx ′
Once σ A (x) and σ B (x) are computed one can easily deduce the expression of σ(x) using (2.4). In the isotropic case σ is symmetric. More precisely, both terms σ A ij + σ A ji and (σ B ij -δ ij σ B pp ) are symmetric. Hence, we can always compute σ ij for i ≤ j and deduce the remaining components by transposition. and the target segment

x 1 x 2 x 3 x 4 x x ′ r t ′ t L ′ L
[ x 3 x 4 ] .
Finite dislocation line Although (2.3) is only defined for a closed dislocation loop, not only is it useful but also necessary to write this expression for a segment, i.e., a finite dislocation line, in order to be able to sum up all contributions of a given discretized dislocation loop. Let us consider a source segment [x 1 x 2 ], i.e., extremity 1 (resp. 2) is located at point x 1 (resp. x 2 ), characterized by a Burger's vector b ′ , a line tangent t ′ and a length L ′ . Thus, the elastic stress created by the finite dislocation line [x 1 x 2 ] reads as

σ ij ( [ x 1 x 2 ] , x) = µ 8π (σ A ij + σ A ji )( [ x 1 x 2 ] , x) + 2 1 -ν (σ B ij -δ ij σ B pp )( [ x 1 x 2 ] , x)
where

σ A ([x 1 x 2 ] , x) and σ B ([x 1 x 2 ] , x) are defined as follows σ A ij ( [ x 1 x 2 ] , x) = ∫ x 2 x 1 r ,mpp (x -x ′ )ε jmk b ′ k t ′ i dx ′ = ∫ L ′ 0 r ,mpp (x -x ′ )ε jmk b ′ k t ′ i dξ ′ σ B ij ( [ x 1 x 2 ] , x) = ∫ x 2 x 1 r ,ijm (x -x ′ )ε nmk b ′ k t ′ n dx ′ = ∫ L ′ 0 r ,ijm (x -x ′ )ε nmk b ′ k t ′ n dξ ′ where ξ ′ ∈ [0, L ′ ] denotes the curvilinear abscissae of x ′ , x ′ = x 1 + ξ ′ t ′ and dx ′ = dξ ′ .
Non-singular formalism The integrals involved in (2.3) actually exist, however they cannot be evaluated numerically using standard quadratures because of the singularity occuring when r approaches 0, i.e., x ′ → x. Attempts at regularizing the integrals have been introduced in Brown [START_REF] Brown | The self-stress of dislocations and the shape of extended nodes[END_REF], that gave the first definition of the core, then Peierls [START_REF] Peierls | The size of a dislocation[END_REF] and Nabarro [START_REF] Frn Nabarro | Dislocations in a simple cubic lattice[END_REF] developed the idea of spreading the Burger's vector around the dislocation line. In this thesis focus is made on the non-singular formulation introduced by Cai et al. [START_REF] Cai | A non-singular continuum theory of dislocations[END_REF]. This method roughly consists in removing the singularity by considering a slightly different Green's function, i.e.,

G ns ij (r a ) = 1 8πµ (δ ij r a,pp - 1 2(1 -ν) r a,ij ),
where the distance is shifted by a small and constant contribution, namely a parameter a called the core width, such that

r a = √ a 2 + r 2 .
The newly defined Green's function G ns is non-singular since for any

(x, x ′ ) ∈ (R 3 ) 2 , r a (x, x ′ ) ̸ = 0.
Not only does this formulation imply very little modifications in terms of implementation compared to the singular variant, but it also relates to the atomistic models * . Cai et al. [START_REF] Cai | A non-singular continuum theory of dislocations[END_REF] also provides analytical expressions for the isotropic elastic stress field and energy using this formalism. These expressions were initially implemented in OptiDis code for the direct computation of the isotropic elastic stress, forces and energy.

Isotropic elastic nodal forces and energy

Sources of motion In the DD model, dislocation lines move under the effect of elastic forces acting on nodes, that divide in 2 contributions: the internal nodal forces f = f int. and the external nodal forces f ext. . While the former is induced by the elastic stress field σ created by dislocations themselves, the latter is the result of an external stress field σ ext. enforced on the boundary of the grain.

The Peach-Koehler force and nodal forces The local force induced by an elastic stress state σ on any point x of a dislocation loop (b, t, L) is called the Peach-Koehler force and it reads as

f P K (x) = (σ(x) • b) × t
On the other hand, the elastic nodal forces f i acting at the extremity i = 3, 4 of a finite dislocation line [x 3 x 4 ] characterized by (b, t, L) are obtained by integration of the Peach-Koehler force on [x 3 x 4 ], i.e.,

f i = ∫ x 4 x 3 f P K (x)N i (x)dx
where N i=3,4 are piecewise linear shape functions, such that N i (x j ) = δ ij . If we denote ξ the curvilinear abscissae of x, then x = x 3 + ξt, dx = dξ and the linear shape functions N i=3,4 read as

N 4 (x) = |x -x 3 | L = ξ L and N 3 (x) = |x 4 -x| L = L -ξ L = 1 - ξ L (2.5)
Hence, the force acting on node 4 and node 3 reads as

f 4 = ∫ L 0 f P K (x) ξ L dξ f 3 = ∫ L 0 f P K (x)dξ -f 4
External forces The external stress state is usually enforced at the boundary of the grain by a simple Dirichlet boundary condition. Therefore, it is considered constant within the grain (far from the boundary). The external nodal forces are computed by integration of the Peach-Koehler force associated to the external stress state σ ext. over the target segment

[x 3 x 4 ], i.e., f ext.,4 = ∫ L 0 (σ ext. • b) × t ξ L dξ = L 2 (σ ext. • b) × t
Hence, external forces are computed directly at a linear cost in n.

Internal forces Similarly, the internal nodal force is computed by integration of the Peach-Koehler force associated to the internal stress state σ int. = σ over the target segment

[x 3 x 4 ], i.e., f 4 i = ∫ L 0 ε ijk σ jp (x)b p t k ξ L dξ, ∀i ∈ {1, 2, 3} (2.6)
Therefore, the non-singular contribution of σ A to the internal forces can be expressed as the following double line integral

f A,4 i = ∫ L 0 ε ijk b p t k ξ L dξ ∫ L ′ 0 r a,mrr (x -x ′ )ε pms b ′ s t ′ j dξ ′ , ∀i ∈ {1, 2, 3} (2.7)
and similarly for the contribution of σ B . Hence, in a naive implementation, the cost of computing the internal forces over the entire network is quadratic in n.

Elastic energy Let us consider the following expression of the isotropic elastic energy of a dislocation network (C) taken from [START_REF] Cai | A non-singular continuum theory of dislocations[END_REF]

E((C)) = - µ 8π (E A + 2 1 -ν (E B -E C + νE D )) (2.8)
where contributions E A , E B , E C and E D are defined as

E A = (C) b i dx i (C) r a,kk b ′ j dx ′ j , E B = (C) b j dx k (C) r a,ij b ′ i dx ′ k E C = (C) b i dx j (C) r a,kk b ′ i dx ′ j , E D = (C) b j dx i (C) r a,kk b ′ i dx ′ j
This expression originates from De Wit et al. [START_REF] De | Some relations for straight dislocations[END_REF] along with various other forms. All existing forms are equivalent up to an integral, that vanishes on a closed loop.

Existing Fast Multipole DD formulations

Since the interactions between dislocations are long-ranged and fast-decreasing (see the expression of the Green's function (2.2)), a common approach is to use the Fast Multipole Method (FMM) in order to reduce the asymptotic cost of computing the elastic fields of a n-segments network from a quadratic to a linear cost in n. While nearfield interactions are evaluated directly using analytical expressions, farfield interactions are approximated in order to be evaluated more efficiently.

Direct evaluation of the nearfield As mentioned by Arsenlis et al. [8] the size of the quadrature required to evaluate (2.6) up to a fixed accuracy grows dramatically as the distance between dislocation lines decreases. Therefore, it is not recommended to evaluate the double line integral involved in the nodal forces numerically. Instead, [8] relies on the non-singular theory [START_REF] Cai | A non-singular continuum theory of dislocations[END_REF], i.e., σ = σ ns , in order to derive analytical expressions of the force field created by a source segment and acting on a parallel or a non-parallel segment, as well as on the segment itself. For the computation of the nearfield contributions, OptiDis implements its own analytical expressions that are derived from the non-singular theory [START_REF] Cai | A non-singular continuum theory of dislocations[END_REF]. For equivalent analytical expressions, please refer to [8] for the stress and force fields (2.7), or [START_REF] Cai | A non-singular continuum theory of dislocations[END_REF] for the energy (2.8).

Approximation of the farfield Several formulations have already been considered for the approximation of the farfield contributions to the elastic stress and energy, such as expansions in spherical harmonics [START_REF] Zhao | A new version fast multipole method for evaluating the stress field of dislocation ensembles[END_REF] or Taylor expansions [START_REF] Wang | Multipole representation of the elastic field of dislocation ensembles[END_REF]. The well-known DD code ParaDis developed at Lawrence Livermore National Laboratory first implemented Taylor expansions of the isotropic elastic fields [8], and only recently incorporated polynomial interpolation based on the bbfmm. The aim of our work was to implement an efficient interpolation scheme into OptiDis, the parallel implementation of NumoDis code (CEA Saclay), that did not involve farfield approximation yet.

3

Randomized low-rank approximation of matrices

The Randomized SVD

Introduction to randomized algorithms

Randomized algorithms for the low-rank approximation of matrices, a.k.a., randomized Low-Rank Approximation (LRA) algorithms, usually provide powerful alternatives to standard matrix factorizations for matrices of relatively low-rank such as the Singular Value Decomposition (SVD), QR Decomposition, Cholesky Decomposition, … Recently, they have gained popularity in the numerical linear algebra community because they are easy to implement, highly parallelizable and most of all they achieve very competitive performance within reasonable accuracy. Moreover, these algorithms often involve very basic matrix computations thus leaving significant room for improvement. They have drawn much attention in scientific fields such as geostatistics [START_REF] Pk Kitanidis | Principal component geostatistical approach for largedimensional inverse problems[END_REF][START_REF] Dehdari | Applications of randomized methods for decomposing and simulating from large covariance matrices[END_REF], machine learning [START_REF] Drineas | On the nyström method for approximating a gram matrix for improved kernel-based learning[END_REF][START_REF] Williams | Using the nyström method to speed up kernel machines[END_REF], genetics [START_REF] Michael | Cur matrix decompositions for improved data analysis[END_REF]. For further references on their applications to data analysis please refer to Mahoney's survey article [START_REF] Michael | Randomized algorithms for matrices and data[END_REF]. Randomized algorithms usually divide in 2 classes: random sampling-and random projection-based algorithms. Both approaches rely on building a sketch version of the input matrix in order to extract relevant information. The random sampling (also called column selection) approach basically consists in selecting a subset of columns and rows, while the random projection more generally consists in multiplying the input matrix by a sketching matrix, that can be either sparse or dense. During this thesis we mainly focused on a dense random projection based LRA algorithm called the Randomized SVD (rSVD), that we present in details in the following subsection.

The Randomized SVD

Standard Decomposition First of all, we recall the formalism of the standard SVD. Let C be an arbitrary m-by-n real matrix, then its SVD reads as

C = UΣV T (3.1)
where U ∈ R m×m , resp. V ∈ R n×n , denotes the set of m left singular vectors, resp. n right singular vectors, associated with the singular values of C, that we denote σ(C). Finally, the diagonal matrix Σ ∈ R m×n contains the singular values of C, i.e.,

Σ ij = δ ij σ i (C), ∀i, j ∈ {1, . . . , n} .
We denote C r the matrix obtained after truncation of the full SVD representation at rank

r, i.e., C ≈ C r = U r Σ r V T r (3.2)
where

U ∈ R m×r , resp. V ∈ R n×r ,

Principle & Algorithm

The rSVD is a very popular random projection-based LRA algorithm introduced and further enhanced in a series of papers [START_REF] Per-Gunnar Martinsson | A randomized algorithm for the approximation of matrices[END_REF][START_REF] Liberty | Randomized algorithms for the low-rank approximation of matrices[END_REF][START_REF] Franco Woolfe | A fast randomized algorithm for the approximation of matrices[END_REF], that provides approximate low-rank representation in SVD form at a quadratic cost. The 2 stage of the rSVD for a symmetric input matrix C ∈ R n×n given a prescribed numerical rank r are summarized in Algorithm 6, they consist in:

• A randomized range approximation: Form Q ∈ R n×r , an approximate length-r basis for the range of C.

-Perform a Gaussian random projection, i.e., application of C to a n-by-r Gaussian random matrix

Ω Y = CΩ (3.4)
-and a subsequent orthogonalization, e.g., a thin QR Decomposition, such that

QR = Y

with an upper triangular matrix R ∈ R r×r and an orthogonal matrix Q ∈ R n×r .

-Hence, Q approximate the range of C and we can get a rank-r approximation of

C in the form C = QQ T CQQ T (3.5)
with tight error bounds that hold with high probability.

• A standard matrix factorization: In order to provide C in SVD form, i.e., C = Ũ Σ ŨT , we perform the following operations:

-We start by assembling the r-by-r matrix B such that

B = Q T CQ (3.6)
-then we compute the full SVD of the small matrix B, i.e.,

B = U B Σ B U T B . (3.7)
-Finally, we obtain an approximate rank-r SVD of C by forming the exact rank-r SVD of C, i.e., Ũ = QU B and Σ = Σ B .

A comprehensive review of the method can be found in Halko et al. [START_REF] Halko | Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions[END_REF]. In this review, we learn that the algorithm is fairly generic and thus easily extends to other factorizations such as the Cholesky Decomposition or the Interpolative Decomposition. Moreover, it provides various improvements such as iterative and pass-efficient variants. The full algorithm is presented in Algorithm 6, please refer to [START_REF] Halko | Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions[END_REF] for the details of the fixed rank and fixed accuracy approximate range finders.

Asymptotic complexity Since the overall cost of the rSVD is dominated by the cost of the r matrix multiplications required in each stage, the algorithm has a O(n 2 ×r) asymptotic complexity in time. While standard LRA algorithms have a cubic complexity, algorithms that produce truncated low-rank presentations (like the Arnoldi SVD or the pivoted QR Decomposition) usually have a quadratic complexity. However, the rSVD is supposed to have a much lower running time than Arnoldi's variant. We verify this assumption on numerical benchmarks in Section 7.4. Moreover, randomized algorithms leave significant room for optimizations and parallelization.

Algorithm 6:

The Randomized SVD Input: symmetric matrix C ∈ R n×n , rank r or accuracy ε, number of power iterations q, oversampling parameter s Output:

[ Ũ, Σ] approximate SVD of C // Stage I: Approximate the range of C if accuracy ε is prescribed then [Q, r] = ARRF (C, n, q, s, ε) O(r × n 2 ) if rank r is prescribed then [Q, ε] = RSI(C, n, q, s, r) (idem) // Stage II: Decompose C = Q(Q T CQ)Q T as Ũ Σ ŨT Build B = Q T CQ ∈ R r×r O(r × n 2 ) Perform SVD of B = U B Σ B U T B O(r 3 ) Form Ũ = QU B and Σ = Σ B O(n × r 2 )
Randomized range finders Algorithm 6 implements the random projection during the randomized range approximation. It can be performed either directly using one large matrixto-matrix multiplication or iteratively by means of multiple small matrix-to-matrix multiplications. More precisely, if the rank of the matrix is known a priori, then the Randomized Subspace Iterations (RSI) algorithm can be used to approximate the range directly given a prescribed rank r, see Algorithm 4.4 in [START_REF] Halko | Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions[END_REF]. However, if the rank is not known in advance and accuracy has to be monitored, the Adaptive Randomized Range Finder (ARRF) can be used to return a near-optimal range approximation given a prescribed accuracy ε, see Algorithm 4.1 in [START_REF] Halko | Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions[END_REF]. Due to their random nature these algorithms have a non-zero chance to fail, however error bounds such as (3.11) generally hold with high probability, e.g., at least 1 -n × 10 -r for the RSI algorithm. For instance, a rank-r approximation of a symmetric matrix of size n = 10 6 has a 1 chance to fail out of 10 10-6 = 10 4 .

Improving accuracy

In order to improve the accuracy of the range approximation the projection can be performed on a slightly larger subspace, namely using a n-by-(r + s) random matrix Ω where s is called the oversampling parameter, and finally keeping the first r columns of Q. Furthermore, the original algorithm does not apply well to matrices with a slow decreasing (or relatively flat) spectrum, as it fails to identify the most important singular values. A common alternative is to stretch this spectrum by using q subspace iterations, i.e., by considering (CC * ) q C instead of C as the input matrix. The resulting C

approximates C within well-established and controlled error bounds, that can be expressed in Frobenius or Spectral norm. First of all, the error of the rSVD, namely ∥C -C∥, cannot be lower than the error of the best rank-r approximation of C provided by the SVD, namely

∥C -C r ∥, which translates to ∥C -C r ∥ ≤ ∥C -C∥ (3.8)
This deterministic lower bound is also called the baseline in [START_REF] Halko | Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions[END_REF] and it reads as

∥C -C r ∥ F = ( n ∑ i=r+1 σ 2 i (C) ) 1/2 (3.9)
in Frobenius norm and

∥C -C r ∥ S = |σ r+1 (C)| (3.10)
in Spectral norm. An upper bound for the average error in Frobenius norm can be expressed in term of the baseline as

E(∥C -C∥ F ) ≤ f F (r, s, q)∥C -C r ∥ F , ( 3.11)
where f F is a polynomial function of the rank r that ensures tighter error bounds in average as s and q grow, see [START_REF] Halko | Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions[END_REF] for a comprehensive discussion on error bounds and explicit expressions of f F and f S , the equivalent in Spectral norm.

Fast Randomized LRA

In this section we give a brief overview of the existing fast algorithms for performing LRA of matrices using randomized techniques.

Fast random projection

The acceleration of randomized LRA algorithms based on random projection usually relies on the ability to apply fast matrix multiplication, hence most research works in this area focus on exploiting the structure and the sparsity of the input matrix C and the random matrix Ω.

Random matrix Random projection algorithms exploiting the nature of the random matrix are presented by Woodruff [START_REF] David | Sketching as a tool for numerical linear algebra[END_REF] in a comprehensive survey within the general framework of matrix sketching. While Sarlós [START_REF] Sarlos | Improved approximation algorithms for large matrices via random projections[END_REF] first mentioned the idea of using structured random (sketching) matrices such as the Fast Johnson Lindenstrauss Transforms (FJLT), Ailon et al. [5] described this algorithm in full details. The idea is based upon Achlioptas's work [START_REF] Achlioptas | Database-friendly random projections: Johnson-lindenstrauss with binary coins[END_REF] that aimed at developing sparse sketching for more efficient random projection. Li et al. [START_REF] Li | Very sparse random projections[END_REF] in turn developed a very sparse random projection-algorithm based on this work. Woolfe et al. [START_REF] Franco Woolfe | A fast randomized algorithm for the approximation of matrices[END_REF] subsequently designed an FJLT algorithm based on the Discrete Fourier Transform, namely the Subsampled Randomized Fourier Transform (SRFT). The ability to compute such transform in quasi-linear time allowed for reducing the computational cost of the random projection from O(n 2 × r) to O(n 2 × log r). Such method can be used to compute LRAs of matrices, in fact in Liberty et al. [START_REF] Liberty | Randomized algorithms for the low-rank approximation of matrices[END_REF] the SRFT is used to speedup the randomized Interpolative Decomposition introduced in Martinsson et al. [START_REF] Per-Gunnar Martinsson | A randomized algorithm for the approximation of matrices[END_REF]. However, it is not clear yet if there exists a fixed accuracy variant for fast transform based projection schemes. Various efficient algorithms for performing dense random projection are studied in Liberty's thesis [START_REF] Liberty | Accelerated dense random projections[END_REF][START_REF] Liberty | Dense fast random projections and lean walsh transforms[END_REF] such as the general FJLT, the SRFT, the Subsampled Randomized Hadamard Transform (SRHT) (based on the Discrete Walsh-Hadamard Transform [START_REF] Boutsidis | Improved matrix algorithms via the subsampled randomized hadamard transform[END_REF]) as well as the Mailman algorithm for matrix multiplication. Finally, Tropp [114] provided an improved error analysis for the SRHT.

Input matrix To our knowledge, very few articles address the acceleration of the matrix multiplications in Gaussian random projection-based algorithm for dense structured input matrices. In the context of low-rank approximations of matrices, Martinsson [START_REF] Martinsson | A fast randomized algorithm for computing a hierarchically semiseparable representation of a matrix[END_REF] describes a way to build HSS matrices using a Fast Multipole Method (FMM)-powered Gaussian random projection algorithm. The FMM is used there in order to accelerate the projection stage, i.e., the product between the input matrix and the Gaussian random matrix (3.4). Then, the factorization in HSS form is handled by the Interpolative Decomposition.

Fast random sampling

Block structure Random sampling techniques are particularly well-suited for applications requiring relatively low accuracy LRA. The most widely used random sampling algorithm for spd matrices is the Nyström method, see [START_REF] Drineas | On the nyström method for approximating a gram matrix for improved kernel-based learning[END_REF][START_REF] Gittens | Revisiting the nystrom method for improved large-scale machine learning[END_REF][START_REF] Zhang | Improved nyström low-rank approximation and error analysis[END_REF]. Applications of the Nyström method to the LRA of covariance kernel matrices can be found in Wang et al. [START_REF] Wang | Efficient sdp inference for fully-connected crfs based on low-rank decomposition[END_REF]. An improved variant of the original Nyström method introduced in Si et al. [START_REF] Si | Memory efficient kernel approximation[END_REF] exploiting the block structure of kernel matrices is known as the Memory Efficient Kernel Approximation (MEKA). It relies on a pre-clustering of the data and a block-wise LRA. Very recently, Wang et al. [START_REF] Wang | Structured block basis factorization for scalable kernel matrix evaluation[END_REF] provided an algorithm based on a similar idea and called the structured Block Basis Factorization (BBF), that does not have stability issues and exhibits smaller standard deviation than the MEKA.

Hybrid approaches

The complementarity of random sampling and random projection algorithms have often been highlighted in the literature [START_REF] Gittens | Revisiting the nystrom method for improved large-scale machine learning[END_REF]. Therefore, the design of an efficient randomized algorithm often relies on the relevant combination of both approaches.

The most obvious example of fast methods for random sampling algorithm is probably the fast approximation of matrix coherence using the randomized SVD proposed by Drineas and Mahoney [START_REF] Drineas | Fast approximation of matrix coherence and statistical leverage[END_REF]. On the other hand, fast random projection techniques such as the SRHT or the SRFT are tightly connected to improved random column selection. Finally, [START_REF] Wang | Structured block basis factorization for scalable kernel matrix evaluation[END_REF] makes use of random sampling, random projection and clustering in order to design an efficient O(n) factorization of a covariance kernel matrix with a low dependence on the ambient dimension. Since it is taylored for the high dimensions involved in machine learning applications, the BBF should perform worse than the FMM in 3D applications.

4

Positioning & Contributions

Efficient interpolation based FMM

Many applications of scientific computing require fast methods for applying large scale dense matricial operators, that can efficiently deal with the dimensionality of the problem. Such methods should be fairly generic and not too intrusive, in order to allow more flexibility and provide a better applicability to a wider range of problems and formulations. In particular, hierarchical methods such as kernel-independent Fast Multipole Method (FMM) allow linear time computations of kernel matrices involved in computational physics, geostatistics,… On the other hand, randomized low-rank approximation techniques provide a powerful generic framework to compute approximate standard decomposition of large low-rank matrices. Here, we describe the context of our research and our motivations for developing faster FMM schemes and randomized Low-Rank Approximation (LRA) techniques.

Efficient interpolation based FMM

As mentioned in Section 1.2, the bbfmm suffers from a suboptimal low-rank representation of the M2L operators [START_REF] Fong | The black-box fast multipole method[END_REF], that is mainly due to the 3D interpolation procedure. Even if M2L operators can be re-compressed [START_REF] Messner | Optimized m2l kernels for the chebyshev interpolation based fast multipole method[END_REF], the algorithms are usually relatively involved and computationally intensive. For instance, they can significantly increase the precomputation time for an arbitrary interaction kernel, as shown by our numerical benchmarks in Section 5.4. Moreover, these optimizations are often kernel-dependent, in particular they depend on certain features of the kernel such as its symmetry or its homogeneity. Finally, these optimizations are based on low-rank approximations, therefore they need to be carefully tuned in order to control the extra error they introduce in the representation of the kernel matrix. Consequently, there is a great need for a generic and exact method that allows for a faster precomputation and application of the M2L with lower memory footprint, i.e., a more optimal low-rank representation of the M2L operators. There is at least 2 options, that we could consider in order to find an interpolation scheme that produces a more optimal rank. The first one would be to use less interpolation points but selecting them more cleverly in order to maintain accuracy, as described in Casenave [START_REF] Casenave | An empirical interpolation based fast summation method for translation invariant kernels[END_REF]. Another option would be to use interpolation grids, that produce structured M2L operators. If those operators can be easily express in low-rank form, then they can be applied faster with a lower memory footprint. In Chapter 5 we present a new efficient interpolation based FMM, called the ufmm, that dramatically reduces the cost of the M2L and the overall running time of the FMM. This method relies on an equispaced grid therefore it becomes instable for large interpolation orders, however many real life applications only require relatively low accuracy and can thus fully benefit from the performance of the ufmm. In the following chapters we address 2 applications requiring fast matrix multiplications for either tensorial kernels or scalar kernels with many right-hand-sides.

Fast Multipole DD simulations

Although most existing Dislocation Dynamics (DD) codes implement a Fast Mutipole evaluation of the elastic fields, the formulations often lack efficiency and genericity. First of all, FMM formulations based on Spherical Harmonics or Taylor expansions usually have slow convergence, therefore they require a large number of modes to provide a relatively low accuracy on the force field evaluation. Second of all, these methods depend heavily on the kernel, therefore switching from one representation of the elastic field to another requires a new implementation and a significant amount of work. In order to provide efficient farfield evaluation of the isotropic elastic fields, we implemented interpolation based FMM in OptiDis, a DD parallel code using ScalFMM as an FMM engine. However, the tensorial nature of the interactions makes the problem inherently expensive to compute. Therefore, the bbfmm may not perform well and we may require more efficient methods such as the ufmm. Additional work should also be considered in order to reduce the effect of computing tensorial farfield interactions. Finally, the main purpose of this work is not to compete with existing DD codes but rather to propose a different alternative and extend the capabilities of OptiDis. Moreover, performing comparisons between the state-of-the-art codes is in our opinion neither easy nor relevant. In fact, all implementations exhibit major differences in term of parametrization and tuning. Additionally, most existing codes are not open-source, including OptiDis, and the few available releases, e.g., ParaDis, are far from having the optimal performance of the development branch.

Fast Randomized LRA

The growing interest for randomized LRA methods and the lack of efficient software solution brought us to develop an optimized code implementing fast numerical methods for random projection based algorithms. In particular, we provided fast algorithms for 2 applications involving the computation of large covariance matrices given as correlation kernels or computed from distance matrices. Let us describe the reasons that motivated us to address such problems in full details.

Gaussian Random Fields First of all, many scientific applications based on Monte Carlo approaches arising from geosciences, cosmology or image processing, rely heavily on intensive covariance matrix computations. A very usual issue encountered in such applications is the generation of large ensembles of multivariate Gaussian Random Variables, a.k.a., Gaussian Random Field (GRF), given an input correlation kernel. In particular, the efficient computation of a square root of the associated covariance matrix cannot be addressed by standard matrix decomposition algorithms. The existing alternatives usually have strong numerical limitations and often do not benefit from the shape of the distribution of points. Randomized LRA techniques provide convenient alternatives, since they benefit from the genericity of the standard matrix decompositions but they simply rely on matrix multipli-cations. Furthermore, a hierarchical method such as the FMM allows for decreasing the asymptotic complexity of the method from quadratic to linear in time, while avoiding the assembly of the full matrix and benefiting from the spatial distribution of points. Finally, the algorithm can benefit from the genericity and performance of the ufmm in order to be applied to any kernel at a low computational cost. Hierarchical methods are becoming more and more popular for performing faster covariance matrix computations [6,7] but are only rarely combined with random projection techniques, we propose to overcome this shortcoming.

Taxonomy Our last contribution is an attempt at extending the capabilities of our library to arbitrary matrices arising from biological applications. The classification of biological species can be addressed using dimensionality reduction algorithms such as the Multidimensional Scaling (MDS). More precisely, MDS allows for visualizing a point cloud in a low dimensional subspace given a distance matrix, by simply computing a square root of the associated covariance matrix, better known as similarity matrix. Since distance matrices are computed from real-life data, namely samples of O(10 5 ) reads coming from new generation DNA sequencing, they result in large but low-resolution covariance matrices. However, most fast MDS implementations rely on column selection based LRA, i.e., on the Nyström method, therefore they are often restricted to spd matrices and may require a significant amount of columns in order to represent the point cloud with sufficient accuracy. On the other hand, a random projection based LRA such as the randomized Singular Value Decomposition (SVD) is a well-suited and competitive algorithm for the computation of such square root, since it also gives useful extra information on the data and the structure of the similarity matrix. In Chapter 7, we further discuss the benefits of using random projection over common algorithms for such applications and its potential complementarity with the state-of-the-art approaches.

Part II Efficient kernel methods for n-body problems

In this chapter, we introduce a new computationally efficient Fast Multipole Method (FMM) that reduces the cost of the M2L operators, namely the ufmm (Algo. 7). We also derive a block low-rank algorithm, called the smooth-ufmm, that is optimized for globally smooth kernels. Then, we evaluate the theoretical complexity and memory requirements of these algorithms and compare them to the bbfmm. Finally, we evaluate the accuracy and running times on artificial numerical benchmarks and we analyze the relative performance of the ufmm and smooth-ufmm compared to optimized variants of the bbfmm.

A new FFT-accelerated FMM, the Uniform FMM

In the present section, we show that interpolating on equispaced grids allows for dramatically decreasing the memory footprint and the computational cost of the M2L operators in the bbfmm, namely from O(p 6 ) to O(p 3 log p) for precomputation and O(p 3 ) for application. First of all, we describe the basics of the method and its usual limitations. Then, we discuss the structure of the associated M2L operators and improvements based on Fast Fourier Transform (FFT). Finally, we present the algorithm as well as optimizations and implementation details.

A new FMM based on equispaced interpolation grids

Uniform interpolation grid For the sake of clarity, let us first consider the 1D case. The interpolation formula based on an equispaced grid still reads as (1.9) but involves different expressions for the interpolation points x and the interpolators S(x, x). More precisely, the interpolation points are defined as

xm = -1 + 2m p , ∀m = 0, . . . , p (5.1)
and the polynomial interpolators, a.k.a., the Lagrange polynomials, take the following form

S n (x) ≡ S(p; x, xn ) = p ∏ m=0 m̸ =n x -xm xn -xm , ∀x ∈ [-1, 1] , ∀n = 0, . . . , p. (5.2)
Addressing Runge phenomenon Interpolation schemes based on equispaced grids usually become unstable for high orders of interpolation p (Runge phenomenon), however as described in a review by Boyd et al. [START_REF] John | Exponentially-convergent strategies for defeating the runge phenomenon for the approximation of non-periodic functions, part i: Singleinterval schemes[END_REF] efficient regularization methods already exist. They are for instance based on optimization techniques such as Tikhonov regularization [START_REF] John | Defeating the runge phenomenon for equispaced polynomial interpolation via tikhonov regularization[END_REF], series expansions such as Gegenbauer regularization [START_REF] John | Trouble with gegenbauer reconstruction for defeating gibbs' phenomenon: Runge phenomenon in the diagonal limit of gegenbauer polynomial approximations[END_REF], subsampling techniques such overdetermined least-squares and Mock-Chebyshev interpolation [START_REF] John | Divergence (runge phenomenon) for least-squares polynomial approximation on an equispaced grid and mock-chebyshev subset interpolation[END_REF]) or even multi-domain approaches [START_REF] John | Exponentially-convergent strategies for defeating the runge phenomenon for the approximation of non-periodic functions, part two: Multi-interval polynomial schemes and multidomain chebyshev interpolation[END_REF]. Moreover, these methods defeat Runge phenomenon while preserving subgeometric convergence. In our approach divergence is not likely to occur since we use interpolation on multiple subdomains and in almost all cases of interest the value of p re-mains relatively small, e.g., p < 20. In particular, machine accuracy was reached with our algorithm for a wide range of non-oscillatory kernels, e.g., 1/r, 1/r 2 , correlation kernels (Section 5.4) and various isotropic elastostatic Green's functions arising in Dislocations Dynamics (Section 6.4.2).

Optimization of the M2L using the Fast Fourier Transform

Structure of the M2L operators The M2L operators associated to a uniform interpolation grid exhibit a well-known structure that can be used to provide a more efficient fast multipole summation scheme. Indeed, since k(•, •) is evaluated on equispaced 1D source (respectively target) grids denoted {x i } i=0...p (respectively {ȳ i } i=0...p ) and only depends on the distance between grid points, i.e.,

k(x i , y j ) = k(x i -y j ) (5.3)
for all i, j = 0, . . . , p, then each diagonal of K contains constant values (see Fig. 5.1 left). We say that the M2L operator K is a Toeplitz matrix. In the case of 2D grids the resulting matrix is block Toeplitz, i.e., the matrix is composed of constant blocks over its diagonals while each block is itself a Toeplitz matrix (see Fig. 5.2 left). In 3D, we introduce an extra level of blocking meaning that the constant diagonal blocks are now block Toeplitz.

Circulant embedding

The embedding of a Toeplitz matrix in a larger circulant matrix is a very common technique that allows for efficiently applying a matrix to vector product as a convolution in Fourier space. Here, we briefly recall this method in the case of a 1dimensional M2L operator K of order p. First of all, since K is Toeplitz, it is fully defined by its first row

R = (R 0 , . . . , R p ) ∈ R p+1 , with R i = k(x 0 , xi ) for i = 0, . . . , p.
and its first column

C = (C 0 , . . . , C p ) t ∈ R p+1 , with C i = k(x i , x0
) for i = 0, . . . , p.

These row and column can be embedded into a row R ∈ R p defined as

R = (R 0 = C 0 , . . . , R p , C p , . . . , C 1 )
with p = (p + 1) + p = 2p + 1. This new row can be used to generate a circulant matrix Ē ∈ R p×p , that embeds K in its upper left corner (see Fig. 5.1 blue frame on the right).

As a result, the application of both matrices to a vector is equivalent, if the last p columns and rows of Ē are masked. On the other hand, if the Toeplitz matrix is symmetric, then it is uniquely defined by its first row R and the associated the embedding is slightly smaller, namely p = p + 1 + (p -1) = 2p. The circulant embedding of a 1D (resp. 2D) symmetric Toeplitz matrix is illustrated on the right-hand-side of Figure 5.1 (resp. Figure 5.2). If M ∈ R p denotes the vector obtained after padding M with p zeros, then for i = 0, . . . , p we have

Conversion to Fourier domain

(L) i = ( Ē M) i = ( F * diag(Λ)F M) i = ( F * [ Λ ⊙ F M]) i = ( F * [ F R ⊙ F M]) i
Hence, the matrix to vector product (5.4) can be performed in the discrete Fourier space in the form of an entrywise product F R ⊙ F M. The same method applies to the 3D case, except that the embedding has to be performed independently in each dimension. In particular, the expansions and the row R need to stored as 3D arrays and transformed by 3D FFT. Transfers back and forth between Fourier and physical domains are done by forward and backward FFT resulting in an asymptotic overall cost of O(p 3 log p) for precomputation and O(p 3 ) for application.

a 0 a 1 a 2 a 4 a 3 a 3 a 0 a 1 a 2 a 4 a 4 a 3 a 0 a 1 a 2 a 2 a 4 a 3 a 0 a 1 a 1 a 2 a 4 a 3 a 0                                   L 0 L 1 L 2 X X                                   M 0 M 1 M 2 0 0                                   apply dense M2L F                                   Λ                                   a 0 a 1 a 2 a 4 a 3                                   buildM2L Figure 5
.3: Le : Mul plica on of a 1D M2L operator for p = 2 (in blue) by a vector, a er embedding it in a circulant matrix. Right: Precomputa on of M2L operator in the Fourier domain.

F                                   M                                   M 0 M 1 M 2 0 0                                   transformM Λ 0 Λ 1 (0) Λ 2 (0) Λ 3 Λ 4                                   L                                   M                                   applyM2L F *                                   L 0 L 1 L 2 X X                                   L                                   transformL

Algorithm and efficient implementation of the ufmm

Let us show how we can efficiently adapt such feature in a Fast Multipole implementation by describing various optimizations and then discuss the associated gains in terms of computational cost and memory footprint. The resulting algorithm, called the ufmm, is presented in Algo 7.

Algorithm 7: ufmm functions

Function precomputeM2L() // Assemble M2L operators (in Fourier domain)

for source cell C (L) y ∈ I(C (L)
x ) do // Compute first row and column of 3D block Toeplitz M2L operators

R β = k(x 0 , ȳβ ), C β = k(x β , ȳ0 ), ∀β/|β| ≤ p // Embed R in the first row R of a 3D block circulant matrix for β/|β| < p = 2p + 1 do if β d > p then Rβ = R β ′ with β ′ d = 2p + 1 -β d and β ′ i = β i for i ̸ = d else Rβ = C β // Apply 3D DFT R = F R with F αβ = e -2iπ p (α•β) , ∀(α, β)/|α|, |β| < p Function applyM2L( L(C (L) x ), C (L)
x ) // Apply M2L operators (in Fourier domain)

for source cell C (L) y ∈ I(C (L) x ) do Lα (C (L) x )+ = δ αβ Rα Mβ (C (L) y ), ∀α/|α| < p Function transformM(M)
// Pad expansion with zeros and transfer to Fourier domain.

M = 0 p Mβ = M β , ∀β/|β| ≤ p M = F M Function transformL( L)
// Transfer back to physical domain and unpad.

L = F -1 L L α = Lα , ∀α/|α| ≤ p
Fast Fourier Transforms In order to design an efficient algorithm we need to perfom as few operations as possible during the M2L step, therefore the forward transformations of the multipole expansions (resp. backward transformation of local expansion) are done during the P2M/M2M (resp. L2P/L2L) steps, as in the compressed-bbfmm (Algo. 4). Furthermore, we only precompute the first row of the circulant embedding of M2L operators and immediately transform and store them in Fourier domain. Consequently, the only operations that remain in the M2L steps are entrywise products between complex valued vectors of size p. However, in order to minimize the cost of storing larger complex valued expansions, it is recommended to apply the backward transformations right after the entrywise product with the M2L operators. Hence, we avoid the storage of the transformed local expansions, which can be crucial for large tree depths (see Section 5.3.3).

Real valued transforms

In the ufmm, the cost of storing and applying the M2L operators decreases from O(p 6 ) to O(p 3 ) with a slightly larger constant that can be divided by 2 using symmetries of the DFT. In 1D, since both the expansions and the M2L operators are real valued arrays, their 1D DFTs return conjugate-even * vectors. Consequently, only the first entry and the first half of the complex valued entries of the conjugate-even vectors need to be stored, i.e., p + 1 complex entries for an input real valued array of size 2p + 1. In multiple dimensions, this data packing only affects the last dimensions, e.g., in 2D the output of the 2D DFT of a (2p + 1)-by-(2p + 1) real valued array is a (2p + 1) × (p + 1) complex valued array. Hence, the storage of the M2L and the multipole expansions is roughly divided by 2.

Moreover, since the application of the M2L is a simple entrywise product, the output as a similar format and the computational time is also divided by 2.

1D transform

The block circulant embedding is a convenient approach, since it is relatively easy to visualize (at least in 1D or 2D). However, there exists slightly different embedding, that may be harder to visualize but lead to a fully circulant embedding matrix.

Let us illustrate such embedding on a simple 2D example with p = 2, where the M2L operator is defined as

K =       a b d e c a f d g h a b i g c a       . (5.5)
Then, the block and fully circulant embedding matrices are respectively defined as

Ē =                    a b c

d e f g h i c a b f d e i g i b c a e f d h i g g h i a b c d e f i g i c a b f d e h i g b c a e f d d e f g h i a b c f d e i g i c a b e f d h i g b c a

                   and Ē′ =                    a b f

d e i g h c c a b f d e i g h h c a b f d e i g g h c a b f d e i i g h c a b f d e e i g h c a b f d d e i g h c a b f f d e i g h c a b b f d e i g h c a

                   .
(5.6) Therefore, it is possible to use the 1D FFT to compute the DFT of the M2L operators. As a result, we can benefit from better performance of the FFT and also further reduce the size of the transformed arrays. In fact, instead of packing in the last dimension, we 5.2. A new block low-rank algorithm for smooth kernels, the smooth-ufmm pack the entire complex arrays. More precisely, after transformation the M2L operators are stored as complex valued rows of size ((2p + 1) 3 -1)/2 + 1 = 4p 3 + 6p 2 + 3p + 1 instead of (2p + 1) × (2p + 1) × (p + 1) = 4p 3 + 8p 2 + 5p + 1, namely 2p(p + 1) extra entries. On the other hand, the bbfmm stores M2L operators as real valued matrices of size (p + 1) 3 × (p + 1) 3 , i.e., (p + 1) 6 real entries.

A new block low-rank algorithm for smooth kernels, the smooth-ufmm

The FMM introduced by Greengard et al. [START_REF] Greengard | A fast algorithm for particle simulations[END_REF] was originally developed for kernels with a strong singularity at the origin, e.g., k(x, y) = 1/|x -y|. Moreover, it relied on spherical harmonics expansions, thus requiring x and y to be sufficiently far appart. Therefore, the algorithm was designed to always ensure well-separation (γ = 2). For globally smooth kernels, e.g., the Gaussian kernel

k(r) = k ∞ (r) = e -r 2 /(2ℓ 2 ) , ∀ℓ ∈ R (5.7)
involved in the application of Chapter 7, we show that this condition can be relaxed in order to define a new efficient block low-rank algorithm.

Extension of the interaction list

Provided the kernel is sufficiently smooth near the origin we want to consider a summation scheme where the nearfield interactions are approximated as well. An algorithm of this nature can be derived from the ufmm algorithm by simply changing the admissibility criterion γ ( L) at the leaf level from 2 to 0. While in the ufmm well-separation of leaves (γ ( L) = 2) is imposed, in the new variant coïnciding leaves form admissible pairs as well (γ ( L) = 0). In fact, since the kernel is sufficiently smooth at the origin, the diagonal subblocks of the kernel matrix associated with coïncident and adjacent cell pairs are relatively low-rank. This variant is denoted smooth-ufmm and can be seen as a block low-rank approximation technique.

Figure 5.5 illustrates the difference between both variants in term of interaction list. In the smooth variant the interaction list at the leaf level I(C L x ) includes all leaf cells whose parent is a direct neighbor of the target parent cell C L-1

x and consequently N (C L x ) = ∅, whereas in the standard variant the direct neighbors of the target leaf cell form the nearfield. The interaction list of a given leaf now includes the direct neighbors and the leaf itself, which means that all interactions will be transferred by means of M2L operations and thus no interaction needs to be computed at the P2P step.

( L) x ufmm k 1/2 (|x -.|) ( L) x smooth-ufmm k ∞ (|x -.|)

Optimal setup

Since the P2P and the M2L do not compete anymore in the smooth-ufmm, these steps do not have to be balanced. Consequently, the concept of level is not crucial anymore but it remains convenient, and the actual tuning parameter becomes the width of the leaf cells. However, we still need to minimize the cost of the M2L by using the fewest number of clusters, i.e., the largest leaf cells. Therefore, our algorithm remains a multi-level scheme, where the depth for the octree is used to control the width of the leaves and thus the overall cost of the algorithm. The optimal setup for the algorithm is the lowest tree depth, that leads to a similar accuracy as the original ufmm scheme.

Comparative cost analysis

ufmm versus bbfmm

The theoretical complexities of the bbfmm and ufmm algorithms are given in terms of computational cost in Table 5.1 and memory footprint in Table 5.2. In the ufmm the precomputation of the M2L operators requires O(p 3 log p) operations (i.e., the cost of the FFT) while their application requires only O(p 3 ) operations (i.e., the cost of an entrywise product). All algorithms scale linearly in n but they optimize the M2L step differently. In fact, the cost of storing and applying M2L operators scales like O(p 3 ) in the ufmm and O(p 6 ) in the bbfmm, therefore we expect significant differences in memory requirements and computational times.

smooth-ufmm versus ufmm

Due to the extension of the interaction list in the smooth-ufmm, the maximum number of M2L operators to store and apply at the leaf level increases from 189 to 189 + 27 = 216 in 3D. However, as shown in the numerical benchmarks presented in Section 5.4, for a given accuracy the optimal number of level is usually slightly lower than the level required for the standard ufmm. Moreover, since low-rank representations replace the direct computations,

Comparative cost analysis

for sufficiently large n the overall computational time of the smooth-ufmm is expected to be lower than that of the ufmm.

bbfmm ufmm smooth-ufmm

P2P n 0 × n n 0 × n 0 P2M/L2P p 3 × n p 3 × n + p 3 log (p) × n/n 0 p 3 × n + p 3 log (p) × n/
× (4p 3 + 6p 2 + 3p + 1) P2M/L2P M/L n/n 0 × (p + 1) 3 n/n 0 × (p + 1) 3 n/n 0 × (p + 1) 3 P2M/L2P M/ L 0 n/n 0 × (4p 3 + 6p 2 + 3p + 1) × 2 n/n 0 × (4p 3 + 6p 2 + 3p + 1) × 2 P2P * K n 0 × n n 0 × n 0 Table 5.2:
Memory footprint at various steps of the ufmm, smooth-ufmm and bbfmm expressed in terms of the number of entries.

We obtain the actual number of Bytes by mul plying the number of entries by the size of the value type, i.e., 8Bytes for real valued entries in double precision arithme cs. The constant n 0 = n/8 L denotes the average number of par cles per leaf in a uniform distribu on of par cles.

Theoretical memory requirements

Although the M2L step is usually the most computationally intensive step of the FMM, because it involves a total of 2 d L × (6 d -3 d ) multipole to local transfers of expansions, where d is the ambiant dimension, only L × (7 d -3 d ) M2L operators have to be stored. On the other hand, the storage of the 2 d L expansions is usually the limitant factor in term of memory consumption. Figure 5.6 representS the memory footprint of the ufmm and the bbfmm w.r.t. the interpolation order p for growing values of L = 3, 5, 7. First of all, for large tree depths, e.g., L ≥ 7, the amount of memory required for the storage of the expansions in both algorithms dominates the overall cost. Most importantly, it even exceeds the ressource available on standard computers for intermediate interpolation orders, namely approximately 50GB for the ufmm and 100GB for the bbfmm with p = 7. In such regime, usually associated with distributions of several million particles, a distributed memory version of the algorithm is usually prefered. For lower tree depths, e.g., L ≤ 5, the M2L of the bbfmm dominates the overall cost due to its O(p 6 ) complexity. On the other hand, in the ufmm the storage of the expansions rapidely dominates the overall memory footprint, namely for L ≥ 5. These values should be qualified by the fact that for heterogeneous geometries the number of non-empty cells can be very low.

For instance, as shown on Figure B.2 for n = 10 6 par cles, the memory required for storing the expansions at level L = 7 could be divided by 100 for a unit sphere and 1000 for a prolate sphere.

Numerical benchmarks

In order to compare the efficiency of our new algorithms to existing variants of the bbfmm, we will consider artificial numerical benchmarks. The distributions of particles considered are presented in Section 5.4.1, they involve different filling densities in the octree. Our comparative analysis is organized as follows:

• We first discuss statistics on the filling of the octree for various geometries.

• Then, we compare the relative accuracy and numerical performance of the bbfmm and the ufmm using the Laplacian kernel.

• Finally, we compare the relative accuracy and numerical performance of the ufmm and the smooth-ufmm using a smooth kernel, namely the Gaussian kernel.

Distributions of points and kernels of interaction.

Distributions Let us consider 3 different geometries for the distribution of particles:

• The unitCube is a volumic distribution, in which the particles are distributed uniformly inside a 2 × 2 × 2 cube.

• The unitSphere is a surfacic distribution, in which the particles are distributed uniformly on a unit sphere.

• The prolateSphere is a surfacic distribution, in which the particles are distributed uniformly on a 0.2 × 0.2 × 2 ellipsoid, a.k.a. prolate sphere with ratio 10.

k(r) = 1/r, (5.8)
and then on a smooth and non-homogeneous kernel, namely the Gaussian kernel k(r) = e -r 2 2ℓ 2 , with ℓ ∈ R.

(5.9)

Tensorial variants of the Laplacian will be used in Chapter 6, while variants of the Gaussian will be used in Chapter 7. Before starting any performance analysis, the accuracy of the algorithm should always be validated by plotting the error ε f mm on the potential w.r.t. interpolation order p.

bbfmm vs. ufmm on the Laplacian kernel

In this subsection we discuss the relative accuracy and sequential performance of the ufmm and optimized variants of the bbfmm using the standard Laplacian kernel (5.8).

Accuracy

The convergence of the ufmm and smooth-ufmm is represented on Figure 5.9 (left) in term of the relative 2-norm of the error on the potential. As expected all variants of the bbfmm have the same accuracy and convergence rate. Moreover, the ufmm is slightly less accurate than the bbfmm due to the near minimax property of the Chebyshev interpolation. Nevertheless, we will show later on that it still performs better for a fixed accuracy.

k(r) = ℓ α r α k(r) = ℓ 2 r 2 , ℓ = 0.1 k(r) = ℓ r , ℓ = 0.1 k(r) = ℓ (r 2 +ℓ 2 ) 1/2 , ℓ = 0.1 k(r) =

Sequential performance

The relative computational costs of the bbfmm, with or without compression of the M2L operators (see [START_REF] Messner | Optimized m2l kernels for the chebyshev interpolation based fast multipole method[END_REF]), and the ufmm are represented on Fig. 5.10 using a fixed number of particles n = 20.000 distributed in the unit cube with L = 3 (i.e., n 0 ≈ 30). The ufmm outperforms the unoptimized bbfmm in terms of both computational time and memory requirements. Let us recall that the symmetric-bbfmm is taylored for symmetric kernels allowing for a massive reduction of the interaction list and that it involves individual compression of the M2L operators, while the compressed-bbfmm only involves a global compression of the M2L operators. As shown on the graphs, the ufmm and the symmetric-bbfmm have similar performance though the ufmm applies to any kernel. More precisely, the ufmm is faster in term of computational time but requires a little more memory than the symmetric-bbfmm.

Detailled timings Figure 5.11 shows the detailled running times of each step of the ufmm and the symmetric-bbfmm for fixed tree depth L = 3 w.r.t. the interpolation order. Consequently, we observe how the cost of the farfield as the accuracy grows, while the cost of the nearfield is fixed. The graphs also illustrate how reducing the cost of the M2L affects the overall application time in the ufmm. Finally, it shows that the other steps (P2M/L2P and M2M/L2L) are more expensive in the ufmm than in the bbfmm if the interpolation order is large, namely p > 7, but cheaper in the opposite case. Performance at fixed accuracy For kernels such as the Laplacian, we may need to increase the interpolation order of the ufmm by 1 in order to get the same accuracy as the bbfmm (Figure 5.9), in particular for accuracy below 10 -8 , i.e., p ≥ 7. In any cases the ufmm beats all bbfmm variants by a factor of at least 2 (for p = 7), that grows with the interpolation order.

ufmm vs. smooth-ufmm on Gaussian kernels

In this subsection we discuss the relative accuracy and sequential performance of the ufmm and the smooth-ufmm using Gaussian kernels (5.9) with length scales ℓ ∈ {0.5, 1.0}.

Accuracy

The convergence of the ufmm and smooth-ufmm is represented on Figure 5.9 (right) in term of the relative 2-norm of the error on the potential. In particular, we observe that the larger ℓ the faster the convergence and the lower the error, which directly results from the smoothness of the function when ℓ grows. For large values of p, this fast convergence allows us to witness another phenomenon. In fact, for any ℓ once the machine accuracy is reached the interpolation scheme may become instable. However, most application require precisions that are far lower than the machine accuracy.

Performance at fixed size The relative computational costs of the ufmm and the smoothufmm are represented on Figure 5.12 using a fixed number of particles n = 20.000 and L = 3, i.e., n 0 ≈ 30. First of all, the extra memory required for storing the extra M2L operators of the smooth-ufmm is negligible compared to the memory required by the storage of the M2L operators of the ufmm. Second of all, the application time grows faster for the ufmm than for the smooth-ufmm as the interpolation order grows. As a result, the application time is 10 times lower for the smooth-ufmm at p = 2, while it is equal to the application time of the ufmm for p = 9. Finally, in this configuration, the ufmm and smooth-ufmm exhibit rather similar performance, which may not be the case in practice. In fact, in practice the optimal tree depth may be lower for the smooth-ufmm, therefore its computational cost may be far lower. Performance at fixed accuracy In order to improve our comparison, we analyze the relative performance of the ufmm and the smooth-ufmm for a given accuracy of about 10 -5 . Figure 7.10 shows the application times for 3 different geometries for n up to 10 6 particles.

The smooth-ufmm unsurprisingly has better performance than the ufmm, since the cost of approximating the nearfield becomes cheaper as n increases. Furthermore, it beats the ufmm by a factor of 10 for sufficiently large n, while the ufmm only beats the bbfmm by a factor of 2. performance and accuracy of the algorithm were tested on artificial benchmarks illustrating the significant benefits of using the ufmm over the bbfmm and its optimized variants. Moreover, a new algorithm called the smooth-ufmm was designed for handling smooth kernels more efficiently. This variant differs substantially from the original FMM in the sense that it implements a different admissibility criterion but it does not depend on the interpolation procedure, therefore it could as well be derived from the bbfmm. In the next two chapters, the applications of the ufmm and the smooth-ufmm to highly expensive simulations will be considered, namely tensorial interactions (Chapter 6) and matrix-to-matrix multiplications (Chapter 7).

6

Fast Multipole DD simulations In this chapter, we present the interpolation formula of tensorial interaction kernels involved in elasticity problems such as DD simulations. Then, we introduce various Fast Multipole Method (FMM) summation schemes for the fast evaluation the isotropic elastic stress, forces and energy. We start by simple and naïve schemes and then we move on to more evolved M2L-optimized schemes. Finally, we evaluate the performance of these schemes on numerical benchmarks.

Fast Multipole computation of the isotropic stress and forces

Introduction

Consistency First of all, we would like to recall that a dislocation loop may cross the edges of the octree, thus part of its segments can contribute to the nearfield interactions while the other segments contribute to the farfield interactions. Since all existing expressions of the elastic stress and energy are only equivalent on a closed dislocation loop, we need to use the same formula for the evaluation of the farfield and for the nearfield, namely (2.4) for the stress and (2.8) for the energy. Therefore, the choice of well-suited expressions is crucial in an efficient implementation and will be discussed here in particular for the evaluation of the energy.

Balance For various computational reasons related to the velocity of the dislocations, the width of the leaf cells cannot be smaller than a few dislocation lengths, namely ωL = 3L max . Consequently, the DD codes implementing the FMM usually impose a maximum depth for the octree. For instance, if L max = 100 with a root bounding box of width ω 0 = 20000 then L ≤ 6. However, the number of nearfield interactions can grow dramatically if the number of segment increases while the depth remains constant. As a result, the balance between nearfield and farfield becomes impossible to maintain. Therefore, the efficient implementation of the nearfield requires a lot of care.

Interpolation of the elastic stress

Let us first define the tensorial interaction kernels k(x, y) as

k i (x, y) = r a,i , k ij (x, y) = r a,ij and k ijk (x, y) = r a,ijk , (6.1)
where r a = √ a 2 + r 2 , r = |r| = √ r i r i and r = xy. The notation r a,i denotes the partial derivatives of r a w.r.t. x i , i.e., r a,i = ∂r a ∂x i . Hence, the analytical expressions associated with the 3 variants of k(x, y) presented above read as

10 -3 10 -2 10 -1 10 0 0 0.1 0.2 0.3 L(ω min ) nearfield farfield a = 3 x = r 1 = r 2 k(r) k ij (r) for a = 3 and x = r 1 = r 2 k 11 (x) = 1 ra -x 2 r 3 a k 12 (x) = -x 2 r 3 a k 11 (x) + k 22 (x) = -
k i = r i r a , k ij = δ ij r a - r i r j r 3 a and k ijk = - δ jk r i + δ ik r j + δ ij r k r 3 a + 3 r i r j r k r 5 a . (6.
2)

The graphs on Figure 6.1 represent certain components of the 2 nd -order interaction tensor in 2D for a = 3 on a representative area, namely a few minimum leaf cell widths, w.r.t. the scaled distance. It shows the maximum variations, that can exhibit the kernel in the nearfield. In particular, we observe that the kernels decrease fast within the first few core widths a and become rather smooth in the farfield, which makes interpolation based FMM well-suited. Although the non-singular kernel is smooth by definition, it still exhibits brutal variations within the nearfield, just like the Gaussian kernel for ℓ = 0.1. Therefore, it cannot be approximated in the nearfield using the smooth-ufmm. Using the compact notations defined in (1.15), we denote k the interpolant of k of order p such that

k(x, y) = ∑ |α|≤p S α (x) ∑ |β|≤p Kαβ S β (y), ∀x, y ∈ R 3 . (6.3)
If we replace the expression (6.3) of the interpolant k in (2.4), we get an approximation of the elastic stress field σ, that reads as

σ ij (x) ≈ µ 8π (σ A ′ ij + 2 1 -ν σB ′ ij ), ∀(i, j) ∈ {1, 2, 3} 2 , ( 6.4)
where for all i, j ∈ {1, 2, 3}

σA ′ ij = σA ij + σA ji σB ′ ij = σB ij -δ ij σB ℓℓ and σA ij (x) = ∑ |α|≤p S α (x) ∑ |β|≤p ( Kαβ ) mℓℓ ε jmk (C ′ ) S β (x ′ )b ′ k t ′ i dx ′ (6.5) σB ij (x) = ∑ |α|≤p S α (x) ∑ |β|≤p ( Kαβ ) mij ε nmk (C ′ ) S β (x ′ )b ′ k t ′ n dx ′ (6.6)
Since (C ′ ) is a dislocation network discretized by the reunion of segments S s = [x 1 s x 2 s ], i.e., (C ′ ) = ∪ n s=1 S s , the contributions σA and σB can also be expressed as the sum of the contributions of each segment, namely

σA ij (x) = ∑ |α|≤p S α (x) ∑ |β|≤p ( Kαβ ) mℓℓ ε jmk ∑ [x 1 x 2 ]∈(C ′ ) b ′ k t ′ i ∫ x 2 x 1 S β (x ′ )dx ′ (6.7) σB ij (x) = ∑ |α|≤p S α (x) ∑ |β|≤p ( Kαβ ) mij ε nmk ∑ [x 1 x 2 ]∈(C ′ ) b ′ k t ′ n ∫ x 2 x 1
S β (x ′ )dx ′ (6.8)

A naïve Fast Multipole summation scheme

Equations (6.7) and (6.8) can be computed using a fast multipole summation scheme, that consists in the following steps:

• For σA , perform successively

(M A β ) mij = ε jmk ∑ [x 1 x 2 ]∈(C ′ ) b ′ k t ′ i ∫ x 2 x 1 S β (x ′ )dx ′ (P2M-SA) (L A α ) ij = ∑ |β|≤p ( Kαβ ) mℓℓ (M A β ) mij (M2L-SA) σA ij (x) = ∑ |α|≤p S α (x)(L A α ) ij (L2P-SA)
• For σB , perform successively

(M B β ) m = (M A β ) mnn (P2M-SB) (L B α ) ij = ∑ |β|≤p ( Kαβ ) mij (M B β ) m (M2L-SB) σB ij (x) = ∑ |α|≤p S α (x)(L B α ) ij (L2P-SB)
• Finally, compute σ using (6.4).

First of all, we notice that the multipole expansions (M B ) m for m ∈ {1, 2, 3} are the traces of { (M A ) mij } i,j≤3 , therefore they can be computed at almost no extra cost. Furthermore, the M2L operators used for σA are the traces of those used for σB , more precisely

( Kαβ ) mℓℓ = tr ( { ( Kαβ ) mij } ij )
Finally, both terms σX for X ∈ {A, B} share the same L2P step, i.e.,

σX = ∑ |α|≤p S α (x)L X α . (6.9)
Nodal forces The computation of the farfield contribution to the nodal forces presented in (2.6) can be computed using a different summation at the L2P step, namely

f 4 i = ∫ L 0 ∑ |α|≤p S α (x 3 + ξt)ε ijk (L α ) jℓ b ℓ t k ξ L dξ, ∀i ∈ {1, 2, 3} (L2P-F4)
Cost estimation In order to perform this summation, we need to store 30 M2L components, namely 3 for σA and 27 for σB . We need to store 27 components of M A β and 3 components of M B β . Then, we need to perform 54 M2L operations in order to compute all 9 components of σA and σB . Indeed, computing (M2L-SA) involves 3 × 9 = 27 M2L operations and computing (M2L-SB) involves 9 × 3 = 27 M2L operations. The summation schemes proposed here is relatively naïve and expensive, therefore it will only be used as a reference for more efficient schemes.

Improved farfield computation of the stress

Here we introduce optimized variants of the previous summation scheme based on symmetries in Section 6.2.1, refactoring of the interactions in Section 6.2.2 and shifted derivatives in Section 6.2.3. Finally, we discuss the relative costs of each variant in Section 6.2.4.

Considering symmetries

Symmetries First of all, in the hypothesis of isotropic elasticity the stress field is symmetric, only 6 components out of 9 need to be computed. Furthermore, since in the expression of the kernel of interactions k the partial derivatives can be applied in an arbitrary order, then ( Kαβ ) ijk is symmetric w.r.t. i,j and k. Consequently, σB is symmetric and only its upper triangular part needs to be assembled during the M2L step. On the other hand, it could be more convenient to assemble σA ′ at the P2M step instead of the L2P, in order to store only symmetric quantities.

Fast Multipole summation scheme Hence, for all (i, j) ∈ {1, 2, 3} 2 with i ≤ j the optimized summation scheme can summarized as • For σA ′ , perform successively

(M A ′ β ) mij = (M A β ) mij + (M A β ) mji (Symm-P2M-SA) (L A ′ α ) ij = ∑ |β|≤p ( Kαβ ) mℓℓ (M A ′ β ) mij (Symm-M2L-SA) σA ′ ij (x) = ∑ |α|≤p S α (x)(L A ′ α ) ij (Symm-L2P-SA)
• For σB ′ , perform successively

(M B β ) m = (M A β ) mnn (Symm-P2M-SB) (L B α ) ij = ∑ |β|≤p ( Kαβ ) mij (M B β ) m (Symm-M2L-SB) σB ′ ij (x) = (σ B ) ij -δ ij (σ B ) ℓℓ (Symm-L2P-SB)
• Finally, compute σ using (6.4).

Cost estimation

In order to perform this summation, we need to store 13 M2L components, namely 3 for σA and 10 for σB . We need to store 3 × 6 = 18 components for M A ′ and only 3 for M B . Then, we need to perform 36 M2L operations in order to compute 6 upper triangular components of σA ′ and σB ′ . Indeed, (Symm-M2L-SA) involves 3 × 6 = 18 M2L operations and (Symm-M2L-SB) involves 6 × 3 = 18 M2L operations.

Refactoring interactions

The number of multipole and local expansions required by the two first schemes is relatively high, namely 39 for the Direct and 27 for the Symmetric scheme. As mentioned in Section 5.3.3, the amount of memory required to store 1 multipole expansion can quickly become prohibitive in interpolation based FMM. Therefore, we propose a new scheme based on expanding and refactoring the tensorial interactions and called the Factorized scheme, that simultaneously minimizes the number of M2L operations and the number of multipole expansions. For the sake of clarity, details of this scheme are presented in Appendix D.1.1.

If the energy needs to be computed as well, this variant is more convenient than using symmetries of the stress, since it implements the same multipole expansions, namely

(M β ) ij = ∑ [x 1 x 2 ]∈(C ′ ) b ′ i t ′ j ∫ x 2 x 1 S β (x ′ )dx ′ .
(Facto-P2M)

Shifting derivatives

Principle As mentioned in Section 5.1, the derivative of a function can be approximated by deriving the interpolated function. Thus, k ijk can be approximated by deriving kij , the interpolant of k ij . The associated interpolation formula reads as

k ijk (x, y) ≈ ∇ x k kij (x, y) = ∑ |α|≤p ∇ k S α (x) ∑ |β|≤p ( Kαβ ) ij S β (y) (6.10)
Since S is a higher order polynomials than ∇ k S for k ∈ {1, 2, 3}, then equation ( 6.3) should provide a more accurate approximation of k ijk than (6.10). On the other hand, (6.10) reduces the number of components to approximate from 10 for all k ijk to 6 for all k ij , while increasing the dimension of the interpolator by 3. Finally, if the derivation is done with respect to y, then it applies to the right interpolator. Moreover, since k(x, y) = k(x -y),

∇ y k kij = -∇ x k kij
and thus

k ijk (x, y) ≈ -∇ y k kij (x, y) = - ∑ |α|≤p S α (x) ∑ |β|≤p ( Kαβ ) ij ∇ k S β (y)
In other terms, for such interpolation formula, shifting derivatives from right to left implies changing the sign in front of the formula.

Fast Multipole summation scheme Shifting the derivatives w.r.t. m and applying the permutation symbol to the right interpolator in (6.7) and (6.8) results in a much cheaper summation scheme in terms of both memory and computational time. Let us show this by writing the new interpolation formula as follows

σA ij (x) = - ∑ |α|≤p S α (x) ∑ |β|≤p ( Kαβ ) ℓℓ ε jmk ∑ [x 1 x 2 ]∈(C ′ ) b ′ k t ′ i ∫ x 2 x 1 ∇ m S β (x ′ )dx ′ (6.11) σB ij (x) = - ∑ |α|≤p S α (x) ∑ |β|≤p ( Kαβ ) ij ε nmk ∑ [x 1 x 2 ]∈(C ′ ) b ′ k t ′ n ∫ x 2 x 1 ∇ m S β (x ′ )dx ′ (6.12)
Let us define the new multipole expansions associated with σA as

(M A β ) ji = ε jmk ∑ [x 1 x 2 ]∈(C ′ ) b ′ k t ′ i ∫ x 2 x 1 ∇ m S β (x ′ )dx ′
Since summation over m can now be performed at the P2M step, the dimension of this multipole expansion is divided by 3. Hence, for all (i, j) ∈ {1, 2, 3} 2 with i ≤ j, the associated fast multipole summation scheme consists in the following steps:

• For σA ′ , perform successively

(M A ′ β ) ji = (M A β ) ij + (M A β ) ji (Shift-P2M-SA) (L A ′ α ) ij = ∑ |β|≤p ( Kαβ ) ℓℓ (M A ′ β ) ij (Shift-M2L-SA) σA ′ ij (x) = - ∑ |α|≤p S α (x)(L A ′ α ) ij (Shift-L2P-SA)
• For σB ′ , perform successively

M B β = (M A β ) nn = tr(M A β ) (Shift-P2M-SB) (L B α ) ij = ∑ |β|≤p ( Kαβ ) ij M B β (Shift-M2L-SB) σB ′ ij (x) = (σ B ) ij -δ ij (σ B ) ℓℓ (Shift-L2P-SB)
• Finally, compute σ using (6.4).

Cost estimation

In order to perform this summation, we need to store 7 M2L components, namely 1 for σA and 6 for σB . Moreover, we need to store 6 components of M A ′ β and 1 components of M B β . On the other hand, one need to apply 12 M2L operations in order to compute all 6 components of σA ′ and σB ′ . Indeed, (Symm-M2L-SA) involves 6 M2L operations and (Symm-M2L-SB) involves 6 M2L operations. This cost reduction should compensate the need for a larger interpolation order induced by the loss of precision.

Comparative analysis of the theoretical cost

Here, we summarize the cost of each variant in two tables, please refer to Table 6.1 for the memory requirements and Table 6.2 for the computational costs. The tables show the constant that needs to be applied to the memory requirements or the computational costs of a given interpolation based FMM (bbfmm or ufmm), for each step and each operator.

Storage Table 6.1 shows that the variant using shifted derivatives requires the less overall memory. Concerning the two other optimized variants, using the symmetries requires approximately half of the memory required when refactoring the interactions for the storage of the expansions. On the other hand, the memory footprint of the M2L does not vary significantly for these 2 approaches. For the shifted variant, one should mention that the computation of the energy does nor require extra memory for the storage of the M2L operators since they share the same kernel of interaction. However different multipole expansions have to be stored for the stress and the energy, respectively 7 and 6.

Application In order to compare the computational cost of each tensorial approach, one should consider the constant in front of the cost of a scalar M2L written in Table 6 Discussion For a given variant, the ufmm is always faster than the bbfmm in theory. This assumption is verified experimentally on the numerical benchmarks presented in Section 6.4. Moreover, the ufmm require less memory for low tree depths. However, as the depth of the tree grows, the expansions require more memory than the M2L and the ufmm becomes more demanding. This may occur for lower depths than in the scalar case, since the ratio between the number of components of the expansions and the M2L exceeds 1 for all approaches. Figure 6.2 shows the total theoretical memory requirement of each variant for either the bbfmm or the ufmm. We observe that for p = 4 (a typical value used in OptiDis, that provides good accuracy) the profile does not change a lot due to the small variations of the ratio (between 1 and 2). The profiles only differ for the low tree depths where the cost of the M2L dominates in the bbfmm. For large tree depth, the memory footprints of the ufmm is always larger since the storage of the expansions dominates the overall memory footprint and it is substantially larger than that of the bbfmm. Based on previous observations, the optimal choice between the symmetric and the refactored variant may vary w.r.t. to the chosen interpolation grid. For instance, the refactored variant seems better suited for a FMM with fast M2L application such as the ufmm. Meanwhile, the symmetric variant can benefit from low memory requirements at the P2M step such as in the bbfmm. Finally, for a fixed interpolation grid, the variant using shifted derivatives has the lowest overall memory footprint. and for a fixed interpola on order p = 4. We consider the maximum amount of memory required for storing the 316 tensorial M2L operators, the mul pole and the local expansions as well as the transformed versions. The results are given in simple precision arithme c. Observa ons: The memory required by the M2L is almost not no ceable for the ufmm, in fact only the cost of storing the expansions is observed. For the bbfmm, the memory required for the storage of the M2L dominates the cost and can thus be observed for the first values of L.

Implementation details 6.3.1 The extended bounding box concept

Extending leaves Inserting the elements of a spatial discretization in an octree structure can sometimes be tedious, in particular when an element crosses an edge of the octree (see Figure 6.3). For instance, in most Fast Multipole BEM implementations [START_REF] Chaillat | Fast Multipole Method for 3-D elastodynamic boundary integral equations[END_REF][START_REF] Messner | Fast Boundary Element Method in Acoustics[END_REF], a boundary element is considered to be in a cluster, if the center of the element lies inside of it. Moreover, since interpolation is involved in our algorithm, it is crucial that the entire element lies inside the interpolation domain. Therefore, as suggested by [START_REF] Messner | Fast Boundary Element Method in Acoustics[END_REF] we define a slightly larger interpolation domain, called the extended bounding box, that contains all segments lying inside a cluster. Since extending the bounding box reduces the minimum distance between points of admissible clusters, it obviously affects the accuracy of the method. Associating segments to a cluster with respect to their center leads to the minimal width extension, i.e., the maximum segment length (see 6.3). Note that the center of the cell remains unchanged after extension of the width. Extending parent cells Once the leaf cells are extended all parent cells need to be extended as well, since interpolation is performed from child to parent and anterpolation from parent to child. The minimal parent's width extension, such that the extended child is contained inside the extended parent, is equal to the child's width extension. Therefore, we use the same width extension for all cells of all levels. Note that since the center remains unchanged, the extensions of adjacent cells intersect.

Numerical integration over segments

Farfield As described extensively in appendix D.3, during the P2M step (resp. L2P step), the integrals over the segments [x 1 x 2 ] (resp. [x 3 x 4 ]) can be computed exactly by mean of a simple Gaussian quadrature rule. For the integral involved in (P2M-SA) or (P2M-SB) the resulting quadrature formula reads as

∫ x 2 x 1 S β (x ′ )dx ′ = Q ∑ q=1 ω ′ q S(x ′ q ),
where Q denotes the quadrature size satisfying 2Q + 1 = p and ω ′ q (resp. x ′ q ) denotes the q-th Gauss weights (resp. points) mapped to the source segment [x 1 x 2 ]. More precisely, if ω q ∈ [0, 1] and θ q ∈ [-1, 1] denote standard Gauss weights and points, then mapped weights reads as ω ′ q = ω q L 2 and mapped points read as

x ′ q = x 1 + ξ(θ q )t = x 1 + (θ q + 1) L 2 t.
Nearfield As mentioned in the introduction 6.1.1, because of the condition imposed on the tree depth, the efficiency of the FMM is lost after a certain problem size. In fact, once the tree depth reaches its maximum value, the cost of computing the direct interactions starts to grow quadratically n. Therefore, it is crucial to have an efficient implementation of the P2P step. First of all, OptiDis already optimizes the evaluation of the nearfield using AVX instructions. Second of all, we implemented a new feature, that allows the direct interactions to be evaluated using numerical quadratures if segments are sufficiently far appart. Indeed, as mentioned in [8], using a quadrature can becomes cheaper than using analytical expressions for very distant segments. For instance, if the target accuracy is 10 -4 , then the required quadrature size equals 2 for d > 5L max and drops to 1 for d > 500L max . Since the maximum distance at the maximum level is d max = 6L max , this represents a small gain at the maximum level. However, this can significantly accelerate the nearfield computation for the low FMM levels if the root bounding box is very large compared to the maximum segment length. For instance, if L max = 100 with a root bounding box of width ω 0 = 100.000 then d max /L max = ω 0 /(2L max ) = 500 at level L = 2.

Enforcing homogeneity

Unlike the singular interaction kernel, the non-singular kernel is not homogeneous (see Section 1.2.4) because of the core width, that creates a shift in the distance. However, since a is very small compared to the minimum cell width ω min , i.e., ω min = 3L max , then both kernels almost coïncides in the farfield, as shown on Figure 6.4. Consequently, using the singular kernel in the farfield computation instead of the non-singular variant should only result in a minor extra error compared to the interpolation error. This assumption is validated in Section 6.4.2. On the other hand, as described in Section 1.1.3, given the homogeneity of the kernel, we can store the M2L operators at only one reference level and thus divide the memory footprint of the M2L step by the number of level L -1. We would like to highlight, that special care must be taken, when combining extended bounding boxes with the scaling of a homogeneous kernel. The kernel values are plo ed w.r.t. the scaled distance r/ω min with ω min = 3L max using a typical maximum segment length of L max = 90. Non-homogeneous kernels are represented for a = L max /30 = 3 (average core width) and a = L max /10 = 9 (large core width). Observa ons: The non-singular kernels coïncide with the singular kernel in the farfield, even for large values of the core width.

Numerical benchmarks

In this section, we present 2 artificial benchmarks in order to validate the accuracy our approach. We first verify the accuracy of the various Fast Multipole summation schemes presented in the previous Section 6.2 and then determine an appropriate interpolation order p, that leads to an accuracy of 10 -4 on the force field. Then, we discuss their respective sequential running times w.r.t. the number of segments using either the ufmm or the bbfmm.

Distributions of loops and visualization of the fields

In order to validate our implementation we do not need a very complex topology. However, it is crucial to respect the order of magnitude of the various parameters, that are involved in the dynamical simulation. In order to ensure the network to be closed and allow easy monitoring of the number of segments, we consider 2 very basic distributions of dislocation loops:

• the first one enforces 1 BCC loop per leaf (Fig. 6.5). As the tree grows, the diameter D and the minimum and maximum segment length decrease. More precisely,

D = 1000/2 L-2 and L = L min + 1 = L max -1 = 200/2 L-2 .
(6.13)

• the second uses a prescribed density d of loops in m/m 3 (Fig. 6.6). As the tree grows, the loop diameter D and segment lengths remain constant but the density grows.

Hence, in both cases the number of segment per leaf remains constant and the tree is always fully populated (worst case scenario) except for lowest tree depth L = 2.

Accuracy

Sources of error In order to verify the accuracy of the various schemes implemented during this thesis, we consider the convergence in term of the relative 2-norm error on the force and energy w.r.t. interpolation order p in the worst case scenario, where all possible sources of errors are taken into account. Therefore, we consider the distribution of loops driven by the specified density for the following reasons, since it is the closest to real life simulations, the network is closed and dislocation can cross the edges of the octree. Finally, the last source of errors is the use of the homogeneous kernel instead of the regularized kernel. Figure 6.7 represents the experimental accuracy of the ufmm and the bbfmm for a direct approach and the variant using shifted derivatives.

∥ f -f ∥ 2 /∥f ∥ 2 Nodal forces module cutoff ufmm -Direct bbfmm -Direct ufmm -Shifted bbfmm -Shifted

Discussion

The accuracy of the cutoff variant is represented on the graphs in order to show the benefit of having a farfield implementation, in this case the error on the force evaluation equals 4.10 -2 . The rate of convergence is similar to the scalar case. As expected, we observe a convergence delay of 1 interpolation order for certain values of p between the direct and the shifted cases. However, for p = 4 all schemes have about the same accuracy of 10 -4 . Therefore, in the range of accuracy targeted by the application, all schemes use the same p.

Sequential performance

Validation procedure First of all, in order to properly evaluate the relative performance of each variant we consider both configurations at comparable tree depths for p = 4. Then, in order to better show the asymptotic complexity of the method, we enforce a constant n 0 and a fully populated octree. Therefore, for a given distribution of loops, we first determine the optimal n for L = 2 and then multiply n by 8 each time we increment L. The sequential running times of the Factorized and the Shifted variants are represented on Figure 6.8.

Discussion

As expected, all variants have a linear complexity with the number of segments, however they do not evaluate the farfield with the same efficiency. The simulations driven by the specified density reach their permanent regime a little slower, since for the first densities the octree is not fully populated. It is important to mention, that by fixing the tree depth we cannot ensure that all variants balance nearfield and farfield computations optimally. In particular, in our numerical benchmarks, only the Symmetric variant with bbfmm is balanced. Nevertheless, the cheapest approach is the Shifted variant with ufmm, where the cost of the farfield computation is 3 times lower than for the Factorized variant.

(3)

(4) (5) (6)
# segments/nodes: n computational time (s)

1 loop per leaf cell

(3)

(4) (5) (6)
# segments/nodes: n computational time (s)

Density of loops

Part III

Randomized low-rank approximations of covariance matrices 7.1 Fundamentals and computational aspects

Generating realizations of Gaussian Random Fields (GRFs) is a key issue in numerous scientific research fields such as cosmology [START_REF] Pen | Generating cosmological gaussian random fields[END_REF][START_REF] Bertschinger | Multiscale gaussian random fields and their application to cosmological simulations[END_REF][START_REF] Carron | On fast generation of cosmological random fields[END_REF], geostatistics [START_REF] Andre | Geostatistics for conditional simulation of ore bodies[END_REF][START_REF] Le Ravalec | The fft moving average (fftma) generator: An efficient numerical method for generating and conditioning gaussian simulations[END_REF], hydrogeology [START_REF] Mjl Robin | Cross-correlated random field generation with the direct fourier transform method[END_REF][START_REF] Bellin | Hydro_gen: A spatially distributed random field generator for correlated properties[END_REF] or even Brownian dynamics [START_REF] Adolfo | Accelerated stokesian dynamics: Brownian motion[END_REF]. A GRF is a multivariate Gaussian Random Variable (GRV), i.e., a vector Y ∈ R n of correlated Gaussian random variables each associated with the points of a grid denoted {x i } i=1,...,n . In the present work we will only consider spatial grids, i.e., x i ∈ R 3 .

Basics

The concept of GRF is extensively presented in [1], in this section we will recall its fundamental concepts and discuss several properties of spatial correlation functions and covariance matrices.

Multivariate Random Variables

In order to accurately represent the first moments of a statistical distribution, namely the mean and the variance, GRFs usually need to be sampled in large ensembles (laws of large numbers). For instance, let X be a univariate GRV with mean 0 and variance 1, as shown on Table 7.1, one need to perform about 10 3 realizations X i in order to reach 1% accuracy in term of sample mean and variance, namely

E(X) = 1 n real n real ∑ i=1 X i and E(X 2 ) = 1 n real n real ∑ i=1 X 2 i
In the multivariate case, let say a length-n GRF, i.e., n GRVs, if n = 10 3 one needs n real = O(10 5 ) realizations to reach errors between 1 and 10% in terms of sample mean E(Y) and covariance E(YY T ). Therefore, efficiently generating many realizations of a GRF given a large spatial grid, i.e., n = O(10 6 ), can rapidly become challenging even on modern computers. i.e.,

n real E(X) -0 roc E E(X 2 ) -1 roc E 2 -1 t X (s) 1 • 10 2 7.68 • 10 -2 0.26 1.62 • 10 -5 1 • 10 3 1.8 • 10 -2 0.
C = {k(r ij )} i,j=1,...,n , where r ij = ||x i -x j || 2 .
The kernel k denotes a correlation function such as the Gaussian kernel

k ∞ (r) = e -r 2 /(2ℓ 2 ) (7.1)
or the exponential kernel

k 1/2 (r) = e -r/ℓ (7.2)
where the length scale ℓ characterizes the decreasing speed of the correlation function. These functions are the extreme cases of Matérn functions [START_REF] Matérn | Spatial variation[END_REF] defined as

k ν (r) = 2 1-ν Γ(ν) ( √ 2νr ℓ ) ν B( √ 2νr ℓ ), (7.3)
where B denotes a modified Bessel function and Γ the gamma function [2]. In particular, as shown on Figure 7.2, k ∞ is globally smooth and k 1/2 has discontinuous first order derivative at the origin. In this chapter, we also consider k 5/2 , which presents an intermediate degree of smoothness and is often used in machine learning applications. Covariance matrices are symmetric positive definite (spd) by definition of correlation kernels [1]. Hence, C admits the following representation

C = AA T , ( 7.4)
where the matrix factor A ∈ R n×n will be called a square root of the covariance matrix C, even though others may call it a Cholesky factor and would expect the square root of an spd matrix to be spd itself.

Sampling from GRF Realizations of a GRF Y ∼ µ(0, C) can be obtained in O(n 2 ) by applying A to a length-n white noise X, i.e., by computing Y = AX, (7.5) where X ∼ µ(0, I n ) is a Gaussian Random Field verifying E(X) = 0 and E(XX T ) = I n .

Limits of standard and alternative techniques

There exist numerous approaches for generating GRFs, that usually differ by the way A is precomputed. Approaches based on standard matrix decompositions such as the Singular Value Decomposition (SVD) or the Cholesky Decomposition [START_REF] Michael | Production of conditional simulations via the lu triangular decomposition of the covariance matrix[END_REF] are the most popular ones, since they provide an exact factorization and rely on well understood and robust algorithms that apply to any spd matrix C. If the decomposition can be truncated at a prescribed numerical rank r, then the resulting low-rank factor A r ∈ R n×r can be applied to a length-r white noise in O(n) operations in order to generate a GRF Y with mean 0 N and covariance close to C, i.e., E(YY

T ) = A r A T r ≈ C.
Alternatives Nevertheless, standard factorization algorithms involve O(n 3 ) operations and thus become computationally prohibitive for large n, i.e., n over a few thousands.

Alternative methods are often considered such as the sequential simulation method [START_REF] Mark | Multivariate statistical simulation: A guide to selecting and generating continuous multivariate distributions[END_REF], the moving average methods [START_REF] Dean | Moving averages for gaussian simulation in two and three dimensions[END_REF], the turning bands method [START_REF] Andre | Geostatistics for conditional simulation of ore bodies[END_REF][START_REF] Mantoglou | The turning bands method for simulation of random fields using line generation by a spectral method[END_REF], continuous [START_REF] Shinozuka | Digital simulation of random processes and its applications[END_REF] or discrete [START_REF] Gutjahr | Fast fourier transform methods for random field generation[END_REF][START_REF] Mjl Robin | Cross-correlated random field generation with the direct fourier transform method[END_REF] spectral methods, combinations of spectral methods with the turning bands [START_REF] Mantoglou | Digital simulation of multivariate two-and three-dimensional stochastic processes with a spectral turning bands method[END_REF] or with moving average [START_REF] Le Ravalec | The fft moving average (fftma) generator: An efficient numerical method for generating and conditioning gaussian simulations[END_REF], improved matrix decomposition [START_REF] Michael | Generating large stochastic simulations-the matrix polynomial approximation method[END_REF]. They usually provide approximate square roots but most importantly they do not always extend to 3D grids.

Circulant Embedding A popular alternative based on matrix decomposition, known as circulant embedding technique, was developed in the early 90's [START_REF] Dietrich | A fast and exact method for multidimensional gaussian stochastic simulations[END_REF][START_REF] Andrew | Simulation of stationary gaussian processes in [0, 1] d[END_REF]. It provides an exact method for computing the product Y = AX, that has a O(n) cost in both running time and storage, but most importantly A is assembled in Fourier domain in O(n log n) operations using the Fast Fourier Transform (FFT). However, it presents significant numerical limitations, especially in 3D, not to mention that it only applies to equispaced grids.

Fast randomized algorithms for generating Gaussian Random Fields

Based on the observation that many covariance matrices of practical interest actually have low-rank, Dehdari and Deutsch [START_REF] Dehdari | Applications of randomized methods for decomposing and simulating from large covariance matrices[END_REF] used the Randomized SVD (rSVD) in order to accelerate the precomputation of A in low-rank form and thus efficiently generate realizations of GRFs.

In this section, we present the details of this randomized algorithm and its numerical limitations. Then, we show that it can be significantly improved by using a H 2 -representation of the covariance matrix such as the Fast Multipole Method (FMM).

Randomized approach

Correlation functions used in practical applications usually lead to globally low-rank covariance matrices, therefore the rSVD (Algorithm 6) can be considered in order to efficiently approximate C by its approximate low-rank representation

C = Ũ Σ ŨT . (7.6)
Since covariance matrices are by definition symmetric positive semi-definite (spsd), a lowrank matrix factor can then be computed as

à = Ũ Σ1/2 ∈ R n×r . (7.7)
In fact, since à provides an exact square root of C, i.e., C = à ÃT , then à provides an approximate square root of C, i.e., C ≈ à ÃT . The matrix à can be applied to a length-r white noise at a O(n × r) cost, which allows for generating realizations of n-variate Gaussian random variables at a linear cost in time. However, this approach still requires the input covariance matrix C to be fully assembled and the cost of precomputing an approximate square root scales quadratically with n. 2.000 particles on the unit sphere 2.000 particles on the unit sphere

k 1/2 , ε = 10 -1 k ∞ , ε = 10 -1 k 1/2 , ε = 10 -2 k ∞ , ε = 10 -2
k 3/2 , ε = 10 -1 k 5/2 , ε = 10 -1 k 3/2 , ε = 10 -2 k 5/2 , ε = 10 -2 Figure 7
.3: Numerical rank r of some covariance kernel matrices w.r.t. the correla on length scale ℓ for 1% and 10% target accuracy in Frobenius norm. We used n = 2.000 par cles distributed on the unit sphere. Observa ons: In this range of accuracy the rank grows quadra cally with the inverse of ℓ, while being inferior to 10 for ℓ = 1. Finally, the rank never exceeds 40% of the full size except for the exponen al kernel k 1/2 with ℓ = 1.0.

Numerical rank

The complexity of the method grows linearly with the rank of the matrix, therefore its applicability is restricted to matrices with a relatively low rank. Before describing our new algorithm we would like to evaluate if the kernel matrices we want to approximate fulfill this condition. For a given accuracy and a given norm, the rank r can be determined numerically by considering the best rank-r error computed from the spectrum of Gaussian k ∞ on the unit sphere Effect of the oversampling s alone (q = 0). Observa ons: The Gaussian kernel is smooth, therefore its spectrum decreases fast and C is rela vely low-rank. Consequently, the rSVD performs well even without power itera ons and with low oversampling (s = 5). An oversampling of s = 50 leads to a near op mal error. Note that C becomes high-rank for small ℓ, e.g., ℓ = 0.1.

∥C -C r ∥ {S,F },rel ∥C -C r ∥ S = |σ r+1 | ∥C -C r ∥ F = |σ 2 r+1 + . . . + σ 2 n | 2 4 6 8
∥C -C∥ F,rel ∥C -C r ∥ F,rel s = 5 -q = 0 s = 50 -q = 0
the input matrix, we call r the numerical rank. The target accuracy required for generating GRFs is usually low, namely between 10 -1 and 10 -2 in either Frobenius or Spectral norm. In this thesis we chose to use the Frobenius norm to determine the numerical rank, since it is more restrictive than the Spectral norm and it better accounts for the shape of the spectrum. In order to allow the computation of the full spectrum within a reasonable time by mean of a full SVD, we chose a relatively low value of n, namely n = 2.000. In fact, as shown by our numerical benchmarks, for sufficiently dense distributions, i.e., n > 1.000 for the unit sphere, the rank of the kernel matrix does not depend on n but rather on the kernel function and the geometry. Oversampling Figure 7.4 illustrates the benefit of using oversampling with the Gaussian kernel k ∞ (7.1) for n = 2.000 particles distributed on the unit sphere. No power iteration is required, since the associated covariance matrix has a fast decreasing spectrum. Hence, we set q = 0 and analyze the effect of the oversampling. First of all, the results show that the fixed rank approach works well for matrices of relatively low rank, e.g. for ℓ ≥ 0.50. For instance, if ℓ ≥ 0.50, r/n ≈ 5% and s = 5, then the relative error in Frobenius norm is already below 2%. On the other hand, for ℓ = 0.1 the matrix is not low-rank, since for a compression rate of 10% the optimal relative error exceeds 10%. Moreover, as shown on Figure 7.4, increasing the oversampling significantly improves the accuracy and even provide near-optimal accuracy on this particular test case. In fact, for ℓ ≥ 0.5, the error for s = 50 almost coïncides with the theoretical lower bound ∥C -C r ∥ F .

Subspace iterations

On the other hand, Figure 7.6 illustrates the benefit of using subspace iterations with the exponential kernel k 1/2 (7.2), that exhibits a relatively flat spectrum. In particular, it shows that with an oversampling of s = 50 the accuracy is still far from the optimal. However, a near-optimal accuracy can be reached with a few subspace iterations, namely q = 1. Since in this case the matrix is not low-rank, it only produces a raw approximation of C in the desired range of compression, i.e., compression rates below 10%. Hence, we will not consider this kernel in the numerical benchmarks presented in Section 7. rank: r best rank-r error: Effect of q power itera ons using a low oversampling, namely s = 5. Observa ons: The exponen al correla on k 1/2 is not smooth, therefore its spectrum is rela vely flat. One power itera on leads to a near-op mal error but the op mal error remains high. In fact, a 1% rela ve error is not reachable for ℓ < 0.5 for a maximum compression rate of 10%.

∥C -C r ∥ {S,F },rel ∥C -C r ∥ S = |σ r+1 | ∥C -C r ∥ F = |σ 2 r+1 + . . . + σ 2 n |
∥C -C∥ F,rel ∥C -C r ∥ F,rel s = 5 -q = 0 s = 5 -q = 1

FMM-powered randomized algorithms

Algorithm & Theoretical complexity In the rSVD, the matrix-to-matrix multiplications involved at each stage, namely CΩ for the random projection (3.4) and CQ for the assembly of B = Q T CQ (3.6), dominate the asymptotic cost of the algorithm. Therefore, the overall cost of the original randomized SVD is quadratic in n in both running time and memory footprint. Since C is given as a kernel matrix it can be applied to another matrix at a linear cost in n without ever assembling the full matrix using a H 2 -method like the ufmm or the bbfmm. The resulting approach remains a dense random projection as Ω stays dense, but it benefits from a data sparse representation of C. Algorithm 8 presents our H 2 -powered randomized SVD and indicates the asymptotic complexity of each stage. The resulting algorithm has an overall O(r 2 × n) computational cost, which is one order less than the O(r × n 2 ) complexity of the original method. However, we will show that for the rank: r best rank-r error: Observa ons: The correla on is smoother than the exponen al, therefore its spectrum decreases faster. On one hand, it seems insufficient to increase the oversampling by a constant value, namely s = 5 or s = 50. On the other hand, one power itera on leads to a near-op mal error but the op mal error remains rela vely high. In fact, a 1% rela ve error is not reachable for ℓ < 0.5 for a maximum compression rate of 10%. grid size and accuracy used in practice, the matrix multiplications still dominate the cost although they have a O(r × n) complexity.

∥C -C r ∥ {S,F },rel ∥C -C r ∥ S = |σ r+1 | ∥C -C r ∥ F = |σ 2 r+1 + . . . + σ 2 n | 2 4 6 8
∥C -C∥ F,rel ∥C -C r ∥ F,rel s = 5 -q = 0 s = 50 -q = 0 s = 5 -q = 1

Algorithm 8: H 2 -powered Randomized SVD

Input: Correlation kernel k(., .), grid size n, positions x ∈ (R 3 ) n , rank r or accuracy ε, number of power iterations q, oversampling parameter s, interpolation order p, octree depth h Output:

[ Ã] approximate square root of C = {k(x i , x j )} i,j=1,...,n // Stage 0: Build H 2 representation of C Build H 2 C = H 2 (k, x, N, p, h) // Stage I: Approximate the range of C if accuracy ε is prescribed then [Q, r] = ARRF (H 2 C, q, s, ε) O(r 2 × n) if rank r is prescribed then [Q, ε] = RSI(H 2 C, q, s, r) (idem) // Stage II: Decompose C r = Q(Q T H 2 CQ)Q T ≈ H 2 C as UΣU T Build B = Q T (H 2 CQ) ∈ R r×r O(r 2 × n) Perform SVD of B = U B Σ B U T B O(r 3 ) // Compute square-root à = C1/2 ≈ C 1/2 as Ũ Σ1/2 B 1/2 = U B Σ 1/2 B O(n × r 2 ) à = QB 1/2 O(n × r 2 )
Accuracy The FMM provide approximate matrix multiplications, therefore it introduces extra errors in each stage of the algorithm. The accuracy of the FMM needs to be carefully the compression rate for n = 2.000 par cles distributed on the unit sphere. Le : For a Gaussian kernel with ℓ = 0.5 the accuracy of the rSVD can be very high even below 10% compression. However, the ufmm only allows for represen ng the covariance matrix up to a fixed accuracy. For a given interpola on order p, the H 2 powered variant has the same accuracy as the dense variant then slightly reaches a plateau. Right: For the exponen al kernel, even with ℓ = 1.0 the accuracy of the rSVD is rela vely low, therefore the ufmm never disturbs the method in this range of compression and p can be set to the lowest value.

monitored in order to avoid perturbing the approximation of the range and the final accuracy of the rSVD. In particular, the optimal order of interpolation should not affect the output rank of the Adaptive Randomized Range Finder (ARRF) or the accuracy of the Randomized Subspace Iterations (RSI). Figure 7.7 shows the effect of ufmm on the accuracy of the rSVD for We found that enforcing one order of magnitude between the error on a fast matrix multiplication ε f mm = ∥ φ -Kw∥ 2,rel and the prescribed accuracy ε of the randomized algorithm, i.e., ε = 10ε f mm , (7.8) preserves the accuracy of the original method, while still providing significantly better performance. Hence, in the context of GRF, the interpolation order of the ufmm never exceeds p = 4.

Efficient Fast Multipole matrix-to-matrix multiplication

In this section, we present an optimized method for performing efficient Fast Multipole matrix-to-matrix multiplication, i.e., FMM with multiple right-hand-sides. In particular, we describe the efficient implementation of the interpolation operators (P2M, M2L & L2P) and the P2P operators in order to minimize the cost of the algorithm in terms of both storage and running times as the number of right-hand-sides grows.

Blocking farfield computations

Generic notations Handling multiple right-hand-sides in a FMM computation consists in computing multivariate potentials using a scalar interactions kernel and multivariate densities. Let r denote the number of right-hand-sides, then the i-th component of the potential generated at the target point x by the i-th component of the density w located at the source point y reads as ϕ i (x) = k(x, y)w i (y), ∀i = 1, . . . , r. (7.9)

The interpolation formula used for the farfield computations reads as

φfar i (x) = ∑ |α|≤p S α (x) ∑ |β|≤p
Kαβ S β (y)w i (y) (7.10)

Summation scheme

The fast multipole summation consists in performing the following steps:

• P2M: Transfer the contributions of the r components of the density w located at the source points y into the r multipole expansions (M β ) i defined on the source interpolation nodes ȳβ

(M β ) i = ∑ y S β (y)w i (y), ∀i = 1, . . . , r (7.11)
• M2L: Transfer the r multipole expansions (M β ) i into the r local expansions (L α ) i defined on the target interpolation nodes xα

(L α ) i = ∑ |β|≤p Kαβ (M β ) i , ∀i = 1, . . . , r (7.12)
• L2P: Transfer the contributions of the r local expansions (L α ) i defined on the target interpolation nodes xα into the r components of the potential φ located at the target points

x φi (x) = ∑ |α|≤p S α (x) (L α ) i , ∀i = 1, . . . , r (7.13)
Memory footprint Although it does not require extra work or memory during precomputation, this summation scheme presents strong limitations in term of storage, in particular when it comes to storing the expansions at run time. Indeed, the r multipole expansions, r local expansions, and the associated transformed expansions have to be stored independently. For instance, as shown on Figure 5.6, for p = 4 and a densely populated octree of maximum level L = 6, the memory required to store the expansions associated to one right-hand-side equals 5 GigaBytes. Therefore, the number r of right-hand-sides should not exceed r max = 10, because this would result in a total memory footprint of 50GB for the expansions only. As a comparison, storing the entire matrix would cost 8 × (10 5 ) 2 = 80GB. Consequently, if the input r exceeds r max , e.g., r is a fraction of the problem size, then we need to perform several sequential but independent multi-rhs FMM. Computational cost In a naïve implementation, the computational cost of this scheme would roughly equal r times the cost of a standard FMM. First of all, our implementation of the P2M and L2P steps decouples the computation of the polynomial interpolators and the multivariate expansions, in order to factorize as much operations as possible. Second of all, the M2L step is implemented in a similar fashion as the tensorial approach presented in Appendix A. In the bbfmm, if r were large enough, the M2L step could be performed in a block way instead of component by component, e.g., through a matrix-to-matrix product based on level-3 Blas routines. However, due to the excessive memory required for large r this approach is not relevant. In the ufmm, there is no point in blocking the operations, since we perform an entry-wise product. Consequently, at the M2L step in both variants, we transfer one rhs at a time using a simple loop over r max . Hence, transferring all rhs costs r times more than transferring one rhs.

K L 1 L 2 . . . L r M 1 M 2 . . . M r M2L S w 1 w 2 . . . w r P2M S T ϕ 1 ϕ 2 . . . ϕ r L2P

Precomputing the P2P operators

As we just demonstrated, there is not much room for improvement in the farfield computation when performing FMM with multiple rhs, since we quickly reach a memory limit. On the other hand, as explained in this subsection, many operations can be factorized during the nearfield computation.

Implementation

The nearfield computations of a single rhs FMM involve the computation of the kernel for all pairs of particle lying in neighbor cells and the direct evaluation of the sum (7.9). However, there is no need to compute a new set of interactions for each right-hand-side, the P2P operators can instead be precomputed and stored. The P2P opera-tors are eventually applied to a set of r rhs at once using Blas::gemm() routines provided by the Intel MKL library. As a result, the application of the P2P is a lot faster than in the original algorithm and the balance is not necessarily obtained for the same tree depth L.

Storage

The memory required for storing the P2P operators

K = { k(x i , y j ) } x i ∈tree,y j ∈N (C ( L) x i ) (7.14)
equals n × 27n 0 × 8 bytes in double precision arithmetic if the tree is fully populated. Since the memory footprint scales like n, i.e., like the number of (non-empty) cells, then it can quickly become prohibitive. Figure 7.9 shows the theoretical amount of memory needed for storing K for a uniform distribution of particles in a cube for various values of n and n 0 . If we consider typical simulations using n 0 = n/8 L = 64 particles per leaf, then for n = 256.000 particles with L = 4, we would already need 3.5GB and for n = 2 × 10 6 particles with L = 5, we would need about 30GB.

Numerical Benchmarks

In this section, we present the performance of our algorithm within the range of accuracy required by the generation of Gaussian Random Fields. We first present the computational running time of our multi-rhs FMM for a smooth kernel. Then, we compare the performance of the H 2 -powered rSVD with the rSVD and the truncated SVD. Finally, we validate our method by generating many Gaussian Random Fields and verifying the convergence of the sample mean and covariance w.r.t. the number of realizations. We also present realizations on the unit sphere with the Gaussian kernel for multiple length scales.

Fast Matrix multiplication multi-rhs

In order to calibrate our multi-rhs FMM, we measure the performance of the ufmm and smooth-ufmm for p = 3, which ensures an accuracy of ε f mm = 10 -3 on the matrix multiplications. Figure 7.10 shows the running times per matrix-to-vector product for r = 10 rhs and particles distributed in the unit cube or on the unit sphere. Since the experiments are very similar to those presented in Section 5.4, we draw similar conclusions concerning the relative performance of the ufmm and the smooth-ufmm, as well as concerning the good performance on heterogeneous grids such as the unit sphere. The computational time are significantly lower since the interpolation order is lower and the P2P operators are precomputed. In order to show the speed-up of the multi-rhs FMM compared to the direct computation, we assemble a full covariance matrix and apply it to r = 10 rhs using MKL Blas routines.

global fft The smooth-ufmm with L = 0 boils down to an FFT-accelerated Lagrange interpolation on the root bounding box. We call this method the global fft as it generalizes the idea of the FFT on an arbitrary grid. The global fft has a theoretical linear complexity in n, that can be verified on Fig. 7.10. Furthermore, for a given root bounding box the performance of the global fft are independent of the shape of the distribution. Although this method requires larger interpolation orders to reach similar accuracies as the ufmm and the smooth-ufmm, that both rely on several subdomains, only one expansion has to be transfered. As shown on Fig. 7.10, even in this range of accuracy, the global approach is always slower than the hierarchical variants even with a fully populated octree (Fig. 7.10 left). Moreover, since hierarchical variants benefit from the heterogeneity of the distribution, they beat the global fft by a significant factor on geometries such as the unit sphere (Fig. 7.10 right). The figures presented in Appendix C confirm these observations for a fixed n and a varying accuracy. However, they show that the global fft exhibits performance that are similar to the hierarchical variants for high values of ℓ, i.e., very smooth kernels.

Fast Randomized SVD

The accuracy and the O(n) complexity of the H 2 -powered rSVD are now illustrated for a Gaussian correlation with ℓ = 0.5 and particles distributed on the unit sphere.

Optimal setup As shown by Figure 7.4 the covariance matrix can be represented by a matrix of rank r ≈ 70 with a precision of about ε = 10 -2 in Frobenius norm using a conventional rSVD. In this case, our experiments showed that the smallest interpolation order that does not affect the randomized algorithm is p = 3 for the ufmm, i.e., ε f mm < 10 -3 . If the desired accuracy equals ε = 10 -3 then r ≈ 100 and p = 4 was found optimal, i.e., ε f mm < 10 -4 . This confirms the optimal setup rule proposed in (7.8). )

O (n ) (2) (3) (4) (5)
# particles: n running time (s) Cube 10 3 10 4 10 5 10 6

10 -3

10 -2 10 -1 10 0 10 1 10 2 O ( n 2 ) O (n ) (2) (3) (4) (5) (6)
# particles: n running time (s) Complexity & Speed-up Figure 7.11 shows the time required to build an approximate square root of C within an accuracy of 10 -2 in Frobenius norm. It confirms the theoretical asymptotic complexities of all approaches, namely cubic for the full SVD (in red), quadratic for the truncated SVD based on Arnoldi's algorithm (in green) and randomized SVD (in black) and finally linear for the H 2 -accelerated randomized SVD. As mentioned in [START_REF] Halko | Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions[END_REF], we observe that the running times of the fixed rank and fixed accuracy variants are very close to each other and that they outperform the truncated SVD by a factor of 6.5. Moreover, the H 2 variants (ufmm and smooth-ufmm) are faster than the original rSVD for n ≥ 10 4 . Not only do these variants allow computation to be performed for n up to a few millions but they also provide significant speedup, e.g., they would be about 10 times faster for n = 10 5 if only a dense matrix multiplication was affordable. This acceleration is even more pronounced for smooth kernels and heterogeneous grids as explained in Section 5.4 and recalled by the results displayed on Figure 7.10. Finally, we notice that even if the cost of computing matrix multiplication has been dramatically reduced, it still governs the overall cost of the algorithm. For instance, for n = 2.592.000 particles and r = 70 the random projection is done in about 360s, while the QR Decomposition only takes 50s. This can be seen directly on Figure 7.10 since the curves corresponding to FMM algorithms have the same shape than those of Figure 7.11.

Realizations of Gaussian Random Fields

Gaussian Random Fields were simulated on particles distributed on the unit sphere using Gaussian correlations with varying length scales ℓ. The approximation is done by mean of a fixed precision (ε = 10 -2 ) rSVD with ufmm-or smooth-ufmm accelerated matrix multiplications. Realizations are displayed on Figure 7.12 for length scales ℓ ∈ {0.25, 0.5} using the smooth-ufmm. The figure shows that the actual correlation length between grid points roughly corresponds to the input ℓ. Smaller length scales computed by mean of a standard Observa ons: As expected, both fixed rank and fixed accuracy variants exhibit similar performance. On the other hand, the graphs confirm the theore cal asympto c costs (linear in blue, quadra c in black and cubic in red). Machine: plafrim/riri -Deca-core Intel Xeon E7-4870 @ 2.40GHz with 1TB ram and 30MB L3 Cache.

ufmm are displayed on Figure 7.13. In particular, it shows that decreasing the accuracy by one order of magnitude is hardly noticeable on such visualizations. Furthermore, since the length scale is quite low, namely ℓ = 0.1, the rank is relatively large but the square root is still precomputed in a reasonable amount of time, namely 90 seconds for 10 -1 accuracy and 175 seconds for 10 -2 .

Accuracy

The sample covariance matrix C real , computed from the realizations as

C real = 1 n real n real ∑ i=1 (Y i -E(Y i ))(Y i -E(Y i )) t (7.15)
provides a good approximation of the experimental covariance of a finite sample. Therefore, we analyze the accuracy of the method using the error between C real and the input covariance matrix C. In order to limit the computational cost required by verifications, the error can be computed on a random subset of the matrices. The accuracy of the method w.r.t. the number of realizations n real is presented in Table 7.1. These results show that the ufmm and the smooth-ufmm can be used to generate correlated random fields within similar accuracies for ℓ ≥ 0.25, while maintaining a convergence rate that is close to the theoretical 1/2.

For lower length scales we recommend using the ufmm in order to better capture the fast variation of the kernel, i.e., the potentially high rank of the matrix.

Running times As shown on Figure 7.11, for a Gaussian correlation with ℓ = 0.50, building an approximate square root with n = 72.000, r = 70 takes about 20 seconds with the ufmm, 8 seconds with the smooth-ufmm while it would take about 5 minutes with dense matrix multiplications. For ℓ = 0.25, the rank r and the computational times are about 3 times larger. The ufmm also provides significant speed-up for matrices of relatively high rank associated with low length scales ℓ or non-smooth correlations. For a given length scale, the performance of the ufmm depends only very weakly on the nature of the correlation function, therefore the computational cost associated with non Gaussian correlation can be easily estimated once the rank is known.

Gaussian k ∞ ℓ = 0.25, r = 210 ℓ = 0.50, r = 70 .000 points distributed on the unit sphere using a Gaussian correla on. The approximate square root was obtained by mean of a ufmm-accelerated rSVD with a fixed accuracy of ε = 10 -2 (le ) and 10 -1 (right). Observa ons: The differences in terms of covariance quality and smoothness of the field are hardly no ceable on a 3D visualiza on. However, since the rank is roughly divided by 2, the method is twice as fast for the lower target accuracy.

Conclusion

Square root algorithm We designed a matrix-free algorithm for computing approximate square roots of kernel covariance matrices in SVD form in linear time. This algorithm relies on a new highly efficient kernel independent FMM that can further benefit from the heterogeneity of the grid and the smoothness of the correlation kernel. The presentation of the algorithm and the experimental results are summarized in a research report [START_REF] Blanchard | Fast hierarchical algorithms for generating gaussian random fields[END_REF]. Many geostatistical applications and in particular Kalman Filtering, a very popular data assimilation technique, involve intensive matrix computations on arbitrary grids with up to several million points. Fast low-rank variant of those algorithms often require the efficient computation of matrix square root, namely Ensemble Kalman Filters. Another useful application of our algorithm called denoising consists in removing correlated noise from a given random field using the inverse square root of the covariance matrix.

ℓ = 0.25 ufmm smooth-ufmm n real ∥C real -C∥ 2,rel roc 2 ∥C real -C∥ ∞,rel roc ∞ 1 • 10 3 2.48 • 10 -1 1.47 • 10 -1 1 • 10 4 7.
n real ∥C real -C∥ 2,rel roc 2 ∥C real -C∥ ∞,rel roc ∞ 1 • 10 3 1.31 • 10 -1 1.12 • 10 -1 1 • 10 4 3.
H 2 -powered Randomized SVD Our approximate SVD algorithm showed good performance even for matrices of relatively high rank and can partially overcome issues related to matrices with relatively flat spectrum. This contribution is a rare attempt at accelerating random projection based Low-Rank Approximation (LRA) algorithms using the block low-rank structure of the input matrix, which proved its efficiency on artificial test cases and should now be applied to real-life applications. Finally, since many covariance matrices arising in geostatistical applications are computed purely algebraically, e.g., from the experimental field, one should consider using the general H 2 format in order to derive linear time factorization algorithm.

Covariance matrix computations are involved in numerous fields of scientific computing such as geostatistics, computational economics, computational biology, machine learning,… In geostatistical applications such as the one addressed in Chapter 7 the covariance matrix is often given as a smooth and fast decreasing kernel with nice low-rank properties, that allows for Fast Multipole Method (FMM)-powered matrix computations, e.g., multiplication, matrix decomposition, square root, inversion,… In the general case where the covariance matrix is only known algebraically, e.g., assembled from the data, it often exhibits a certain structure, that can be exploited in order to improve the performance of the matrix computations. In this chapter, we focus on a dimensionality reduction algorithm known as Multidimensional Scaling (MDS), that allows for visualizing a point cloud in a low dimensional subspace by the only knowledge of the set of pairwise distances. In particular, we want to make this method tractable to large biological datasets, where distances are computed as dissimilarities between DNA sequences. This application was addressed for the following reasons:

• Extend the capabilities of our randomized LRA library to arbitrary covariance matrices.

• Design an operational tool for performing classification on large real life datasets provided by New Generation Sequencing (NGS).

• Provide new materials and design new metrics for investigating new ways of characterizing the biodiversity.

We first describe the standard MDS approach as well as some specificities of biological datasets. Then, present a new approach based on random projection and discuss its complementarity with the state-of-the-art method. Finally, we provide some visuals in order to better understand what useful outputs our method can provide to perform classification on large datasets.

Multidimensional Scaling on biological datasets

In this section we present the fundamentals of MDS and discuss some key issues related to error estimation. Then, we put the emphasis on the computational aspects and we highlight the limits of the standard fast MDS algorithm. Finally, we discuss some specific issues related to the operational chain addressed in this thesis.

Multidimensional Scaling

The method implemented here is the metric MDS, please refer to [START_REF] Izenman | Modern Multivariate Statistical Techniques[END_REF] for a recent survey and [START_REF] Cox | Multidimensional Scaling -Second edition[END_REF] for a seminal monograph.

Metric MDS Let us have a set of n items i with i ∈ E = {1, . . . , n}, and a distance D ij between items i and j, then (E, D) is a finite metric space. Let r ∈ N and R r have a bilinear form B with an associated quadratic form q. Then, (E, D) is embedded into (R r , B) if there exists a map

E : i ---→ X i ∈ R r such that D 2 ij = q(X i -X j ). (8.1)
If (R r , B) is an Euclidean space, this is the classical

D 2 ij = ∥X i -X j ∥ 2 2 . ( 8.2)
For a given dimension r ∈ N, the metric MDS consists in finding a map such that ∥X i -X j ∥ 2 is as close as possible to D ij . The solution is well-known when B is definite positive [START_REF] Cox | Multidimensional Scaling -Second edition[END_REF] and it is implemented in three steps:

1. Pre-treatment: We denote S, the similarity matrix associated with B,i.e.,

S ij = B(X i , X j ) (8.3)
whose entries can be computed from the distances only as

B(X i , X j ) = - 1 2 ( D 2 ij - 1 n ∑ k D 2 kj - 1 n ∑ k D 2 ik + 1 n 2 ∑ k,ℓ D 2 kℓ ) (8.4)
2. Factorization: Compute the eigenvalue decomposition of S S = UΛU T (8.5)

3. Components: Represent the point cloud by the components of

X = UΛ 1/2 (8.6)
It is straightforward that X verifies S = XX T , in fact

S = UΛU T = (UΛ 1/2 )(UΛ 1/2 ) T = XX T (8.7)
If all eigenvalues are strictly positive, then (E, D) can be embedded isometrically in a Euclidean space, and a best low-dimensional approximation can be computed. The accuracy of this map is usually quantified by the error committed on the distance reconstruction, which reads as

δ(X) = ∥∆ 2 (X)∥ F ∥D 2 ∥ F , ( 8.8)
where the entries of the distortion matrix ∆ 2 (X) read as

∆ 2 ij (X) = D 2 ij -∥X i -X j ∥ 2 2
, ∀i, j = 1, . . . , n. (

Negative eigenvalues However, in general, not all eigenvalues of S are strictly positive, since the input distance is not Euclidean. This means that √ Λ may have imaginary parts, which may induce extra distortion in the point cloud and affect the quality of the reconstruction. In fact, it is possible to show that there exists an isometry into a pseudo-euclidean space, i.e., a vector space with a bilinear symmetric form, with positive and negative eigenvalues (see [START_REF] Goldfarb | A unified approach to pattern recognition[END_REF] and [START_REF] Pekalska | The dissimilarity representation for pattern recognition[END_REF] with references therein). Then, there exists two Euclidean spaces E + and E -of dimensions r + and r -respectively if the bilinear form is of signature (r + , r -), such that

D 2 ij = ∥X + i -X + j ∥ 2 2 -∥X - i -X - j ∥ 2 2
. Let Λ + be the diagonal matrix of strictly positive eigenvalues, and -Λ -of strictly negative ones. Let X + be the n × r + matrix of images in E + , and X -of images in E -. Then, X + = U + (Λ + ) 1/2 and X -= U -(Λ -) 1/2 . Hence, the distortion matrix can be expressed on both E + and E -as

∆ 2 ij (X + , X -) = D 2 ij -∥X + i -X + j ∥ 2 2 + ∥X - i -X - j ∥ 2 2 , ∀i, j = 1, . . . , n (8.10)
or on E + only as

∆ 2 ij (X + ) = D 2 ij -∥X + i -X + j ∥ 2 2 , ∀i, j = 1, . . . , n. (8.11)
Finally, we can show that for an exact representation of the similarity matrix, i.e., X i ∈ R n , the distortion associated to the representation in E + defined as

δ(X + ) = ∥∆ 2 (X + )∥ F ∥D 2 ∥ F (8.12)
is related to the energy of the matrix associated to negative eigenvalues ε -by the following expression

δ 2 (X + ) = ε - ε + -ε -, with ε -= ∑ r - i=1 (λ - i ) 2 ∑ n i=1 (λ i ) 2 , ( 8.13)
In the case of a rank-r approximation this relation is more complicated and will be analyzed experimentally in this chapter.

Other factorizations Many matrix decomposition techniques (Cholesky Decomposition, Singular Value Decomposition (SVD), …) can provide a valid factor X, i.e., real valued components, as long as the matrix is spsd. Since S is in general not strictly positive definite, the Cholesky Decomposition does not apply. However, S is symmetric, thus a factorization equivalent to (8.5) is achievable using an SVD, which takes the following form S = UΣV T (8.14) where V is also an orthonormal matrix and the singular values σ i need be positive and verify

σ i = |λ i |, ∀i = 1, . . . , n (8.15)
If λ i < 0 then the i-th column of V integrates this negative sign. Hence, if the factorization is obtained in SVD form then the conversion to Eigen Value Decomposition (EVD) form is straightforward.

Limits of standard techniques and alternatives

As far as calculation is concerned, the most demanding step is step 2, i.e., computing the eigenvectors and eigenvalues of S. As mentioned numerous times in this thesis, the cubic cost of direct approaches such as full SVD is intractable to large sample sizes. Although there exist iterative algorithms for computing EVD, e.g., Arnoldi, randomized low-rank approximation techniques usually perform better and offer more flexibility, while reaching near optimal accuracy for most matrices involved in practical applications.

Column selection

The usual work around in metric MDS is to consider subsets of the rows and columns of the input matrix, then perform the MDS on the resulting subsample and finally use a matrix factorization technique known as the Nyström method [START_REF] Williams | Using the nyström method to speed up kernel machines[END_REF] to reconstruct the entire point cloud. The standard Nyström approach consists in first sampling c columns from S and building an approximation of S as S ≈ CW † C T (8.16) where C ∈ R n×c gathers the sampled columns and W † is the Moore-Penrose pseudo-inverse of W the matrix of sampled rows and sampled columns. Hence, Nyström based MDS provide an approximate map, that reads as X ≈ CW †/2 (8.17)

Many fast MDS variants derive from that approach, e.g., Landmark MDS [START_REF] De | Global versus local methods in nonlinear dimensionality reduction[END_REF][START_REF] De | Sparse multidimensional scaling using landmark points[END_REF], FastMap [START_REF] Faloutsos | FastMap: A fast algorithm for indexing, datamining and visualization of traditional and multimedia datasets[END_REF] or MetricMap [START_REF] Wang | Evaluating a class of distance-mapping algorithms for data mining and clustering[END_REF]. More precisely, each approach differs from the other by the variant of the Nyström method it implements, namely the way it selects the columns and/or builds the factorization. Nyström based MDS is usually very convenient for biological applications since it offers interpretability, i.e., easy connection between components and individuals, and very competitive performance. However, it presents some significant drawbacks, that we discuss in the following.

Limitations First of all, the Nyström factorization only applies to spsd matrices. As stated in [START_REF] Gisbrecht | Metric and non-metric proximity transformations at linear costs[END_REF], although the extension of Nyström approach to non-spd similarity matrices was addressed in [START_REF] Belongie | Spectral partitioning with indefinite kernels using the nyström extension[END_REF], the alternatives usually lack stability. Moreover, the design of efficient non-negative embeddings is not the purpose of this work, but the reader will find a fair amount of articles on that topic under the name non-metric MDS. Second of all, as shown in [START_REF] Platt | Fastmap, metricmap, and landmark mds are all nystrom algorithms[END_REF], certain variant like Landmark MDS and FastMap may produce negative eigenvalues. Finally, it is usually very involved to predict or even just monitor the number of selected rows/columns to achieve satisfying accuracy on the rank r approximation of a given matrix, not to mention the variability of the approximation error. An improved Nyström factorization was proposed in [START_REF] Wang | Improving cur matrix decomposition and the nyström approximation via adaptive sampling[END_REF] and designed to provide better accuracy, that reads as S ≈ (CC † )S(CC † ) T (8.18)

However, as mentioned in [START_REF] Gittens | Revisiting the nystrom method for improved large-scale machine learning[END_REF], this variant is just another instantiation of the random projection based Low-Rank Approximations (LRAs) described in [START_REF] Halko | Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions[END_REF]. In particular, it has a similar computational cost without the many benefits of random projection described in the present chapter.

Random projection Although the concept of random projection is intrinsically related to dimensionality reduction through Johnson-Lindenstrauss Lemma [START_REF] William | Extensions of lipschitz mappings into a hilbert space[END_REF], to our knowledge there are very few contributions to MDS that involve neither sparse or dense random projection if any. Moreover, since random projection belongs to a broader class of sketching techniques, including random sampling to some extent, it offers a large variety of fast alternatives. In fact, many random projection techniques such as the original Fast Johnson Lindenstrauss Transforms (FJLT) [5], rely on sparse or very sparse projection [START_REF] Achlioptas | Database-friendly random projections: Johnson-lindenstrauss with binary coins[END_REF], which is in essence very close to column selection and performs similarly. On the other hand, dense random projection based LRAs [START_REF] Halko | Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions[END_REF][START_REF] Drinea | A randomized singular value decomposition algorithm for image processing applications[END_REF][START_REF] Vempala | The Random Projection Method[END_REF] represent significant alternatives to iterative EVD algorithms, as they share the same complexity, namely O(rn 2 ), but usually perform significantly better. Despite involving more intensive computations than sparse projection or random sampling they have many benefits such as better accuracy, robustness and low variability. Furthermore, they provide factorization in many practical forms such as SVD, Cholesky Decomposition or Interpolative Decomposition, and they can be tuned to reach near optimal accuracy at a relatively low extra cost. Finally, since the efficiency of the random projection based LRA algorithms mostly depend on the ability to apply fast matrix multiplication, it can be improved by exploiting the structure of the matrix and it can be parallelized straightforwardly.

Origin and nature of the data

Workflow We integrated our new algorithm in an operational chain that consists in the three following steps:

• Compute distance matrices from sequences provided by NGS and store them on iRods. Please refer to [START_REF] Frigerio | diagno-syst: a tool for accurate inventories in metabarcoding[END_REF] for further details on the pairwise distance computations.

• Transfer the distance matrix D associated with a sample to plafrim2, perform MDS on the full sample or subsamples and store associated components X.

• Post-treat components on a local or a cluster computer and provide some data to be analyzed, namely statistics and visuals.

We mostly contributed to the improvement of the second step, while focusing on providing relevant information for the third step that can be analyzed easily. In this section, we define the goals of our contribution in terms of performance and storage improvement, error estimation and visualization. In particular, we discuss the specificity of the data involved in this chain and illustrate it on small subsets of full samples for which exact full SVD is affordable.

Sample compute distances transfer and convert subsample perform MDS (Lake Geneva) S-W on Turing iRods→plafrim 10 4 reads full SVD 10 5 reads 4h 20min + 1h40 30s 20min

Table 8.1: Some average running mes related to expensive opera ons involved in the opera onal chain. Observa ons: The large amount of me required to convert distance matrices from text files to binaries actually represent a significant gain compared to the ini al state of the chain (text files read from python). In fact, since opera ons such as subsampling may be repeated numerous mes for a given sample, it is recommended to store the full distance matrix on the disk in a format that allows for a fast loading of the full sample in memory.

Samples

The samples studied in this thesis were collected in Lake Geneva at various time t = 1, . . . , 10 of the year, that we denote L t . They correspond to diatoms scraped from 5 stones each month between April 2012 and March 2013 in order to investigate a seasonal dynamics. Some individuals have previously been identified while other remain unknown.

Distances & Similarities

As opposed to the previous applications considered in the thesis, where the point cloud is known a priori (or generated randomly) and the interactions are computed using kernel functions, here the data are given as raw distance matrices, a.k.a. dissimilarity matrices. There are several ways of comparing DNA sequences (see [START_REF] Gusfield | Algorithms on strings, trees, and sequences[END_REF][START_REF] Yang | Computational Molecular Evolution. Oxford Series in Ecology and Evolution[END_REF]), the dissimilarity used here is computed from a local alignment score [61, p. 232] using Smith-Waterman algorithm [START_REF] Smith | Identification of common molecular subsequences[END_REF]. The entries of the distance matrices have a very low resolution, namely 0.05 for a maximum value below 650. In particular, the distance matrices can be multiplied by 100 and stored as short unsigned integers. Then, the similarity matrix S can be assembled very efficiently using single precision arithmetics and then stored in memory using single float. For n = 10.000 the similarity matrix is computed in 0.2s and stored using 400MBytes, therefore a full sample takes about 40GBytes in memory, which fits the 128GBytes available on plafrim2/miriel nodes and leaves room for storing a large number of components, namely up to r = n.

Best rank-r relative error: ∥S -S r ∥/∥S∥

Spectrum

In the following we will mostly focus on the properties and the low-rank structure of the similarity matrix S, that will often be refer to as a similarity or a covariance matrix, since it represents correlations between sequences of the sample. In order to design an efficient LRA technique, we first need to ensure that the similarity matrix fulfills the hypothesis of low numerical rank, i.e., r << n. The numerical rank is usually defined in term of the reconstruction error on the covariance and cannot be lower than the best rank-r error provided by the SVD. As shown on Figure 8.1, for subsamples up to 10.000 reads selected from a full sample at various times, the reconstruction error reach a value of 10 -3 in Spectral norm and about 10 -2 in Frobenius norm for a rank of about r = 200, i.e., 2% compression. However, one needs to set the rank to r ≈ n/2 in order to reach a relative error of 10 -3 in Frobenius norm. Moreover, the best rank-r relative errors of all subsamples (of any size and origin) almost coïncide for the first r = 100 eigenvalues. In particular, the similarity of all subsamples can be approximated with r = 10 within 10% of accuracy in Frobenius norm and r = 40 within 0.1% of accuracy in Spectral norm. These remarks on the numerical rank will be confirmed in Section 8.2 on larger subsamples and full samples. Finally, despite being rather flat, the spectrum of the similarity matrices associated to our samples decreases sufficiently fast (Figure 8.1 left) to be well approximated by random projection.

Distortion A key issue in our application is to ensure a certain accuracy when reproducing the point cloud, that we quantify in term of the distortion δ(X) , i.e., the representativity error of the computed components on the distance. First of all, the numerical rank should ideally be determined w.r.t. to the distortion, thus providing an application-oriented criterion. Second of all, the relation between the distortion and the accuracy of the low rank approximation of the similarity matrix should be investigated in order to set the target accuracy of our LRA algorithm w.r.t. to the target distortion. Figure 8.2 shows the evolution of the distortion w.r.t. to the rank r for the first 1.000 entries of a full sample. We can observe that the distortion decreases then reaches a plateau at about r = 10, where it stays slightly inferior to 10%, and then decreases again until it reaches its minimum value, namely about 0.1 to 0.01%. If 10% distortion seems sufficient then only 20 components can be used to represent the point cloud, for better accuracy one should consider much larger r, namely between 200 and 500. On the other hand, if we look at the fraction of the energy of the similarity matrix associated to negative eigenvalues, we observe that it grows fast once the first negative eigenvalue occurs. However, the relative contribution of the negative eigenvalues to the Frobenius norm, defined as

ε - r = ∑ r - i=1 (λ - i ) 2 ∑ r i=1 (λ i ) 2 , ( 8.19)
stops increasing after r = 200, where it reaches its maximum value, namely 5 to 6%. This illustrates the correlation between distortion and negative singular values mentioned earlier on.

Coherence A crucial information when performing subsampling is the coherence of the matrix, please refer to [START_REF] Talwalkar | Matrix coherence and the nystrom method[END_REF][START_REF] Mohri | Can matrix coherence be efficiently and accurately estimated?[END_REF] for nice introductions to the concept of matrix coherence and its relation to the performance of the Nyström method. If each column contains the same amount of information on the matrix, i.e., information is spread equally throughout the columns, then the matrix is incoherent. As a result, columns can be sampled uniformly at random in order to accurately represent the full matrix. On the other hand, if only certain columns concentrate valuable information on the full sample, then the matrix is coherent and columns should be selected more carefully. In order to quantify the importance of a column we usually associate the i-th column with a leverage score ℓ r i for i = 1, . . . , n defined as the Euclidean norm of the rows of the singular vectors of the input matrix, namely

ℓ r i = ∥U i ∥ 2 = √ U 2 i1 + . . . + U 2 ir . (8.20)
In particular, the leverage scores quantify the correlation between the columns and the basis formed by the first r eigenvectors. A common technique to select columns properly for coherent matrices is based on random sampling with probability proportional to the leverage scores [START_REF] Michael | Cur matrix decompositions for improved data analysis[END_REF][START_REF] Gittens | Revisiting the nystrom method for improved large-scale machine learning[END_REF], hence the higher the score of a column the higher the probability to sample that column. However, computing the leverage scores of a matrix requires computing the rank-r SVD of the input matrix, which is in general not affordable. A common alternative is approximating those scores by mean of a Randomized SVD (rSVD) [START_REF] Drineas | Fast approximation of matrix coherence and statistical leverage[END_REF], which make the design of an efficient rSVD even more relevant for MDS applications. On the other hand, there exist other methods to select relevant columns in a matrix that perform relatively good in practice, e.g., deterministic approach based on the pivoted QR Decomposition [START_REF] Engquist | A fast directional algorithm for high frequency acoustic scattering in two dimensions[END_REF], semi probabilistic approaches on k-cluster partitioning [START_REF] Aggarwal | Adaptive sampling for k-means clustering[END_REF] or combination of randomized and deterministic algorithms [START_REF] Boutsidis | An improved approximation algorithm for the column subset selection problem[END_REF]. The matrices studied here are expected to have low coherence since all individuals are supposed to have similar significance, as a matter of fact their presence in the sample is inherently random. This assumption will be verified in Section 8.2 since the computation of the coherence of the full sample can only be achieved by mean of a randomized SVD.

Visualization

The main goal of our approach is to achieve an efficient and handy visualization technique that provides relevant information on some real-life samples of about 100.000 sequences at various time of the year, namely about 10 measures corresponding to different months. First of all, we would like to be able to characterize the shape of the point clouds in the first few dimensions, identify potential clustering and relate clusters to identified species. Second of all, we want to compare these shapes between various times individually or using cross correlations. Finally, we would like to provide handy visualization technique to characterize clustering in larger dimensions.

Random projection aided MDS

In this section, we present our main contribution, namely a random projection based MDS algorithm, that allows for treating samples up to several hundred thousand reads on a single node in reasonable time, namely within a few minutes. The algorithm also allows for :

• computing distortion compared to the input distance,

• approximating the spectrum of the matrix,

• and evaluating the coherence of the matrix to be used for improving a subsequent subsampling.

We first discuss the benefits of random projection and present our approach. Then, we setup the parameters of the algorithm and we analyze the distortion on small subsamples. Finally, we analyze the performance of the algorithm on full samples.

Fast MDS based on random projection

Here we introduce a new random projection based MDS algorithm and show how it can contribute to the design of more efficient randomized MDS algorithm in general. Then, we perform simulations on subsamples with an intermediate size, namely 10, 20 or 30 thousand reads, in order to calibrate the input rank and the tune up parameters.

Random Projection

Our prior analysis of the data (Section 8.1.3) showed that similarity matrices involved in our application could be well-treated by random projection in term of both accuracy and performance. In particular, we highlighted the relatively low numerical rank of the similarity matrix S given the low target distortion, namely about 10%. Furthermore, within this range of accuracy, the spectrum associated to each sample appeared to decrease sufficiently fast and behave independently of the time t or subsample size. Therefore, we have implemented an MDS algorithm relying on a random projection LRA technique, namely the rSVD (Algorithm 6) introduced in Chapter 3. In this approach, the similarity matrix is approximated as follows S ≈ (QQ T )S(QQ T ), (8.21) namely the projection on a rank-r approximation Q of its range, itself obtained by dense Gaussian random projection (3.4). While being inherently related to dimensionality reduction, this approach presents many numerical benefits compared to random sampling technique, such as strong theoretical error bounds, near optimal accuracy in practice and good performance, while leaving significant room for improvement. Although the theory behind this approach is relatively involved, the method is easily implemented, stable and rather robust, which makes it both attractive and accessible to computational biologist for instance. Furthermore, the factorization in EVD form provided by the rSVD gives useful information on the structure of the similarity matrix, e.g., coherence, rank and any quantity that can be computed from the approximate eigen-pairs. In particular, we can estimate the reconstruction error in Spectral norm, compute the energy of the approximate similarity or even detect negative eigenvalues and evaluate the associated distortion.

Subsampling As discussed in Chapter 3, the design of an efficient randomized algorithm mostly relies on a good knowledge of the input data and matrices (structure, clustering,…) and it can often benefit from efficient combinations of both random sampling and random projection techniques. In our application, since working on the full sample or on a random subset should in essence lead to similar conclusions for a given sample, it can be very useful to perform a few computations on subsamples before performing computations on the full dataset. This allows us to draw preliminary conclusions on the structure of the similarity matrix in order to ultimately tune up our algorithm. Since our random projection algorithm can be used to efficiently approximate the coherence of the matrix [START_REF] Drineas | Fast approximation of matrix coherence and statistical leverage[END_REF], we can verify by analysis of the leverage scores whether or not subsets of the full sample can be selected uniformly at random. If it is not the case, then the algorithm provides useful information to select a more relevant subset. In that sense, our contribution not only provides an alternative to column selection algorithms with many extra benefits, but it could as well help design a more efficient column selection approach.

Fixed rank or accuracy Our new approach involves very little modification on the original MDS algorithm, as it only relies on a new way to perform the EVD of the similarity matrix S. In that sense, the target accuracy or the rank of the algorithm is set in the same fashion as a standard tSVD. In order to validate and calibrate our algorithm it is more convenient to prescribe the number of components that we want to visualize. Therefore, we prefer the fixed rank rSVD. Since we do not want to visualize more than a few hundred components, we will determine the prescribed rank to be the rank that leads to a few percent distortion on a length-n subset of the full sample such that r << n. The associated reconstruction error on the covariance matrix will ultimately be taken as input of the fixed accuracy rSVD.

Oversampling and Power iterations

As opposed to the tSVD, the rSVD provides approximate rank-r factorization, therefore it needs to be carefully tuned in order to make this approximation near optimal at a small extra cost. Figure 8.3 shows the accuracy of the rSVD on the spectrum of the similarity matrix for the sample L 6 . The rSVD accurately estimates the first 100 eigenvalues even without oversampling (s = 0), but it may fail to identify the last few eigenvalues. The figure also shows that significant oversampling (s = r) can substantially improve the accuracy. Although the accuracy still decreases slightly on the last part of the spectrum, the associated eigenvalues hardly contribute to the reconstruction error (Figure 8.1). For instance, for n = 10.000 random reads from L 6 , the reconstruction error achieved for r = 200 equals 2.1% in Frobenius norm while the best rank-r is 1.6%. Finally, simulations showed that this issue can be fixed by a single power iteration, nevertheless we prefer to avoid that and tolerate this small difference in order to develop a faster algorithm.

Validation of the method

Here we validate our method on various subsamples of a given L t sample, namely L 6 . We first show the benefits of using random projection to select a relevant subset of the full sample. Then, we discuss the accuracy of our approach on various subsamples and compare it with full samples.

Subsampling First of all, we want to verify whether or not uniform sampling of the columns applies well to our datasets. Let us analyze the leverage scores associated to the sample L 6 . Figure 8.4 shows some statistics such as the mean, variance and extrema of the leverage scores up to r = 1.000. In particular, we observe that for an artificial sample with n = 10.000 (Figure 8.4 left) the scores are spread equally around the mean and have a rather small variance. This shows that the similarity matrix associated with this sample is incoherent, i.e., all columns contain about the same amount of information on the sample. However, the full sample exhibit a different behavior (Figure 8.4 right) in the sense that some leverage scores are significantly larger than the rest and make the distribution deviates from Random reads from L 6 Le : First 500 eigenvalues of a random subsample of L 6 compared to approximate values obtained by rSVD without (s = 0) or with (s = r/2) oversampling. Right: First 200 approximate eigenvalues obtained by rSVD with significant oversampling s = r on larger subsamples: n = 10.000, n = 20.000 and n = 30.000. We use the result of a rSVD with r = 1.000 as an overkill to compare the approximate spectrum to a near op mal spectrum, but also to show that the rSVD delivers a result with very low variability. Observa ons: The approximate spectrum differs slightly from the exact one only for the smallest values, namely r > 100. Furthermore, as the size of the subsample n increases, the spectrum gets shi ed upwards but s ll has the same shape, i.e., same reconstruc on error, but the distor on may differ. its mean. On the other hand, the variance remains relatively small. Hence, uniform sampling may still be relevant for selecting columns in a Nyström based approach or simply selecting a relevant subset of the sample in a subsampled random projection approach. But a finer analysis of the leverage scores may lead to a more accurate representation at a lower cost.

Distortion In order to properly set the rank r of the rSVD we need to quantify the distortion on subsamples with intermediate number of reads, namely n = 10.000 to n = 30.000. Although r = 200 components lead to a reconstruction error of the similarity matrix of about 10 -3 in Spectral norm and 10 -2 in Frobenius norm (Figure 8.3), it only leads to a distortion of about 12% for 10.000 reads and 17% for 30.000 reads (Figure 8.5). Moreover, for n = 10.000 we need to use 5 as much components to reach 8% distortion, i.e., r = 1000. Hence, based on the distortion criterion the ranks of our matrices are not extremely low, namely r ≈ n/10.

Performance on real-life samples

Here we discuss the computational cost of the method in terms of memory and running time. We also discuss the time required to compute various statistics on the point cloud and the similarity matrix.

Memory footprint

Obviously, since the distance and similarity matrices are stored and applied in a dense way, the cost is quadratic in n. In particular, storing the distance matrix the full sample (right). The graphs represent the mean, variance, minimum and maximum values of the approximate leverage scores provided by rSVD w.r.t. the input rank r. Observa ons: For n = 10.000 the minimum and maximum values of the leverage scores are rela vely close to the mean (±20%), while the variance is very small, which shows that leverage scores are equally spread around the mean and rarely deviate from this value. Therefore, the matrix has low coherence and sampling with probability propor onal to the leverage scores boils down to sampling uniformly. On the other hand, the leverage scores associated with the full sample are distributed in a much more complex fashion, which would imply sampling with probability propor onal to them. 8% while the propor on of eigenvalues is always around 50%. On the other hand, for a given rank, the distor on grows with the problem size.

using short integers takes about 40GigaBytes for a full sample with 100.000 individuals and the similarity matrix 20GigaBytes. Since the computers used for our simulations, namely plafrim2/miriel, have 128GB RAM, the matrices associated with the computation of a full sample are relatively large but still fit in memory and leave significant room for the storage of the components X even for large r.

Running times Figure 8.6 shows various experimental running times measured during the computation of a subsample with n = 10.000 reads and a full sample of about 100.000 reads. First of all, reading the distance matrix and computing the similarity is rather cheap. In particular for the full sample, it requires as much time as computing the rSVD with r = 200, namely about 100s. Second of all, Figure 8.6 confirms that the cost of the rSVD grows linearly with the rank and quadratically with n. Third of all, as expected, the computation of the actual distortion of the point cloud representation has a similar complexity as the rSVD, namely O(rn 2 ). However, its magnitude is 10 times higher than the rSVD, which makes its computation intractable to large r, e.g., 6hours for r = 2000. Hence, a better understanding of the distortion is required here in order to develop estimators for this application before addressing faster factorization techniques. This work will be addressed in a forthcoming article, that introduces our random projection based MDS algorithm. Finally, the computation of the coherence is performed on the full sample by mean of a rSVD and only introduces a small O(nr) extra cost for computing the leverage scores. The plot shows the me required to read the distance and compute the similarity matrix S (do ed lines), the running me of the rSVD (solid lines) and the me required to compute the distor on δ (dashed lines).

Visualization of the point cloud

As shown in Section 8.2, enforcing a few percent distortion as a target threshold of our method can result in relatively high ranks. Although this can still be addressed efficiently using random projection, it is rather demanding in term of error computation. However, the analysis of the components of the point cloud is currently not done in a large number of dimensions. Until high dimension analysis is made possible, visualizing data in small dimensions, e.g., r ≤ 10, still provides numerous crucial information on the datasets. Here, we show representations of the full samples in 3D for easy visualization of the point clouds and basic comparison of their shapes. Then, we discuss the visualization of a cloud obtained by our method in higher dimensions.

Representation in 3D

The easiest way of visualizing a point cloud is to project the value of the components in 3D. This can be done by mean of three 2D projections on the first feature pairs, namely (F 1, F 2), (F 1, F 3) and (F 2, F 3).

Subsample

We first consider the cloud associated with a subsample, namely 10.000 random reads out of L 6 . Figure 8.7 represents the position of each individual on the first feature pairs. First of all, we observe that certain regions of the domain associated to a feature pair are more concentrated than others. Second of all, these regions can have various shapes and be well separated from the others, for instance in the plane (F 3, F 2). Another specificity of this cloud is that a small cluster lies very far from the rest of the cloud in the first dimension F 1 and is surrounded by all the other individuals in (F 2, F 3). Finally, such clustering of the data can be better shown by coloring the identified individuals w.r.t. their species. Hence, association between known species and features can be better understood. In particular, the presence of individuals from the species Nitzichia dissipita in the isolated cluster can probably help understand the meaning of feature F 1.

Full Sample

Let us now consider the cloud associated with the full sample L 6 computed by mean of a rank r = 50 rSVD. Figure 8.8 represents the position of each individual on the first feature pairs, while Figure 8.9 represents the concentration of the cloud. First of all, without both coloring and concentration, it would be harder to identify clusters than for the smaller subsample (Figure 8.7). Despite being 10 times as dense, the cloud associated with the full sample has a similar shape than that of the subsample. This can be observed more clearly on Figure 8.9, where the concentration of the population is represented in logarithmic scale. In particular, this figure shows that only very small regions of the full domain contain most of the individuals and are surrounded by regions of fast decreasing population. Second of all, the previous remarks on clustering of known individuals w.r.t. their species hold even stronger for the full sample. Moreover, a large part of the cloud contains unknown individuals but still exhibits a certain structure and a significant density, which suggests a deeper analysis of the data. In particular, a combined analysis of the concentration and the identified individuals seems crucial in order to better characterize regions of high concentration. Variations in time Let us recall that samples L 1 to L 10 represent monthly sampling over 10 months at the same location, namely Lake Geneva. Comparison between various samples can be discussed in 3D, where the point clouds exhibit significant differences along time but similarities as well. In particular, we still observe clustering and highly concentrated areas. However, the magnitude of the features, i.e., the size of the domain containing the cloud, may vary by a factor of 2 to 3 in certain dimensions. For instance, L 1 and L 2 spread over a similar domain, namely a box of size roughly equal to 60 and centered in (0, 0, 0), while L 3 and L 4 spread over a smaller box that is slightly shifted from the origin. Moreover, rotation or translation of the point clouds should be considered in the analysis of the shape as the result of MDS is only defined up to a combination of such transformations. Observa on: The shape of the cloud varies significantly along the me, while s ll presen ng similar characteris cs such as clear clustering and highly concentrated areas.

L 1 L 2 L 3 L 4

Representation in many dimensions

Parallel Coordinates In order to visualize the point cloud in many dimensions a common alternative to 2D plots is the representation in parallel coordinates such as the one displayed on Figure 8.11. This technique offers various advantages as it decouples the features and displays them on a single axis. As shown on Figure 8.11 (top left), the observation of all reads on the same plot without characterization of the concentration cannot be exploited, which suggests the use of a different plotting software with better rendering than matplotlib or a finer pretreatment of the components. However, the expected clustering of the reads can be confirmed by observing that individuals of the same species (top right and bottom left), follow a similar path along the first few dimensions, while spreading only on a fraction of the full domain. This observation is made more clear by the observation of each species separately (bottom left). As the index of the feature grows the associated eigenvalue of the similarity matrix decreases, however it has almost no impact on the first 10 features and can only be noticed on more than 100 features.

Conclusion

Conclusion

Conclusion

Let us summarize our various achievements and contributions one topic after another.

Fast Multipole Method

We developed a highly competitive interpolation based Fast Multipole Method (FMM), called the ufmm, relying on a very simple idea and benefiting from the performance of the Fast Fourier Transform (FFT). It provides a kernel independent M2L-optimized approach that may suffer from instabilities at very high accuracies but performs significantly better in all cases of interest than the state-of-the-art approaches, namely the bbfmm and optimized variants. This method was implemented in the open-source parallel HPC library ScalFMM and it is already used by several users for diverse applications such as Molecular Dynamics or Boundary element methods for elastostatics. In this thesis, we used the ufmm in various computationally intensive applications with rather low target accuracies in order to demonstrate its good performance and illustrate its genericity.

Dislocation Dynamics

We implemented the isotropic elastic force and energy farfield computation based on the ufmm and bbfmm within OptiDis, an operational parallel Dislocation Dynamics (DD) code developed in collaboration with the CEA Saclay. In particular, it already allowed for simulating the phenomenon of clear band channels in Zirconium alloys, as presented in Arnaud Etcheverry's thesis [START_REF] Etcheverry | Simulation de la dynamique des dislocations à très grande échelle[END_REF]. An efficient summation scheme was proposed along with several variants minimizing the cost of certain operators. We adapted certain features of the code to accurately and efficiently handle segments within an octree structure, exact integration of the farfield as well as the specific nature of the interaction kernel (tensorial, homogeneous, singular,…). In particular, this lead us to develop a generic formalism for tensorial interactions within ScalFMM library. Additionally, this allowed us to illustrate the rather high memory requirements of the ufmm, compared to the bbfmm, while proposing alternative formulations that could compensate this phenomenon. Moreover, this issue arise for high tree depths and only affects the storage of the expansions that is usually spread over several nodes in a shared memory parallelism.

Sampling from GRF A radically different application was addressed, that consisted in accelerating the matrix multiplication involved in a Randomized Low-Rank Approximation (LRA) algorithm by mean of the ufmm and thus provide low-rank square root of covariance matrices. This allowed us to extend the capabilities of our algorithm to matrix-to-matrix multiplication using a multi-rhs FMM implementation. The matrix studied here were covariance kernel matrices used for sampling from Gaussian Random Field (GRF), i.e., generating multivariate Gaussian Random Variables. Several correlation kernels were tested, such as Gaussian, Exponential and several other Matérn correlations, in order to generate GRFs on artificial heterogeneous distributions of points. We designed a smooth variant of the ufmm in order to represent smooth kernels more efficiently. This was tested on the Gaussian kernel but this should have ideally been tested on a less nice kernel, for which there does not al-ready exist fast algorithms. To our knowledge, this contribution is one of the rare attempts at combining randomized LRA with the data sparse representation of a matrix. The results were promising and setting the FMM accuracy w.r.t. to the global low-rank approximation turned out to be extremely easy.

Random projection aided MDS

We addressed another shortcoming in randomized numerical linear algebra, namely the lack of fast Multidimensional Scaling (MDS) algorithm based on random projection. Application of a Randomized Eigen Value Decomposition (EVD) on the similarity matrix allowed for deriving an efficient MDS algorithm that provides useful information on the structure of the matrix. This method is complementary with the state-of-the-art method, namely Nyström factorization, while being applicable to any similarity matrix. We provided a discussion and several experiment to relate the distortion with the presence of negative eigenvalue in the spectrum and a best or near-best low-rank representation. The method was applied on biological datasets coming from real-life samples and an operational workflow was organized during the thesis in order to provide visualization data and thus help design a geometric view on biodiversity. The algorithmic and numerical aspects of our random projection based approach to MDS is the topic of a forthcoming paper, that will also present the idea behind this geometric view on biodiversity, with strong references on the theory of dimensionality reduction.

Perspectives

A significant amount of work is still ongoing in order to consolidate these contributions, however our current point of view allows us to draw a considerable number of perspectives.

Uniform FMM The main drawback of the ufmm should be dealt with in a near future, namely the relatively high amount of memory required for the storage of the expansions in Fourier domain. A first idea would be to transfer expansions to Fourier domain only when it is actually required. Then, an interesting problem would be the application of the ufmm in a directional Fast Multipole Boundary Element Formulation [START_REF] Traub | A kernel interpolation based fast multipole method for elastodynamic problems[END_REF][START_REF] Messner | Fast Boundary Element Method in Acoustics[END_REF] where the cost of the M2L can be critical but the interaction lists can be very large as well. Another issue that we need to address would be the regularization of the equispaced interpolation involved in the ufmm, as instability are very likely to occur for oscillatory kernels for instance. Additionally, the significant improvement of the M2L step should allow the implementation of much complex FMM schemes such as a variable order FMM [START_REF] Tausch | The variable order fast multipole method for boundary integral equations of the second kind[END_REF] or H 2 -methods [START_REF] Börm | Approximation of integral operators by variable-order interpolation[END_REF]. Most importantly, the ufmm and its smooth variant, the smooth-ufmm, are in essence so close from the general H-matrix formats, therefore we hope that our modest contribution can motivate the design of a unified framework for hierarchical low-rank matrix approximations behind the H 2 format that can benefit to both the FMM and the H-matrix communities.

In our opinion, a first step would be for us to compare the performance of our method with recent interpolation based FMM [START_REF] Casenave | An empirical interpolation based fast summation method for translation invariant kernels[END_REF] and commonly used interpolation based H 2 -methods [START_REF] Hackbusch | H2-matrix approximation of integral operators by interpolation[END_REF][START_REF] Dölz | An interpolation-based fast multipole method for higher-order boundary elements on parametric surfaces[END_REF].

Anisotropic DD First of all, our fast approach to isotropic DD simulations should allow to perform the accurate simulation of massive ensembles of dislocations and thus better understand the plastic behavior of complex crystalline materials. Although our contribution to the isotropic force evaluation was in our opinion not novel enough, we think that writing an article on an interpolation based FMM formulation for both isotropic and anisotropic elastic fields would represent a significant and novel contribution to the domain of Multiscale Materials Modeling. In fact, anisotropic models should allow to better simulate and understand the behavior of materials that exhibit large anisotropic coefficients like Fe. During this thesis, an efficient representation of the anisotropic force field based on the Stroh formalism and expansions in spherical harmonics were also addressed, that should allow for an efficient fast multipole implementation in a similar fashion as the recent methodology presented in [START_REF] Aubry | Use of spherical harmonics for dislocation dynamics in anisotropic elastic media[END_REF]. Special care should be taken in order to derive an equivalent non-singular formalism for anisotropy as well as analytic expressions for the integrals of the expansion in spherical harmonics. The ufmm should efficiently address the high dimensional interactions involved in anisotropic DD simulations and make the comparison of various formulations easier.

Applications to Data Assimilation

We hope that our H 2 -powered randomized Singular Value Decomposition (SVD) algorithm can help solving large scales problems of geostatistics and that it can be beneficial to the computational geosciences community. Therefore, we are considering applying the method to a more active topic of geostatistics, namely data assimilation via Kalman Filtering. Recent advances in Kalman Filtering have brought ensemble variants and thus low-rank square root algorithms in front of the scene, namely Ensemble- [START_REF] Evensen | The ensemble kalman filter: Theoretical formulation and practical implementation[END_REF], Unscented-Kalman Filters [START_REF] Simon | Unscented filtering and nonlinear estimation[END_REF] and reduced order variants [START_REF] Padilla | An adaptive-covariance-rank algorithm for the unscented kalman filter[END_REF][START_REF] Chandrasekar | Reduced-rank unscented kalman filtering using cholesky-based decomposition[END_REF]. As these methods rely either on efficient computation of a covariance matrix square root or generation of correlated noises, they provide an ideal framework for the application of our algorithm. However, designing a relevant application to emphasize the benefits of our approach is a rather challenging issue. During this thesis we met or exchanged with people from several teams of the data assimilation community, namely hydrologists, climatologists and oceanographists from either CERFACS or Stanford Civil Engineering Department, but only recently came up with an application that could benefit to a broader community.

Linear scaling MDS We are planning on developing a linear scaling factorization technique in a near future in order to address very large biological datasets in a more efficient way. This approach should benefit from the clustering phenomenon occurring on the point clouds and will benefit from recent advances in k-means algorithms and graph partitioning. Moreover, we would like to be able to better characterize and compare the shape of the clouds between various time of the year. In particular, we wish to analyze the point cloud associated to several samples at the same time. Another crucial perspective consists in designing new metrics as well as efficient application-oriented visualization tools in the spirit of Chris Johnson's work, e.g., [START_REF] Rhyne | Information and scientific visualization: Separate but equal or happy together at last[END_REF].

Concluding remarks Randomized Numerical Linear Algebra (NLA) is a relatively young research area compared to the FMM or H-matrices and it is most definitely much wider. However, recently the combination of both randomized and hierarchical algorithms seems to be a topic of broad interest [START_REF] Wang | Structured block basis factorization for scalable kernel matrix evaluation[END_REF][START_REF] Victor Minden | Fast spatial gaussian process maximum likelihood estimation via skeletonization factorizations[END_REF]. A crucial issue here is to develop new hierarchical methods with lower dependence in the ambient dimension, a topic that we wish to address very soon. We hope that our open-source library for the randomized low rank approximation of matrices can benefit from the variety of topics addressed in this thesis and from our upcoming publications. Finally, we hope that the hierarchical and randomized NLA communities can keep growing and that we can continue exchanging on these beautiful topics, their complementarity as well as their numerous applications.

B

Distribution of particles on various geometries Here we give advanced statistics on various distributions of points w.r.t. to the geometry and the depth of the octree.

B.1 Geometries with various heterogeneity

B.2 Statistics on the octree

The evolution of number of non-empty cells w.r.t. the number of level for n = 10 6 particles is represented on Figure B.2. In particular, we observe that the unit cube is fully populated up to L = 6 then the filling rate decreases with the tree depth. Meanwhile filling rates of the unit and prolate spheres have similar decreasing speed, but as expected the prolate sphere has more empty leaves than the unit sphere. Figure B.3 represents the total and average number of nearfield and farfield interactions w.r.t. to the tree depth L. When the particles are distributed in the unit cube, the number of points is not high enough to have a sufficient density and maintain the filling rate. a surface, e.g., the unit sphere or the prolate sphere, the tree is very sparse (low filling rate). Indeed, the number of non-empty leaves differs significantly from the maximum number of leaves, i.e., 8 L. Right: Average number of par cles per leaf w.r.t. to L for various geometries with n = 10 6 par cles. In prac ce, leaves contain between 10 and 500 par cles in average (green area) depending on the applica on, i.e., on the nature of the kernel, the accuracy and the expansion technique. L. Observa ons: The average number of operators stays rela vely constant for all geometries, in par cular we quickly reach the maximum number of operators for the unit cube, namely 27 P2P and 189 M2L, while the unit sphere only has 15 P2P and 50 M2L in average. The prolate sphere exhibit 2 regimes, namely L ≤ 4 or L ≥ 5, whether farfield interac ons are considered between par cles lying on the small sec on or not.

C.1 Numerical benchmarks

Here we present comparative results on the convergence of the multi-rhs ufmm and smoothufmm w.r.t. the interpolation order p for various geometries: the unit sphere, a cube and a prolate sphere. In particular we analyze the effect of the length scale on the error of the Fast Multipole Method (FMM) for the Gaussian correlation kernel for ℓ = 0.5 or ℓ = 1.0. All computations were performed on a cluster computer, namely plafrim/mirabelle: Hexa-core Westmere Intel Xeon X5670 @ 2. C.2 Performance of the ufmm and the smooth-ufmm.

The ufmm and the smooth-ufmm exhibit similar performances if the tree depth is set to the same value. However, as the depth decreases, the smooth-ufmm becomes significantly cheaper. Hence, the results presented here corroborate the one discussed in Section 5.4.

Precomputing the P2P operators takes about 3 to 8 seconds depending on the geometry and the tree depth, which is approximately the time required to apply the ufmm to 1 rhs with an intermediate accuracy. Observa ons: ufmm and smooth-ufmm have approximately the same cost when they share the same tree depth. The smoothufmm is significantly faster if we choose an op mal tree depth, e.g., L = 3 represented in solid brown lines. If the Gaussian decreases sufficiently slow, i.e., the rank is sufficiently low, then the cost of the global is similar to the op mal smooth-ufmm. However, the global exhibits an instability for the highest interpola on order.

Prolate Sphere

D.1 M2L-Optimized stress computation

In this section we start by introducing an existing approach used to optimize the number of M2L operators and operations per interaction in the case of isotropic elastic forces computation. Then, we present the principle of our new method based on refactoring the M2L operations. Finally, we apply our method on the evaluation of the isotropic elastic energy.

D.1.1 Existing approach

Arsenlis approach In [8] the isotropic elastic stress field (2.4) is written as a combination of permutation symbols and a 5 th -order tensor G ijklm defined as

G ijklm (x) = C R ,klm (x -y)b i dy j (D.1)
Therefore, the contributions to the elastic stress read as

σ A ij (x) = ε jmk G kimpp (x) (D.2) σ B ij (x) = ε nmk G knmij (x) (D.3)
Let us determine the number of independent components of G that are actually involved in the computation of σ. First of all, given the symmetries of R ,klm only the following 10 components (instead of 3 3 = 27) of G ij... need to be computed for all (i, j) ∈ {1, . . . , 3} 2

G ij111 G ij122 G ij133 G ij123 G ij211 G ij222 G ij233 G ij311 G ij322 G ij333
Besides only the following components of G ...lm for all (l, m) ∈ {1, . . . , 3} 2 are actually required when the permutation symbols are applied:

σ A 11 = -G 213pp + G 312pp , σ A 12 = -G 311pp + G 113pp , σ A 13 = + G 211pp -G 112pp , σ A 21 = -G 223pp + G 322pp , σ A 22 = -G 321pp + G 123pp , σ A 23 = + G 221pp -G 122pp , σ A 31 = -G 233pp + G 332pp , σ A 32 = -G 331pp + G 133pp , σ A 33 = + G 231pp -G 132pp , σ B lm = +G 231lm -G 321lm -G 132lm + G 312lm + G 123lm -G 213lm ,
The other components are multiplied by 0. If we combine these 2 optimizations we can reduce the set of computed components from • The computation of σ A ij for all (i, j) ∈ {1, . . . , 3} 2 involves the following (3×6)×3 = 54

D.2 M2L-Optimized energy computation

Here, we present our implementation of the farfield contributions to the isotropic elastic energy, whose expression is recalled in Section 2.2.2. In order to remain consistent with formula already implemented in OptiDis, we used the formula from [START_REF] Cai | A non-singular continuum theory of dislocations[END_REF]. However, the efficient fast multipole implementation is more challenging than that of the stress field, therefore we propose to use a different formula in order to reduce the cost of computing the energy using Fast Multipole Method (FMM).

D.2.1 Fast Multipole computation of the isotropic energy

Interpolation If we replace the expression (6. These 4 contributions can also be expressed as the sum of the contributions of each segment composing the network.

Fast Multipole summation scheme Equation (D.4) can be computed using a fast multipole summation scheme, that consists in the following steps:

• First, aggregate segments contributions at interpolation nodes

(M β ) ij = ∫ [x 1 x 2 ] S β (x ′ )b ′ i t ′ j dx ′ (P2M-E)
• Then, transfer multipole expansion to local expansion for

E C or E D L A α = ∑ |β|≤p+1 ( Kαβ ) kk (M β ) ii (M2L-EA) (L B α ) jk = ∑ |β|≤p+1 ( Kαβ ) ij (M β ) ik (M2L-EB) (L C,D α ) ij = ∑ |β|≤p+1 ( Kαβ ) kk (M β ) ij (M2L-ECD)
• Finally, accumulate contributions at target points

E A ( [ x 1 x 2 ] , [ x 3 x 4 ] ) = ∑ |α|≤p+1 ∫ [x 3 x 4 ]
S α (x)b j t j L A α dx (L2P-EA)

E B ( [ x 1 x 2 ] , [ x 3 x 4 ] ) = ∑ |α|≤p+1 ∫ [x 3 x 4 ]
S α (x)b j t k (L B α ) jk dx (L2P-EB)

E C ( [ x 1 x 2 ] , [ x 3 x 4 ] ) = ∑ |α|≤p+1 ∫ [x 3 x 4 ]
S α (x)b i t j (L C,D α ) ij dx (L2P-EC)

E D ( [ x 1 x 2 ] , [ x 3 x 4 ] ) = ∑ |α|≤p+1 ∫ [x 3 x 4 ]

Cost estimation

In order to perform this summation, we need to store 7 M2L components, namely 6 for ẼB and 1 for the other terms. Moreover, we need to store 9 components of M β . On the other hand, one need to apply 47 M2L operations in order to compute all terms of the energy, namely 1 for ẼA , 27 for ẼB and 9 for ẼC,D .

D.2.2 Efficient implementation

Direct Approach In order to evaluate the isotropic elastic energy via FMM we need to compute 9 multipole expansions at the P2M step, namely

M ij = (C)
S(x ′ )b i dx ′ j , for i, j = 0 . . . 2

Then, apply R ,kk to each of these expansions, i.e., 9 M2L operations, and store the results in

L B 0 = R 0 M 0 + R 1 M 3 + R 2 M 6 , L B 3 = R 1 M 0 + R 3 M 3 + R 4 M 6 , L B 6 = R 2 M 0 + R 4 M 3 + R 5 M 6 , L B 1 = R 0 M 1 + R 1 M 4 + R 2 M 7 , L B 4 = R 1 M 1 + R 3 M 4 + R 4 M 7 , L B 7 = R 2 M 1 + R 4 M 4 + R 5 M 7 , L B 2 = R 0 M 2 + R 1 M 5 + R 2 M 8 , L B 5 = R 1 M 2 + R 3 M 5 + R 4 M 8 , L B 8 = R 2 M 2 + R 4 M 5 + R 5 M 8 ,
In conclusion, for the approach based on Cai's formula, the local expansions are stored in a vector of size 18 and involve 36 M2L operations with 9 different multipole expansions:

L ′ I =
M ij = S(x ′ )b ′ j dx ′ i L AB nj = R kk M nj and L C i = R ij ε jmn M nm
In this variant, 2 different local expansions have to be stored and contribute to the 3 terms of the energy. The first one is equal to the L A in Cai's variant and can thus be written as

L AB I = R 6 M I
Let us denote M ′ i = -ε imn M mn then the second kind of local expansion reads as

L C 0 = R 00 M ′ 0 + R 01 M ′ 1 + R 02 M ′ 2 L C 1 = R 10 M ′ 0 + R 11 M ′ 1 + R 12 M ′ 2 L C 2 = R 20 M ′ 0 + R 21 M ′ 1 + R 22 M ′ 2
In conclusion, for the approach based on Mura's formula, the local expansions are stored in a vector of size 12 and involve 18 M2L operations with 12 different multipole expansions:

L ′ 0 = R 0 (M 7 -M 5 ) + R 1 (M 2 -M 6 ) + R 2 (M 3 -M 1 ) L ′ 1 = R 1 (M 7 -M 5 ) + R 3 (M 2 -M 6 ) + R 4 (M 3 -M 1 ) L ′ 2 = R 2 (M 7 -M 5 ) + R 4 (M 2 -M 6 ) + R 5 (M 3 -M 1 ) L ′ 3+I = R 6 M I , I = 0 . . . 8

D.3 Numerical/analytical integration over segments

In this subsection we discuss the evaluation of the integral over the segments. The elastic stress field (2.3), respectively the force field (2.7) and energy (2.8) involves 1, respectively 2 successive, integrals over a line in R 3 . The evaluation of the integrand is done via polynomial interpolation, hence the numerical integration can be achieved exactly using a Gaussian quadrature with a well adjusted size. Let us recall the integral of an interpolated function in 1D over a subinterval S p (Φ -1 (x), xm )dx

A usual method to compute such integrals numerically is to use a Gaussian quadrature. This requires to map the interval of integration to [-1, 1] first Once this is done the integral is computed by summing the results of evaluating the integrand at Gauss point X q weighted by corresponding Gauss weights w q such that ∫

I p = p ∑ m=0 f (Φ(x m )) ∫ [-1,1] S p (Φ -1 ( x 2 -x 1 2 X + x 1 + x 2 2 ), xm ) x 1 + x 2 2 dX = p ∑ m=0 f (Φ(x m )) ∫ [-1,1] S p (Φ -1 ( L 2 X + c S ),
[-1,1] g(X)dX ≈ Q ∑ q=1 w q g(X q ) (D.7)
where Q denotes the size of the quadrature. Gaussian quadrature of size Q can achieve exact integration for polynomial integrands of order at most 2Q + 1 which make them very competitive.

Index, Glossary and Bibliography

Acronyms

  .p par un polynôme S d'ordre p. Cette technique revient à approximer le noyau d'interaction par la formule suivante k(x, y) ≈ ∑ a≤p S(x, xa ) ∑ b≤p k(x a , ȳb )S(ȳ b , y)

Figure 1 :

 1 Figure 1: Appari ons de bandes claires dans les alliages de Zirconium (Zr) simulées avec Op Dis.

  to be in farfield interactions with C (L) x , while the other are in nearfield interactions. The neighbor list N (C (L) x ) is defined as the set of cells in nearfield interactions with C (L) x . For γ = 2, N (C (L) x ) contains at most 3 d cells for any L, i.e., 27 in 3D. The interaction list I(C (L)x ) is defined as the set of non-empty cells C

Figure 1 . 1 :

 11 Figure 1.1: Por on of a 1D cluster tree at level L and L -1. The interac on lists I(C (L) x ) (green) and neighbor lists N (C (L) x ) (red) of a target par cle x is represented for γ = 2. The contribu on of the source par cle y a is computed analy cally as part of the nearfied while contribu ons of y b and y c are approximated as part of the farfield at level L and L -1 respec vely.

Figure 1 . 2 :

 12 Figure 1.2: The por on of a quadtree at level L and L -1 used to precompute the M2L operators. The interac on list I(C (L) x ) (green) and neighbor list N (C (L) x ) (red) of a target par cle x are represented at the leaf level for γ = 2.

  j , ȳβ )w(y j ), ∀β/|β| ≤ p else // M2M:

2

 2

Figure 2 . 1 :

 21 Figure 2.1: Illustra ons of a material along the various steps of the mul scale modelling. The axis show the typical me and length scales associated to the methods.

Figure 2 . 2 :

 22 Figure 2.2:A disloca on loop discre zed in 6 segments. All segments share the same Burger's vector b as the loop but they have various line tangents t depending on their orienta ons.

Figure 2 . 3 :

 23 Figure 2.3: Parametriza on of 2 interac ng disloca on segments, the source [ x 1 x 2 ]

Contents 3 . 1

 31 The Randomized SVD . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.1 Introduction to randomized algorithms . . . . . . . . . . . . . . . 3.1.2 The Randomized SVD . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Fast Randomized LRA . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.1 Fast random projection . . . . . . . . . . . . . . . . . . . . . . . 3.2.2 Fast random sampling . . . . . . . . . . . . . . . . . . . . . . . .

  denotes the first r left singular vectors, resp. right singular vectors, associated with the first r singular values σ(C). According to Eckart-Young theorem, C r provides the best rank-r approximation of C. In particular, if the rank of C equals r, then C = Cr . For self-adjoint positive definite matrices C ∈ R n×n , the SVD boils down to C = UΣU T , (3.3) a.k.a. the Eigen Value Decomposition (EVD). If C has negative eigenvalues, then Σ contains their absolute value. The full EVD can be obtained at an O(n 3 ) cost using optimized LA-Pack routines (e.g., Intel MKL wrapper), while the rank-r truncated EVD can be obtained iteratively at an O(rn 2 ) cost using the Arnoldi's algorithm provided by ARPack.

Contents 4 . 1

 41 Efficient interpolation based FMM . . . . . . . . . . . . . . . . . . . . 4.2 Fast Multipole DD simulations . . . . . . . . . . . . . . . . . . . . . . 4.3 Fast Randomized LRA . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 5 . 1 :

 51 Figure 5.1: Circulant embedding (right) of a symmetric 1D Toeplitz matrix (le ), namely a Gaussian kernel matrix evaluated on a 1D grid of 16 points. The Toeplitz matrix (blue frame) is embedded in the upper le corner of the circulant matrix (right).

Figure 5 . 2 :

 52 Figure 5.2: Circulant embedding (right) of a symmetric 2D Toeplitz matrix (le ), namely a Gaussian kernel matrix evaluated on a 1D grid of 4-by-4 points. The small Toeplitz matrices (blue frames) are embedded in the upper le corner of the small circulant matrices (right), while blocks are also embedded in the upper le corner (red frame) of a circulant block structure (right).

  The discrete convolution theorem implies that the set of eigenvectors of any circulant matrix Ē coïncides with the Discrete Fourier Transform (DFT) operator, i.e., F = { p-1/2 e -2iπmn/p } m,n=0,...,p , while the DFT of the first column of Ē yields the vector Λ containing the eigenvalues of Ē, i.e., Λ = F R. Let us consider a multipole expansion M ∈ R p+1 and a resulting local expansion L ∈ R p+1 , such that L = KM(5.4)

Figure 5 . 4 :

 54 Figure 5.4: Le : Padding and Forward DFT of the mul pole expansion. Center: Applica on of the M2L operator in the Fourier domain. Right: Backward DFT and unpadding of the local expansion.

Figure 5 . 5 :

 55 Figure 5.5: Illustra on of both ufmm variants on a 1D tree. The target par cle x (blue) lies in the leaf cell C L x . The interac on

Figure 5 . 6 :

 56 Figure 5.6: Theore cal memory requirements of the interpola on based FMM in the worst case scenario, namely no empty cell.These values should be qualified by the fact that for heterogeneous geometries the number of non-empty cells can be very low.

Figure 5 . 7 :

 57 Figure 5.7: Uniform distribu on of par cles in the unit cube (le ), on the unit sphere (center) and on a prolate sphere (right) with propor ons (0.1, 0.1, 1.0).

Figure 5 .

 5 Figure 5.7 shows examples of such particle distributions for n = 10 5 particles. All theses geometries share the same bounding box, i.e., the 2 × 2 × 2 box, however they fill the octree differently. More precisely, the unitCube has a fully populated octree, while the unitSphere and the prolateSphere have respectively an intermediate and a low number of non-empty cells. For advanced statistics on the filling of the octree and on the number of nearfield and farfield operators please refer to Appendix B.

Figure 5 . 8 :

 58 Figure 5.8: Shape of the kernel w.r.t. the distance r if the width of the root bounding box equals 2. The limit between nearfield and farfield is represented for L = 2. Le : Nega ve powers of r and some regularized variants. Right: Exponen al and Gaussian kernels represent the extreme cases of smoothness in the family of Matérn func ons.

Figure 5 . 9 :

 59 Figure 5.9: Experimental accuracy of the ufmm and variants of the bbfmm w.r.t. the interpola on order p. We considered a uniform distribu on of 20.000 par cles inside the unit cube with L = 3. Le : Homogeneous Laplacian kernel k(r) = 1/r. Right: Non-homogeneous Gaussian kernel k ∞ (r, ℓ) with ℓ ∈ {0.5, 1.0}. The computa on is done on simple desktop computer.

Figure 5 . 10 :

 510 Figure 5.10: Experimental memory requirements and computa onal me of the ufmm and variants of the bbfmm w.r.t. the interpola on order p. We use the homogeneous Laplacian kernel k(r) = 1/r on a uniform distribu on of 20.000 par cles inside the unit cube with L = 3. The computa on is done on simple desktop computer.

TotalFigure 5 . 11 :

 511 Figure 5.11: Detailled computa onal me for the symmetric-bbfmm (le ) or the ufmm (right) w.r.t. the interpola on order p. We used n = 20.000 par cles distributed uniformly in the unit cube with L = 3.

Figure 5 . 12 :

 512 Figure5.12: Experimental memory requirements and computa onal me of the ufmm and smooth-ufmm w.r.t. the interpola on order p. We use a smooth and non-homogeneous Gaussian kernel k ∞ (r, ℓ) with ℓ = 0.5 and a uniform distribu on of 20.000 par cles inside the unit cube with L = 3.

2r 2 a -r 2 r 3 aFigure 6 . 1 :

 361 Figure 6.1: Representa on of the components k 11 , k 12 and k 11 + k 22 in 2D, i.e., r =√ r 2 1 + r 2 2 , for r 1 = r 2 ∈ [0, 2].Values are plo ed w.r.t. the scaled distance x = r/ω min for a = 3 and a typical segment length L max = 90, i.e., a = L max /30. We consider the cri cal case, where the leaf cell width equals ω min = 3L max .

Figure 6 . 2 :

 62 Figure 6.2: Theore cal memory requirements at the leaf level w.r.t. the tree depth in the worst case scenario (no empty leaf)

Figure 6 . 3 :

 63 Figure 6.3: Illustra on of the bounding box extension in the worst case scenario, where the center of a segment [ x 1 x 2 ] lies on the edge of a leaf cell.

9 Figure 6 . 4 :

 964 Figure 6.4: Comparison of the homogeneous kernel k(r) = 1/r with the non-homogeneous kernel k a (r) = 1/ √ r 2 + a 2 for 2 representa ve values of a. The kernel values are plo ed w.r.t. the scaled distance r/ω min with ω min = 3L max using a typical maximum segment length of L max = 90. Non-homogeneous kernels are represented for a = L max /30 = 3 (average core width) and a = L max /10 = 9 (large core width). Observa ons: The non-singular kernels coïncide with the singular

(a) 4 × 4 × 3 Figure 6 . 5 :

 44365 Figure 6.5: 3 configura ons of BCC loops with b = [111]. The number of segments is driven by the number of loops, i.e., the number of leafs, while the number of segments per loop remains constant. The exact nodal force modules are represented on the disloca on lines for each configura on.

(a) d = 5 .

 5 10 19 m -2 , n = 11.250, L = 4 (b) d = 10 20 m -2 , n = 90.090, L = 5

Figure 6 . 6 :

 66 Figure 6.6: 3 configura ons of BCC loops with b = [111]. The number of segments is driven by the density d, while the segment length L = 20 and loop diameter D = 100 remain constant. The exact nodal force modules are represented on the disloca on lines for each configura on.

  order: p relative L2 error:

Figure 6 . 7 :

 67 Figure 6.7: Accuracy of the nodal forces (le ) and energy (right) evalua on for different interpola on schemes (ufmm and bbfmm, with and without shi ed deriva ve) w.r.t. the interpola on order p. The farfield core width is set to 0. Configura on (b) of Fig. 6.5.

Figure 6 . 8 :

 68 Figure 6.8: Computa onal me required by the nodal forces evalua on for different interpola on schemes w.r.t. the number of segments n. The tree depth L is displayed under the marks of the cutoff variant. Le : Configura ons of Fig. 6.5, ensures aconstant number of segments per leaf and a fully populated octree. Right: Configura ons of Fig.6.6 is more realis c, however the asympto c complexity is harder to emphasizes, since there exists empty leaf cells. Observa ons: Only the ufmm-Facto and the bbfmm-Symm variants have well-balanced nearfield and farfield. For a given n, the bbfmm-Facto is 3 mes slower and the ufmm-Shift is about 2 mes faster. Machine: plafrim2/miriel 2 Dodeca-core Haswell Intel Xeon E5-2680 @ 2,5 GHz RAM 128Go (DDR4 2133MHz).

  63 4.64 • 10 -2 0.74 1.26 • 10 -4 1 • 10 4 8.89 • 10 -3 0.31 1.84 • 10 -2 0.4 5.07 • 10 -4 1 • 10 5 1.76 • 10 -3 0.7 5.18 • 10 -3 0.55 5.79 • 10 -3

Figure 7 . 1 :Figure 7 . 2 :

 7172 Figure 7.1:Convergence of a univariate GRV in terms of sample mean error E(X) -0 and sample variance error E(X 2 ) -1 w.r.t. the number of realiza ons n real . We used a Mersenne Twister random number generator with mean 0 and variance 1. The rates of convergence (roc) are defined as the ra o of the logarithm of two consecu ve values of an error. Finally, t X denotes the me required to generate the random variable and compute the errors. We observe the theore cal 1/2 convergence rate for the mean and variance.

Figure 7 . 4 :

 74 Figure 7.4: Accuracy of the fixed rank rSVD w.r.t. the compression rate for a Gaussian correla on and n = 2.000 par cles distributed on the unit sphere. Le : Op mal rela ve error in Spectral and Frobenius norms w.r.t. to r = 1, . . . , n. Right: Accuracy of the rSVD w.r.t. prescribed rank r = 20, . . . , 200.Effect of the oversampling s alone (q = 0). Observa ons: The Gaussian kernel is smooth, therefore its spectrum decreases fast and C is rela vely low-rank. Consequently, the rSVD performs well even without power itera ons and with low oversampling (s = 5). An oversampling of s = 50 leads to a near op mal error. Note that C becomes high-rank for small ℓ, e.g., ℓ = 0.1.

Figure 7 . 5 :

 75 Figure 7.5: Accuracy of the fixed rank rSVD w.r.t. the compression rate for an exponen al correla on and n = 2.000 par cles distributed on the unit sphere. Le : Op mal rela ve error in Spectral and Frobenius norms w.r.t. to r = 1, . . . , n. Right: Accuracy of the rSVD w.r.t. prescribed rank r = 20, . . . , 200.Effect of q power itera ons using a low oversampling, namely s = 5. Observa ons: The exponen al correla on k 1/2 is not smooth, therefore its spectrum is rela vely flat. One power itera on leads to a near-op mal error but the op mal error remains high. In fact, a 1% rela ve error is not reachable for ℓ < 0.5 for a maximum compression rate of 10%.

Matérn k 5

 5

Figure 7 . 6 :

 76 Figure 7.6: Accuracy of the fixed rank rSVD w.r.t. the compression rate for a Matérn ν = 2.5 correla on and n = 2.000 par cles distributed on the unit sphere. Le : Op mal rela ve error in Spectral and Frobenius norms w.r.t. to r = 1, . . . , n. Right: Accuracy of the rSVD w.r.t. prescribed rank r = 20, . . . , 200. Observa ons: The correla on is smoother than the

5 Figure 7 . 7 :

 577 Figure7.7: Accuracy of the fixed rank H 2 -powered rSVD w.r.t. the compression rate for n = 2.000 par cles distributed on the unit sphere. Le : For a Gaussian kernel with ℓ = 0.5 the accuracy of the rSVD can be very high even below 10% compression. However, the ufmm only allows for represen ng the covariance matrix up to a fixed accuracy. For a given interpola on order p, the H 2 powered variant has the same accuracy as the dense variant then slightly reaches a plateau. Right: For the exponen al kernel, even with ℓ = 1.0 the accuracy of the rSVD is rela vely low, therefore the ufmm never disturbs the method in this range of compression and p can be set to the lowest value.

Figure 7 . 8 :

 78 Figure 7.8: Schema c view of the FMM with mul ple (r) right-hand-sides.

Figure 7 . 9 :

 79 Figure 7.9: Theore cal memory requirements for the nearfield computa on in the worst case scenario (no empty cell).

2 O ( n 2

 22

MKLFigure 7 . 10 :

 710 Figure7.10: Running me (s) per matrix-to-vector product w.r.t. the number of par cles n for a total of r = 10 vectors and target accuracy of about 10 -3 . In order to compare all algorithms we use a smooth kernel, namely the Gaussian kernel with ℓ = 0.5. Par cles are distributed uniformly in the unit cube (le ) and on the unit sphere (right). The tree depth L of the ufmm (wri en below blue squares) ensures a rela vely constant average number of par cles per leaf: n 0 = 74 (le ) and n 0 = 62 (right). Machine: plafrim/riri -Deca-core Intel Xeon E7-4870 @ 2.40GHz with 1TB ram and 30MB L3 Cache.

3 Figure 7 . 11 :

 3711 Figure7.11: Time for compu ng a randomized SVD using either the fixed rank or fixed accuracy algorithm with q = 0 and s = 10. Matrix mul plica ons are computed either in a dense fashion (filled squares) or by means of the ufmm (empty squares) or smooth-ufmm (empty circles) with p = 3, par cles are randomly distributed on the unit sphere and the correla on is Gaussian with ℓ = 0.5. Dense matrix mul plica ons, SVD and QR Decomposi on are performed using sequen al MKL Blas rou nes.

Figure 7 . 12 :

 712 Figure 7.12: Realiza on of a Gaussian random field on n = 72.000 points distributed on the unit sphere using a Gaussian correla on with ℓ = 0.25 (le ) and ℓ = 0.50 (right). The approximate square root was obtained by mean of a smooth-ufmm- accelerated rSVD with a fixed accuracy of ε = 10 -2 .

Figure 7 . 13 :

 713 Figure7.13: Realiza on of a Gaussian random field on n = 72.000 points distributed on the unit sphere using a Gaussian correla on. The approximate square root was obtained by mean of a ufmm-accelerated rSVD with a fixed accuracy of ε = 10 -2 (le ) and 10 -1 (right). Observa ons: The differences in terms of covariance quality and smoothness of the field are hardly

L 6 -L 6 -

 66 n = 1.10 3 L 8n = 1.10 3 L 6n = 5.10 3 L 8n = 5.10 3 L 6n = 10 4 L 8n = 10 4 n = 1.10 3 L 8n = 1.10 3 L 6n = 5.10 3 L 8n = 5.10 3 L 6n = 10 4 L 8n = 10 4

Figure 8 . 1 :

 81 Figure 8.1: Best rank-r error w.r.t. the rank in term of Frobenius (right) and Spectral (le ) norm for a subset of a full sample from Lake Geneva at various mes. S r is obtained by mean of a full SVD, therefore we can only afford to compute small subsamples, namely up to n = 10 4 .

8 Figure 8 . 2 :

 882 Figure 8.2: Le : Eigenvalues of the matrix associated with the first reads of sample L 1 (top) and L 3 (bo om) with highlight on the nega ve eigenvalues (blue diamonds). Right: Distor on (top) and contribu on of the nega ve eigenvalues to the Frobenius norm (bo om) w.r.t. number of components or rank r of the Truncated SVD (tSVD) for various sample L t . Observa on:The correla on between nega ve eigenvalues and distor on is illustrated by the regime r = 10, . . . , 100. In fact, the distor on first decreases with the rank then as nega ve eigenvalues occur the distor on is affected and reaches a plateau before decreasing fast again.

  largest eigenvalues: λ k n = {1, 2, 3} .10 4 and s = r n = 1.10 4r = 1000 n = 1.10 4r = 200 n = 2.10 4r = 200 n = 3.10 4r = 200

Figure 8 . 3 :

 83 Figure 8.3:Le : First 500 eigenvalues of a random subsample of L 6 compared to approximate values obtained by rSVD without (s = 0) or with (s = r/2) oversampling. Right: First 200 approximate eigenvalues obtained by rSVD with significant oversampling s = r on larger subsamples: n = 10.000, n = 20.000 and n = 30.000. We use the result of a rSVD with r = 1.000 as

Figure 8 . 4 :

 84 Figure 8.4: Analysis of the coherence of the similarity matrix associated with 10.000 reads from the sample L 6 (le ) and

Sample L 1 ,

 1 Fraction of the energy: εr Sample L 1 , rSVD with s = r rSVD with s = r full sample n = 3.10 4 n = 1.10 4

Figure 8 . 5 :

 85 Figure 8.5: Analysis of the point cloud distor on w.r.t. the input rank r of the rSVD for 10.000 reads from the sample L 1 or the full sample. Le : Contribu on of the nega ve eigenvalue to the energy of the approximated similarity. Right: Distor on of the point cloud. Observa ons: The contribu on of the nega ve eigenvalues to the energy of the representa on is s ll around

Sample L 1 ,Figure 8 . 6 :

 186 Figure 8.6:Running mes of the rSVD (in seconds) w.r.t. the input rank r for the full L 1 sample (n ≈ 70.000) and a subsample with n = 10.000. The plot shows the me required to read the distance and compute the similarity matrix S (do ed lines), the running me of the rSVD (solid lines) and the me required to compute the distor on δ (dashed lines).

Figure 8 . 7 :

 87 Figure 8.7: Representa on of the point cloud associated to 10.000 random reads of the sample L 6 on the first 3 features F 1, F 2, F 3 using 2D plots. The components were computed using rSVD with s = r = 200. Iden fied species are represented using large colored markers, while unknown species are represented in black. Observa on: Such representa on allows for easily verifying that individuals of the same iden fied species gather in clusters. However random subsampling does not allow to represent all individuals and in par cular all the previously iden fied ones.

Figure 8 . 8 :

 88 Figure 8.8: Representa on of the point cloud associated the full sample L 6 on the first 3 features (F 1, F 2, F 3) using 2D plots.The components were computed using rSVD with s = r = 200. Iden fied species are represented using large colored markers, while unknown species are represented in black. Observa on: Such representa on allows for easily verifying that individuals of the same species gather in clusters. Moreover, individuals sharing the same root like the Encyonema are themselves located in cluster or at least small bounded domains. Finally, in such low dimensions, the cloud has a similar shape than the subsample with n = 10.000 although the distor on is supposed to be a lot higher because of the larger size and the smaller input rank of the rSVD.

Figure 8 . 10 :

 810 Figure 8.10: Popula on of the point clouds associated to various samples represented on the first 3 features (F 1, F 2, F 3) using 2D plots. The components were computed using rSVD with s = r = 50 on the full samples: L 1 (top le ), L 2 (top right), L 3 (bo om le ) and L 4 (bo om right). Observa on: The shape of the cloud varies significantly along the me, while s ll

Figure 8 . 11 :

 811 Figure 8.11: Representa on of the point cloud associated to the full sample L 6 on the first 10 features (F 1, . . . , F 10) using parallel coordinates from matplotlib. The components were computed using rSVD with s = r = 50.

Contents B. 1 B. 2 B. 1 .

 121 Geometries with various heterogeneity . . . . . . . . . . . . . . . . . . 142 Statistics on the octree . . . . . . . . . . . . . . . . . . . . . . . . . . 142 Geometries with various heterogeneity

Figure B. 1 Figure B. 1 :

 11 Figure B.1 represents actual distribution of points for all geometries considered in this thesis. It also shows associated 2D distributions inside a quadtree with L = 3, 4, 5. This gives an estimation of the filling of the tree in each case, as well as the number of nearfield and farfield interactions. 3D grid L = 3 L = 4 L = 5

2 Filling 5 Figure B. 2 :

 252 Figure B.2: Le : Filling rate at the leaf level w.r.t. to tree depth L for various geometries. When the points are distributed on

Figure B. 3 :

 3 Figure B.3: Total and average number of leaf-level M2L and P2P operators for various distribu ons of par cles (n = 10 6 ) w.r.t.L. Observa ons: The average number of operators stays rela vely constant for all geometries, in par cular we quickly reach the maximum number of operators for the unit cube, namely 27 P2P and 189 M2L, while the unit sphere only has 15 P2P and 50 M2L in average. The prolate sphere exhibit 2 regimes, namely L ≤ 4 or L ≥ 5, whether farfield interac ons are considered

  93GHz with 96GB ram and 12MB L3 Cache. Results are presented on Figure C.1 for the cube, Figure C.1 for the unit sphere and Figure C.1 for the prolate sphere. They show the computational time of a fast multipole matrix multiplication per rhs for a total of 10 rhs.

Cube

2 4 smooth 2 Figure C. 1 :

 2421 Figure C.1:Running me of a mul -rhs FMM w.r.t. the rela ve error magnitude using various algorithms: ufmm and smoothufmm with p = 2 . . . 11, and global with p = 6 . . . 14. We used 72.000 par cles randomly distributed in a cube.

2 6 smooth 2 Figure C. 3 :

 2623 Figure C.3:Running me of a mul -rhs FMM w.r.t. the rela ve error magnitude using various algorithms: ufmm and smoothufmm with p = 2 . . . 11, and global with p = 6 . . . 14. We used 72.000 par cles randomly distributed on a prolate sphere (ra o 1:1:10). Observa ons: ufmm and smooth-ufmm have approximately the same cost when they share the same tree depth. The smooth-ufmm is significantly faster if we choose an op mal tree depth, e.g., L = 4 represented in brown. If the Gaussian decreases sufficiently slow (i.e., the rank is sufficiently low), then the global performs rela vely well compared to the hierarchical variants.

3 5 =

 5 243 to n = n A + n B -n AinterB = 18 × 3 + 6 × 6 -(6 × 3 + 6) = 66. In fact,

  3) of the interpolant k in (2.8) we get an approximation of the elastic energy field σ, that reads asE((C)) ≈ Ẽ((C)) = -µ 8π ( ẼA + 2 1 -ν ( ẼB -ẼC + ν ẼD )) (D.4)where contributions ẼA , ẼB , ẼC and ẼD are defined asẼA = ∑ |α|≤p+1 (C)S α (x)b i dx i

S

  α (x)b j t i (L C,D α ) ji dx (L2P-ED)All terms share the same multipole expansions, namelyM β = C S β (y)b ⊗ tdx ′ , (D.6)while the L2P, resp. M2L, step differs for each term except E B and E C , resp. E C and E D .

9

  local expansions in order to compute E A , E C and E D at the L2P step. For E B , we need to apply each R ,ij to the M jk expansions, i.e., (3 × 3) × 3 = 27 M2L operations, and store the results in (3 × 3) = 9 local expansions. Let us use the superscript A for all operators used for computing E A,C,D and B for E B . Hence, i, j = 0 . . . 2 the local expansions reads asL A ij = R kk M ij and L B jk = R ij M jk If we denote I ij ∈ [0, 6] the index mapping from [0, 2] × [0, 2]to the set of symmetric components indices plus one index for the trace of R, namely R kk , then L A I = R 6 M I and

  L B I and L ′ 9+I = L A I = R 6 M I , I = 0 . . . 8 Mura's variant Let us recall the equation of the energy derived from Mura's formula

[x 1 , x 2 ]

 12 of a larger interval [a, b] ∫ [x 1 ,x 2 ] f (x)dx ≈ I p = ∫ [x 1 ,x 2 ] p ∑ m=0 S p (Φ -1 (x), xm )f (Φ(x m ))dx =

2 X

 2 + c S ) -c B ), xm ) L 2 dX

  2 and Table 1.1 gather detailed theoretical estimations of its complexity and memory footprint.

	Interpolators During the P2M/L2P steps, part of the expansion can be precomputed, e.g.,
	the polynomial basis, and stored at a O(p) cost. The rest of the expansion is computed on the
	fly at a O(p 3 ) cost. The M2M/L2L step can be slightly reworked in order to be performed in
	O(p 4 ) operations instead of the naive variant, that involves O(p 6 ) operations and a memory
	footprint of O(p 3 ). The optimized variant only requires the interpolator S to be assembled
	on a single fictitious pair of child and parent 1D interpolation grid at each level, which
	results in a memory footprint of O(p 2 ).
	M2L operators In the framework of interpolation based FMM, the (p + 1) 3 -by-(p + 1) 3 matrices K = { Kαβ } |α|,|β|≤p are called the M2L operators. These operators are used to
	transfer multipole expansions M, i.e., expansions w.r.t. the sources y, to local expansions
	L, i.e., expansions w.r.t. the targets x. The M2L operators are known explicitly, therefore
	they are precomputed at each level for each possible pair of interacting cells. Their number
	reduces to a maximum of 316 per level in 3D thanks to the translation invariance. The O(p

Table 1 .1: Asympto

 1 c memory footprint of each step of the bbfmm w.r.t. the interpola on order p. In order to get the exact number of bytes one needs to mul ply by the value type size, i.e., 4 for single floats and 8 for double floats. The complexi es of the non-homogeneous M2L and the M2M/L2L is mul plied by the number of levels, namely L -1.

			Before Computation		During Computation
	Step	P2M/L2P	M2M/L2L		M2L		global
				homogeneous non-homogeneous
	Operator	T (x)	S(x, x′ )	K		Kref	M/L
	# operators	8 L	3 × 8( L -1)	316( L -1)	316	n cells
	# entries	(p + 1)	(p + 1) 2	(p + 1) 6	(p + 1) 6	(p + 1) 3
							For the P2M/L2P the
	complexi es are mul plied by the number of leaves. Finally the order-p 3 expansions stored during the computa on is mul plied by the number of non-leaf cells n cells = ∑L -1 L=2 2 dL = (8 L -8 2 )/7.
		Step	P2P	P2M/L2P	M2M/L2L	M2L
		# flops 27 × 8 L × n 2 0	n × p 3	8n cells × 3p 4 189(n cells + 8 L) × p 6

Table 1 .2: Asympto

 1 c complexity of each step of the bbfmm w.r.t. the interpola on order p. The complexity of the M2L is mul plied by the number of cells n cells + 8 L and the complexity of the M2M/L2L by the number of non-leaf cells n cells .

Algorithm 4: compressed-bbfmm functions

	Function precomputeM2L()
	// For all levels in tree, or 1 level if kernel is homogeneous
	for L = 2, . . . , L do
	// Assemble all possible interactions, i.e., 316.
	for i-th source cell C y ∈ I f ull (C (L) x ) do (L)
	Ki αβ (L) = k(x α , ȳβ ), ∀α, β/|α|, |β| ≤ p
	// Compress the entire set of interactions.

Table 5 . 1 :

 51 Asympto c complexi es of the ufmm variants compared to the bbfmm. The complexi es of the M2L and M2M/L2L are given per level and per non empty cell, they should be mul plied by the number of non-empty cells (max. 8 L ) and then summed over levels. For the P2P and P2M/L2P the complexi es are given globally. The constant n 0 = n/8 L denotes the average number of par cles per leaf.

	n 0

  Distributions of loops and visualization of the fields . . . . . . . . 6.4.2 Accuracy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.4.3 Sequential performance . . . . . . . . . . . . . . . . . . . . . . .

	6.4.1

Contents 6.1 Fast Multipole computation of the isotropic stress and forces . . . . . . 6.1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.1.2 Interpolation of the elastic stress . . . . . . . . . . . . . . . . . . 6.1.3 A naïve Fast Multipole summation scheme . . . . . . . . . . . . . 6.2 Improved farfield computation of the stress . . . . . . . . . . . . . . . 6.2.1 Considering symmetries . . . . . . . . . . . . . . . . . . . . . . . 6.2.2 Refactoring interactions . . . . . . . . . . . . . . . . . . . . . . . 6.2.3 Shifting derivatives . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2.4 Comparative analysis of the theoretical cost . . . . . . . . . . . . 6.3 Implementation details . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3.1 The extended bounding box concept . . . . . . . . . . . . . . . . 6.3.2 Numerical integration over segments . . . . . . . . . . . . . . . . 6.3.3 Enforcing homogeneity . . . . . . . . . . . . . . . . . . . . . . . . 6.4 Numerical benchmarks . . . . . . . . . . . . . . . . . . . . . . . . . . 6.1. Fast Multipole computation of the isotropic stress and forces

Table 6 . 1 :

 61 Number of FMM operators required for the computa on of the stress and the energy, that represent the constants in front of the theore cal memory requirements for each approach.quently, the variant using shifted derivatives should be the fastest one in theory, since it has the lowest constant. The cost of the P2M and L2P steps is similar for all optimized variants, except that the P2M of the symmetric variant is far more expensive than that of the other optimized variants.

	Quantity Approach L2P	M2L	P2M
	Stress	Direct	9	27 + 27 = 54	27
		Symm	6	18 + 18 = 36	27
		Facto 1	6	46	9
		Facto 2	6	18 + 33 = 51	9
		Shift	6	6 + 6 = 12 6 + 1 = 7
	Energy	Direct	19	47	9

Table 6 .2: Constants

 6

in front of the theore cal complexity for each approach w.r.t. to the FMM steps.

  real ∥C real -C∥ 2,rel roc 2 ∥C real -C∥ ∞,rel roc ∞

						n 1 • 10 3	2.48 • 10 -1		1.57 • 10 -1
		88 • 10 -2	0.5	4.53 • 10 -2	0.51	1 • 10 4	7.88 • 10 -2	0.5	5.68 • 10 -2	0.44
	1 • 10 5	2.63 • 10 -2	0.48	1.57 • 10 -2	0.46	1 • 10 5	2.67 • 10 -2	0.47	1.84 • 10 -2	0.49
	1 • 10 6	1.02 • 10 -2	0.41	8.53 • 10 -3	0.27	1 • 10 6	1.09 • 10 -2	0.39	1.09 • 10 -2	0.23
					ℓ = 0.50
			ufmm				smooth-ufmm

Table 7 . 1 :

 71 real ∥C real -C∥ 2,rel roc 2 ∥C real -C∥ ∞,rel roc ∞ Error (and rate of convergence, roc) on the covariance matrix C (n = 72.000) w.r.t. the number of realiza ons n real for a Gaussian correla on with ℓ = 0.25 (top) and ℓ = 0.50 (bo om). The approximate square root used to generate realiza ons is computed using a fixed accuracy randomized SVD (ε = 10 -2 ) powered by ufmm (le ) or smooth-ufmm (right) with p = 3. The sample covariance C real is computed from the n real realiza ons as(7.15). Observa ons: The convergence rate is close to the theore cal rate, i.e., 1/2. The rate slightly decreases as n grows, since C fails to represent C with sufficient accuracy, though the sample covariance accuracy is already very sa sfying, namely below 10%.

						n 1 • 10 3	1.32 • 10 -1		1.21 • 10 -1
		75 • 10 -2	0.54	4.56 • 10 -2	0.39	1 • 10 4	3.72 • 10 -2	0.55	3.95 • 10 -2	0.48
	1 • 10 5	1.41 • 10 -2	0.43	1.76 • 10 -2	0.41	1 • 10 5	1.39 • 10 -2	0.43	1.35 • 10 -2	0.47
	1 • 10 6	4.72 • 10 -3	0.47	5.95 • 10 -3	0.47	1 • 10 6	4.45 • 10 -3	0.49	5.69 • 10 -3	0.37

Table 8 .

 8 2 shows the size of the different samples involved in this chapter, that go from about 70.10 3 to 140.10 3 reads. Therefore, samples can have quite different sizes. We expect the point clouds to have different shapes while always present clusters and comparable numerical rank.

	sample name	L 1	L 2	L 3	L 4	L 5	L 6
	#sequences	72083 98492 72897 136450 75218 99594

Table 8 . 2 :

 82 Size of diatoms samples from Lake Geneva at various me of the year. Observa ons: The number of sequences, a.k.a., reads, can vary by a factor of 2, while the average number of reads is approximately 100.000.

  ACA Adative Cross Approximation. 27 ARRF Adaptive Randomized Range Finder. 40, 95 BBF structured Block Basis Factorization. 43 DD Dislocation Dynamics. 4, 47, 129, 131, 167 EVD Eigen Value Decomposition. 39, 109, 110, 116, 117, 130 FFT Fast Fourier Transform. 25, 27, 28, 52, 54, 57, 59, 90, 99, 129, 165 FJLT Fast Johnson Lindenstrauss Transforms. 42, 110 FMA Fast Multipole Algorithm. 10-14, 16, 173 FMM Fast Multipole Method. 3-5, 10, 11, 14, 16, 17, 20, 21, 23-25, 35, 42, 43, 46-48, 52, 58, 60, 61, 66, 67, 70, 71, 74, 76-80, 90, 94-100, 102, 106, 129-132, 136, 137, 146-148, 154, 155, 165, 167-169, 171, 173 GRF Gaussian Random Field. 5, 47, 88-90, 92, 95, 129 GRV Gaussian Random Variable. 88, 168 LRA Low-Rank Approximation. 4, 5, 38-43, 46-48, 103, 110, 113, 116, 129, 130 MDS Multidimensional Scaling. 48, 106, 107, 109-111, 113, 115-117, 120, 125, 130 MEKA Memory Efficient Kernel Approximation. 43 NGS New Generation Sequencing. 106, 110 NLA Numerical Linear Algebra. 3-5, 132 RSI Randomized Subspace Iterations. 40, 41, 95 rSVD Randomized SVD. 38-41, 90-95, 98-100, 102, 113, 116-126, 168, 169 SRFT Subsampled Randomized Fourier Transform. 42, 43 SRHT Subsampled Randomized Hadamard Transform. 42, 43 Acronyms SVD Singular Value Decomposition. 27, 38-41, 43, 48, 90, 92, 93, 98, 100-103, 108-113, 115, 131 tSVD Truncated SVD. 114, 117

* See OptiDis code's home page at http://optidis.gforge.inria.fr/

* This theory is obtained by spreading the Burger's vector smoothly over the core of the dislocation.

* The first value is real and the rest is complex. Moreover, the second half of the rest is the conjugate of the first half.

Acknowledgments

https://gitlab.inria.fr/solverstack/ScalFMM.

)

(

# particles: n application time (s) )

(2)

(3)

# particles: n application time (s) )

(3)

# particles: n application time (s) -5 . In order to compare all algorithms we use a smooth kernel, namely the Gaussian kernel with ℓ = 0.5. Par cles are distributed uniformly in the unit cube (le ), on the unit sphere (center) and on the prolate sphere (right).

Conclusion

We introduced a new interpolation based FMM with efficient formulation of the M2L step in Fourier domain. Please refer to our research report [START_REF] Blanchard | Fast hierarchical algorithms for generating gaussian random fields[END_REF] for a shorter introduction. The L 6

Appendices

A

Tensorial interpolation based FMM

A.1 Nature of tensorial interactions

Although the Fast Multipole Method (FMM) was originally designed for particle interactions with k(r) = 1/r [START_REF] Greengard | A fast algorithm for particle simulations[END_REF], i.e., a scalar valued kernel, many n-body problems formulated in terms of Boundary Integral Equations [START_REF] Bonnet | Boundary integral equation methods for solids and fluids[END_REF][START_REF] Nishimura | Fast multipole accelerated boundary integral equation methods[END_REF] have arisen in the past decade that involve tensorial valued kernels [START_REF] Traub | A kernel interpolation based fast multipole method for elastodynamic problems[END_REF][START_REF] Schanz | Wave propagation in viscoelastic and poroelastic continua: a boundary element approach[END_REF][START_REF] Chaillat | Fast Multipole Method for 3-D elastodynamic boundary integral equations[END_REF] and that can be very well treated by the FMM. The tensorial nature of the kernel appears as soon as we consider that densities are multivariate quantities evaluated on grid points x, i.e.,

where d w denotes the dimension of w(x). For instance, in most BEM formulations used for isotropic elasticity the Green's function is given as

for the infinite elastic space [START_REF] Chaillat | Fast Multipole Method for 3-D elastodynamic boundary integral equations[END_REF], where µ and ν denote the Lamé coefficients, whereas for acoustics it is given as a scalar valued function,

namely Laplace kernel, or Helmholtz kernel in the frequency domain [START_REF] Messner | Fast Boundary Element Method in Acoustics[END_REF][START_REF] Sylvand | La méthode multipôle rapide en électromagnétisme : Performances, parallélisation, applications[END_REF]. In elasticity problems, G is tensorial because it represents the relation between a multivariate density w (namely a point force) and the displacement u. This tensorial nature can become problematic for problems involving large number of parameters, e.g., viscoelasticity and poroelasticity [START_REF] Schanz | Wave propagation in viscoelastic and poroelastic continua: a boundary element approach[END_REF], highly oscillatory kernels, e.g., elastodynamics in time or frequency domain [START_REF] Traub | A kernel interpolation based fast multipole method for elastodynamic problems[END_REF], or multi-layered spaces [START_REF] Chaillat | A new fast multipole formulation for the elastodynamic half-space fundamental solution[END_REF]. In particular, the author of the present contributed to the design of an efficient interpolation based directional Fast Multipole BEM formulation for elastodynamics in Laplace domain and results are discussed in Thomas Traub's PhD thesis [START_REF] Traub | A kernel interpolation based fast multipole method for elastodynamic problems[END_REF]. In the present thesis we used tensorial interactions in the context of Dislocation Dynamics (DD) simulations (see Chapter 2), a special case of elastostatic interactions between line segments.

A.2 A generic formalism

Let us first define a generic formalism for tensorial interactions using a direct formulation. Then, we describe the associated fast multipole summation scheme. For the sake of clarity and genericity we denote d the ambient dimension, although most of the time d = 3.

A.2.1 Direct formulation

First of all, multivariate potentials (resp. densities) are written in bold characters as p (resp. w) and defined by d p (resp. d w ) components identified by a multi-index I (resp. J) of order o p (resp. o w ), such that |I| = max i≤op (I i ) ≤ d and d w = d ow . Then, the kernel of interaction is denoted k and defined for any grid points x and y as

Finally, the I-th component of the potential generated at the target point x by all components of the densities located at the source points y reads as

where

A.2.2 Fast Multipole formulation

Computing such interactions on a large number of points can be done in O(n) operations using FMM. The steps of the associated Fast Multipole summation scheme are defined by the following equations:

• P2M: Aggregate contributions of the d w components of density w located at source points y into d w multipole expansions defined on source interpolation nodes ȳβ

A.4 Optimizations and numerical complexities A.4.1 Multidimensional interpolators

The L2P, resp. P2M, steps must be applied for d p , resp. d w , potentials, resp. density. In a naïve implementation, this would result in multiplying the cost of these steps by their associated dimension. Most importantly, the steps described above require that d w multipole expansions and d p local expansions are stored during the algorithm. We recall that storing an expansion costs O(p 3 ) bytes, more precisely 8 × (p + 1) 3 for a real valued expansion in double precision. Each cell of the octree contains at least * 2 expansions, therefore it must be multiplied by the number of cells, e.g., # cells (L) = ( 23 ) L = 8 L for an arbitrary level L. Hence, the cost of storing the expansions grows dramatically with the depth of the octree.

A.4.2 Tensorial M2L

The tensorial nature of the interactions also affects the critical M2L step. In fact, all components k IJ need to be evaluated, stored and applied to the local expansions. Therefore, the cost of the precomputation is roughly multiplied by the number of components d k . If some components depend on each other by a constant factor (e.g., symmetric tensors), then less components need to be stored. We denote d min k the number of independent components of the tensor of interactions k that actually need to be stored. Finally, in (A.7) the summation over J is done outside the loop over the interpolation grid in order to preserve the optimized application of the M2L operators.

C.3 Performance of the global fft

The global fft will always have the same cost at a given interpolation order, since this method is oblivious of the shape of the distribution. On the other hand, the cost of the hierarchical methods may vary significantly from one distribution to another. Let us for instance consider a Gaussian correlation with ℓ = 0.5. If the particles are distributed in the cube (i.e., an homogeneous distribution) the cost of the global fft lies between the ufmm and the smooth-ufmm with optimal L, see Observa ons: ufmm and smooth-ufmm have approximately the same cost when they share the same tree depth. The smoothufmm is significantly faster if we choose an op mal tree depth, e.g., L = 3 represented in brown. If the Gaussian decreases sufficiently slow (i.e., the rank is sufficiently low), then the cost of the global is slightly lower than the ufmm but the op mal smooth-ufmm s ll performs be er. The associated octree has larger depths ( L = 6) than that of the unit sphere ( L = 4), however in the lowest level a large number of cells are empty. Consequently, the running times are only slightly larger than for the unit sphere.

C.4 Other distributions

• While the computation of σ B ij for all (i, j) ∈ {1, . . . , 3} 2 with i < j involves the following

• The 24 components shared by σ A ij and σ B ij are colored in grey.

We could not recover the factor of 6 mentioned in [8], i.e., the 40 components out of 3 5 = 243 using this approach. Even considering σ A ij + σ A ji instead of σ A ij does not seem to change the result. Moreover, Arsenlis et al. [8] affirms that none of the elements of G iiklm , G ijkli and G ijklj are required in the computation, but we found that it is only true for the computation of σ B . On the other hand, G ijill is indeed not involved in the computation of either σ A or σ B . The aforementioned paper is known to contain a few typos, therefore we developed our own approach. Besides we do not think this is the most relevant way of optimizing M2L operations, especially when considering polynomial interpolation.

D.1.2 Refactoring M2L operations

Principle In order to define an efficient summation scheme, we want to identify the optimal definition of the multipole expansions that leads to the minimum number of M2L operations. Our method consists in rewriting the full expression of each component of σ and refactoring it w.r.t. to the components of R. Let us first denote G ijklm = G IJ with

Computing σ A ij and summing with the transposed matrix may be a good idea in order to obtain

Similarly, for σ B we have

Finally, the expression of the elastic stress field σ

Using the previous forms of σ A and σ B the applications of the M2L operators are optimal in the sense that the minimum number of R ... -to-m .. products is performed, namely 46, and the minimum number of multipole expansions is stored, namely the 9 components of m. Until now this is the cheapest method in term of memory requirements, however it requires assembling the term in front of each component of R just before the application of the M2L operator, i.e., approximately 189 times per cell for a full octree.

Final formulation A good compromise can be found in order to accelerate the M2L step, while minimizing the memory requirements. First of all, if we transfer the contributions to