## A Fast Summation Method for translation invariant kernels

Fabien Casenave ∗

August 25, 2015

## Abstract

We derive a Fast Multipole Method (FMM) where a low-rank approximation of the kernel is obtained using the Empirical Interpolation Method (EIM). Contrary to classical interpolationbased FMM, where the interpolation points and basis are fixed beforehand, the EIM is a nonlinear approximation method which constructs interpolation points and basis which are adapted to the kernel under consideration. The basis functions are obtained using evaluations of the kernel itself. We restrict ourselves to translation-invariant kernels, for which a modified version of the EIM approximation can be used in a multilevel FMM context; we call the obtained algorithm Empirical Interpolation Fast Multipole Method (EIFMM). An important feature of the EIFMM is a built-in error estimation of the interpolation error made by the low-rank approximation of the far-field behavior of the kernel: the algorithm selects the optimal number of interpolation points required to ensure a given accuracy for the result, leading to important gains for inhomogeneous kernels.

## 1 Introduction

We consider the following problem: compute approximations of many sums of the form

$$f ( \bar { x } _ { i } ) = \sum _ { j = 1 } ^ { N } \sigma _ { j } K ( \bar { x } _ { i }, \bar { y } _ { j } ), 1 \leq i \leq N,$$

for a set of potentials σ j , and target points ¯ x i ∈ Ω and source points ¯ y j ∈ Ω, where Ω ⊂ R 3 , with a controllable error glyph[epsilon1] . The domain Ω and N are fixed for all computations, but the potentials σ j and the points ¯ x i and ¯ y j can vary from one computation to another. The many-query situation can come from a temporal N-body problem (in which case the set of source points equals the set of observation points), or the iterative resolution of a dense linear system.

Various algorithms have been proposed to tackle such problems. Some methods only using numerical evaluations of K have been proposed, see [1, 8, 15]. Many methods, e.g. the mosaic skeleton method [27], hierarchical matrices [20], the Fast Multipole Method (FMM) [17] (and its algebraic generalization H2-matrices [18, 19]), the panel clustering method [21], rely on local low-rank approximations of K :

$$K ( x, y ) & \approx \sum _ { l, m = 1 } ^ { d } \Delta _ { l, m } g _ { l } ( x ) h _ { m } ( y ), x \in I, y \in J,$$

where I and J are subsets of Ω. Such approximations can be analytic, like interpolation (for instance the Chebyshev interpolation as in [14, 20]), or multipole expansions. FMM for general kernels have been developed leading to kernel-independent procedures, see [24, 28]. In [13], another kernelindependent FMM algorithm, based on the Chebyshev interpolation to obtain the low-rank approximation (2) of K , is developed and called Black-Box FMM. Since the algorithm we present here is in the same family as the one developed in [13], our numerical experiments contain comparisons with this method. For other kernel-independent interpolation based procedures, see [10, 11]. Adaptive Cross Approximation (ACA) [3, 4] directly constructs approximations from the evaluations of the kernel on the target and source points K x , y (¯ i ¯ ), by selecting some columns and rows of the matrix ( j K x , y (¯ i ¯ )) j i,j .

∗ IGN LAREG, Univ Paris Diderot, Sorbonne Paris Cit´ e, 5 rue Thomas Mann, 75205 Paris Cedex 13, France.

Approximating a matrix by extracting some of its rows and columns has been investigated in [16]. With ACA methods, (2) is such that g l depends on J and h m depends on I , which leads to a storage complexity of dN log( N ). With FMM and H2-matrices, g l and h m are independent of respectively J and I , which allows additional factorization and the use of nested basis to lower the storage complexity to dN . An improved ACA exploiting the nested basis property has been proposed in [5].

In this work, a low-rank approximation of K is obtained with the Empirical Interpolation Method (EIM), a procedure that represents a two-variable function in a so-called affine dependent formula and that is used in the Reduced Basis community, [2]. The EIM provides a non-linear and functiondependent interpolation: the basis functions and interpolation points are not fixed a priori but are computed during a greedy procedure, which differs from the Chebyshev interpolations. The procedure should capture any irregularity of the kernel better than a fixed basis interpolation procedure, since the basis functions of the EIM are based on evaluations of the kernel itself. Theoretical results on the convergence of EIM interpolation error in the general case are limited, to the author knowledge, to an upper bound of 2 d -1 of the Lebesgue constant, where d is the number of interpolation points, see [2, Proposition 3.2]. Various numerical studies illustrate the efficiency of the EIM [12, 22, 23, 26], even though a theoretical analysis validating this efficiency is not available yet. Depending on the kernel (and on the level in the tree for inhomogeneous kernels), the EIM selects an optimal number of terms that ensures a given accuracy of the approximation. Such a variable-order interpolation procedure have been used in [6] in the H2-matrix context with Lagrange polynomials.

A direct application of the EIM algorithm to the kernel K to obtain a local low-rank approximation (2) is such that g l depends on J and h m depends on I , and therefore is not compatible with a multilevel FMM: a nested basis procedure cannot be derived. To correct this, we make use of the translation-invariant hypothesis of K and derive a local-global low-rank EIM approximation of K , in the sense that one variable is restricted to a local small subdomain while the other variable is taken in almost all the domain. Then, we derive recursive formulae to produce a multilevel FMM for this modified approximation: this is detailed in Section 3.3. The modified approximation uses four summations, instead of two summations for the Black-Box FMM. However, we will see that the M2L step, the most expensive one, has the same complexity as the Black-Box FMM one with respect to the number on interpolation points. Our algorithm requires a precomputation step, during which the greedy step of the EIM at each level is carried-out, at well as compression of some operators. In our simulations, we did not notice any speedup on homogeneous kernels compared to Black-Box FMM. For inhomogeneous kernels, the approximation problem gets easier as we get higher in the tree: corresponding EIM require much less interpolation points and speedups are measured. In 10 6 points test-cases, depending on the kernel, the EIFMM gets more interesting if one has to compute more than 2 to 18 queries of the summation. In a 10 8 points test-case, for the kernel e -r 2 , the precomputation step becomes negligible, and a single query is computed twice as fast for a nearly one order of magnitude better accuracy.

In Section 2, the EIM is recalled. In Section 3 is presented the new FMM algorithm, called the Empirical Interpolation Fast Multipole Method (EIFMM): first, a four-summation formula based on the EIM is derived in Section 3.1, that can be used in a multilevel FMM procedure, then a monolevel version of the EIFMM is proposed in Section 3.2, and a multilevel one in Section 3.3. In Section 4, the overall execution time is reduced using some precomputation and classical compression of the M2L step. Finally, numerical experiments are presented in Section 5, and conclusions are drawn in Section 6.

## 2 The Empirical Interpolation Method

Consider a function K x,y ( ) defined over D x × D y assumed to be real-valued for simplicity, where D x , D y ⊂ Ω. Fix an integer d . The Empirical Interpolation Method provides a way to approximate this function in the following form:

$$K ( x, y ) \approx ( I _ { d } K ) ( x, y ) \coloneqq \sum _ { m = 1 } ^ { d } \lambda _ { m } ( x ) q _ { m } ( y ),$$

where λ m ( x ) is such that

$$\sum _ { m = 1 } ^ { d } B _ { l, m } \lambda _ { m } ( x ) = K ( x, y _ { l } ), \quad \forall 1 \leq l \leq d.$$

The functions q m ( ) · and the matrix B ∈ R d × d , which is lower triangular with unity diagonal, are constructed in Algorithm 1, where δ d = Id -I d and ‖· ‖ D y is a norm on D y , for instance the L ∞ ( D y )-or the L 2 ( D y )-norm. In practice, the argmax appearing in Algorithm 1 is searched over finite subsets of D x and D y , denoted respectively by D x, trial and D y, trial and called training sets. Notice that Algorithm 1 also constructs the set of points { y l } 1 ≤ ≤ l d in D y used in (4), and a set of points { x l } 1 ≤ ≤ l d in D x . These points are interpolations points for the approximation of K (see Property 2.1) and are different from the target and source points ¯ x , y i ¯ . j

## Algorithm 1 Greedy algorithm of the EIM

- 1. Choose d &gt; 1
- 2. Set k := 1
- 3. Compute x 1 := argmax x ∈D x ‖ K x, ( · ) ‖ D y
- 4. Compute y 1 := argmax y ∈D y | K x ,y ( 1 ) |

$$\begin{array} { c } 5. \text{ Set } q _ { 1 } ( \cdot ) \coloneqq \frac { K ( \lambda _ { 1 }, \cdot ) ^ { \nu } } { K ( x _ { 1 }, y _ { 1 } ) } \\. \text{ c} \sim D \quad. = 1 \end{array}$$

- 6. Set B 1 1 , := 1
- 7. while k &lt; d do
- 8. Compute x k +1 := argmax ( x ∈D x ‖ δ k K )( x, · ) ‖ D y
- 9. Compute y k +1 := argmax ( y ∈D y | δ k K x )( k +1 , y ) |

$$\begin{array} { c c c } 1 0. & \text{Set $q_{k+1}(\cdot) \coloneqq \frac{(\tilde{d}_{k} \tilde{K})\tilde{(x_{k+1},\cdot)}$} } { ( \delta _ { k } K ) ( x _ { k + 1 }, y _ { k + 1 } ) } \\ \text{C++ $r$} \end{array}. \begin{array} { c } \text{Set $q_{k+1}(\tilde{d}_{k} \tilde{K})\tilde{(x_{k+1},\cdot)}$} } { ( \delta _ { k } K ) ( x _ { k + 1 }, y _ { k + 1 } ) } \\ \text{C++ $r$} \end{array}. \begin{array} { c } \text{Set $q_{k+1}(\tilde{d}_{k} \tilde{K})\tilde{(x_{k+1},\cdot)}$} } \\ \text{C++ $r$} \end{array}.$$

- 11. Set B k +1 ,m := q m ( y k +1 ), for all 1 ≤ m ≤ k +1
- 12. k ← k +1
- 13. end while

The following assumption is made:

(H) The dimension of Span ( x ∈D x K x, ( · ) , ) is larger than d .

From [ (H) ], the functions { K x , ( l · ) } 1 ≤ ≤ l d are linearly independent (otherwise, ( δ k K )( µ k +1 , x k +1 ) = 0 for some k in Algorithm 1).

It is shown in [7] that the approximation (3) can be written in the form

$$\left ( I _ { d } K \right ) ( x, y ) = \sum _ { l, m = 1 } ^ { d } \Delta _ { l, m } K ( x, y _ { l } ) K ( x _ { m }, y ),$$

where ∆ = B -t Γ -1 , and where the matrix Γ is is upper triangular and constructed recursively in the loop in k of Algorithm 1 in the following way:

- · k = 1:

$$\Gamma _ { 1, 1 } = K ( x _ { 1 }, y _ { 1 } ),$$

- · k → k +1:

$$\Gamma _ { k + 1, k + 1 } & = ( \delta _ { k } K ) ( x _ { k + 1 }, y _ { k + 1 } ), \\ \Gamma _ { m, k + 1 } & = 0, & \forall 1 \leq m \leq k, \\ \Gamma _ { k + 1, m } & = \alpha _ { m }, & \forall 1 \leq m \leq k,$$

where the vector α is such that ∑ k m =1 B l,m α m = K x ( k +1 , y l ), for all 1 ≤ l ≤ k .

[First interpolation point]

[First basis function]

[Initialize matrix B ]

[( k +1)-th interpolation point]

[( k +1)-th basis function]

[Increment matrix B ]

[Increment the size of the decomposition]

The particular form (5) is the key to the use of the EIM approximation in a fast multipole context. We recall the interpolation property of I d K ; see [23, Lemma 1]:

Property 2.1 (Interpolation property) . For all 1 ≤ m ≤ d ,

$$\begin{cases} \, ( I _ { d } K ) ( x, y _ { m } ) = K ( x, y _ { m } ), & \text{for all } x \in \mathcal { D } _ { x }, \\ \, ( I _ { d } K ) ( x _ { m }, y ) = K ( x _ { m }, y ), & \text{for all } y \in \mathcal { D } _ { y }. \end{cases}$$

An important aspect is that the error made by the EIM approximation over the training set D x, trial ×D y, trial , denoted glyph[epsilon1] k and defined as glyph[epsilon1] k := max ( y ∈D y | δ k K )( x k +1 , y ) , is monitored during the greedy | procedure. In practice, a stopping criterion is used on glyph[epsilon1] k instead of on the number of interpolation points d , which was used to simplify the presentation of the EIM. Hence, depending on the regularity of the function K , the EIM will automatically select the number of interpolation points to certify the error over D x, trial ×D y, trial .

## 3 The Empirical Interpolation Fast Multipole Method

The FMM algorithm is based on a low-rank approximation of the far-field behavior of the kernel of the form (2). Such approximations enable some factorizations in the computation of key-quantities and the propagation of information through a tree structure, resulting in lowering the complexity of the computation of (1). For a general presentation of the FMM algorithm, see the original article [17, Section 3].

## 3.1 A suitable approximation

For the simplicity of the presentation, the algorithm is described in 1D. It can be readily extended in 2D and 3D. Consider a kernel K x,y ( ), ( x, y ) ∈ Ω × Ω, where Ω = ( -L 2 , L 2 ) is the complete domain on which the FMM computation is carried out. The kernel is supposed translation invariant, i.e. such that K x,y ( ) = K x ( -a, y -a ), ∀ a ∈ R . The domain Ω is partitioned in 2 κ intervals, κ ∈ N ∗ , of length 2 l κ , with l κ := 2 - -κ 1 L . Consider two intervals I and J , well-separated in the sense that min ( x,y ) ∈ × I J | x -y | ≥ 2 l κ . To anticipate the generalization to higher dimensions, I and J are called 'boxes' in what follows. We denote c I and c J the center of the boxes I and J . Consider an EIM approximation constructed over I × J :

$$K ( x, y ) \approx \sum _ { l, m = 1 } ^ { d } \Delta _ { l, m } ^ { I, J } K ( x, y _ { l } ^ { J } ) K ( x _ { m } ^ { I }, y ), \quad x \in I, \ \ y \in J,$$

where I and J superscripts have been added to underline the box-dependent quantities (as deduced from Algorithm 1). In (6), K x,y ( J l ) depends on J , the box containing y and K x ( I m , y ) depends on I , the box containing x : as stated in the introduction, this approximation is then not suitable for a multilevel FMM procedure.

Define now I 0 := { x -c ˆ J , x ∈ ˆ I } for all box ˆ in the partitioning of Ω and all box J ˆ well separated I from ˆ , J and define J 0 the box of length 2 l κ located at the center of Ω (notice that J 0 is not a box from the partitioning of Ω). In our setting, J 0 = ( -l κ , l κ ) and I 0 = ( -L + l κ , -3 l κ ) ∪ (3 l κ , L -l κ ), see Figure 1. In other words, J 0 is a source box translated to the center of the domain, I 0 is the union of possible target boxes, translated in the same fashion.

Since K x,y ( ) = K x ( -c J , y -c J ), with x -c J ∈ I 0 and y -c J ∈ J 0 , it is inferred that

$$K ( x, y ) \approx \sum _ { l, m = 1 } ^ { d } \Delta _ { l, m } ^ { I _ { 0 }, J _ { 0 } } K ( x - c _ { J }, y _ { l } ^ { J _ { 0 } } ) K ( x _ { m } ^ { I _ { 0 } }, y - c _ { J } ),$$

where the EIM is carried out on D x = I 0 and D y = J 0 . An important aspect is that I 0 and J 0 are fixed domains, independent from the boxes I and J containing respectively x and y . Consider the

Figure 1: Representation of the partitioning of Ω in 1D for κ = 3 and of the corresponding boxes I 0 and J 0 .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0210v3-figure-2803f6d4-cc25-4654-9b2f-e1f7158121e3.jpg)

first evaluation of the kernel in (7). There holds K x ( -c J , y J 0 l ) = K x ( -c I , y J 0 l + c J -c I ), where x -c I ∈ J 0 and y J 0 l + c J -c I ∈ I 0 . It is then inferred that

$$K ( x - c _ { J }, y _ { l } ^ { J _ { 0 } } ) \approx \sum _ { l ^ { \prime }, m ^ { \prime } = 1 } ^ { d } \Delta _ { l ^ { \prime }, m ^ { \prime } } ^ { J _ { 0 }, I _ { 0 } } K ( x - c _ { I }, y _ { l ^ { \prime } } ^ { I _ { 0 } } ) K ( x _ { m ^ { \prime } } ^ { J _ { 0 } }, y _ { l } ^ { J _ { 0 } } + c _ { J } - c _ { I } ), \quad \ 1 \leq l \leq m,$$

where a second EIM has been carried out on D x = J 0 and D y = I 0 (with the same number of terms d for simplicity). Notice that if the kernel is symmetric, ∆ I 0 ,J 0 = (∆ J ,I 0 0 ) , t x I 0 m = y I 0 m for all 1 ≤ m ≤ d , and y J 0 m = x J 0 m for all 1 ≤ m ≤ d . Injecting (8) into (7):

$$K ( x, y ) & \approx \sum _ { l ^ { \prime } = 1 } ^ { d } K ( x - c _ { I }, y _ { l ^ { \prime } } ^ { J _ { 0 } } ) \sum _ { m ^ { \prime } = 1 } ^ { d } \Delta _ { l ^ { \prime }, m ^ { \prime } } ^ { J _ { 0 }, J _ { 0 } } \sum _ { l = 1 } ^ { d } K ( x _ { m ^ { \prime } } ^ { J _ { 0 } }, y _ { l } ^ { J _ { 0 } } + c _ { J } - c _ { I } ) \sum _ { m = 1 } ^ { d } \Delta _ { l, m } ^ { I _ { 0 }, J _ { 0 } } K ( x _ { m } ^ { I _ { 0 } }, y - c _ { J } ), x \in I, y \in J. \\ \tau _ { n } ( \in \tau _ { n } ) \tau _ { n } ( \tau _ { n } ) \tau _ { n } ( \tau _ { n } ) \tau _ { n } ( \tau _ { n } ) \tau _ { n } ( \tau _ { n } ) \tau _ { n } ( \tau _ { n } ) \tau _ { n } ( \tau _ { n } ) \tau _ { n } ( \tau _ { n } ) \tau _ { n } ( \tau _ { n } ) \tau _ { n } ( \tau _ { n } ) \tau _ { n } ( \tau _ { n } ) \tau _ { n } ( \tau _ { n } ) \tau _ { n } ( \tau _ { n } ) \tau _ { n } ( \tau _ { n } ( \tau _ { n } ) \tau _ { n } ( \tau _ { n } ) \tau _ { n } ( \tau _ { n } ) \tau _ { n } ( \tau _ { n } ) \tau _ { n } ( \tau _ { n } ) \tau _ { n } ( \tau _ { n } ) \tau _ { n } ( \tau _ { n } ) \tau _ { n } ( \tau _ { n } ) \tau _ { n } ( \tau _ { n } ) \tau _ { n } ( \tau _ { n } ) \tau _ { n } ( \tau _ { n } ) \tau _ { n } ( \tau _ { n } ) \tau _ { n } ).$$

In (9), the only term depending on the boxes I and J (besides x and y themselves) are c I and c J , the center of the boxes I and J . Hence, in the low-rank approximation, the x -dependent function K x ( -c I , y I 0 l ′ ) no longer depends on J and the y -dependent function K x ( I 0 m , y -c J ) no longer depends on I . As a consequence, the non-trivial manipulation leading to (9) are required for the use of EIM interpolation in a multilevel FMM context.

Notice that in all the kernel evaluations in (9), the first and second variables are always separated by a distance of at least 2 l κ .

## 3.2 A monolevel multipole method

Consider a set of N target points { ¯ x i } 1 ≤ ≤ i N , N source points { ¯ y j } 1 ≤ ≤ j N in Ω, and N potentials { σ j } 1 ≤ ≤ j N , with N glyph[greatermuch] 1. We look for a low complexity approximation of

$$f ( \bar { x } _ { i } ) = \sum _ { j = 1 } ^ { N } \sigma _ { j } K ( \bar { x } _ { i }, \bar { y } _ { j } ), \quad 1 \leq i \leq N.$$

We say that two boxes I and J are well-separated if min x,y ∈ × I J | x - | ≥ y 2 l κ . The set of boxes wellseparated from I is denoted F ( I ), its complement in Ω is the neighborhood of I and is denoted N ( I ). The FMM aims at compressing the 'far' interactions in the following fashion:

$$f _ { \mathcal { F } } ( \bar { x } _ { i } ) = \sum _ { J \in \mathcal { F } ( I ) } \sum _ { \bar { y } _ { j } \in J } \sigma _ { j } K ( \bar { x } _ { i }, \bar { y } _ { j } ), \quad \text{ for all box } I \text{ and all } \bar { x } _ { i } \in I.$$

A monolevel multipole method can be directly derived from (9) to approximate f F (¯ ), see Algox i rithm 2. The obtained approximation of f F (¯ ) is denoted x i f F , FMM (¯ ). x i An approximation of f (¯ ) is x i then given by

$$f ( \bar { x } _ { i } ) \approx f _ { \mathcal { F }, \text{FRM} } ( \bar { x } _ { i } ) + f _ { \mathcal { N } } ( \bar { x } _ { i } ),$$

Table 1: Complexity of each step of Algorithm 2 and of the near interactions term.

| Step number   | Complexity     |
|---------------|----------------|
| 1             | dN             |
| 2             | d 2 (2 D ) κ   |
| 3             | d 2 (2 D ) 2 κ |
| 4             | d 2 (2 D ) κ   |
| 5             | dN             |
| f N           | N 2 (2 D ) - κ |

where the 'near' interactions are defined by

$$f _ { \mathcal { N } } ( \bar { x } _ { i } ) = \sum _ { \bar { y } _ { j } \in \mathcal { N } ( I ) } \sigma _ { j } K ( \bar { x } _ { i }, \bar { y } _ { j } ).$$

## Algorithm 2 Monolevel EIFMM

- 1. Compute W I m := ∑ ¯ y j ∈ I σ K x j ( I 0 m , y ¯ j -c I ), for all box I and all 1 ≤ m ≤ d
- 2. Compute ˆ W I l := ∑ d m =1 ∆ I 0 ,J 0 l,m W I m , for all box I and all 1 ≤ l ≤ d
- 3. Compute l I m ′ := ∑ J ∈F ( I ) ∑ d l =1 K x ( J 0 m ′ , y J 0 l + c J -c I ) ˆ W J l , for all box I and all 1 ≤ m ′ ≤ d
- 4. Compute ˆ l I l ′ := ∑ d m ′ =1 ∆ J ,I 0 0 l ′ ,m ′ l I l ′ , for all box I and all 1 ≤ m ′ ≤ d
- 5. Compute f F , FMM (¯ ) := x i ∑ d l ′ =1 K x (¯ i -c I , y I 0 l ′ ) ˆ l I l ′ , for all box I and all x i ∈ I

In Table 1 is given the complexity of each step of Algorithm 2. For step 3 and the computation of the near interactions to be in the same order in N , one must impose (2 D ) 2 κ ∝ N 2 (2 D ) -κ , leading to a overall complexity of the sum of order N 4 3 . As classically done in the literature, the derivation of a multilevel multipole method allows to lower the complexity in N of the algorithm to linear.

## 3.3 A multilevel multipole method

The domain of computation Ω = ( -L 2 , L 2 ) is now organized as a binary tree with κ +1 levels. Level 0 is the root of the tree: there is one box equal to Ω, level 1 contains the partitioning of Ω into two boxes of length L 2 , and the following levels are obtained by partitioning each box of the previous level in the same fashion, up to level κ .

A superscript k is now added to the quantities previously defined, to indicate that they depend on the level k . Two EIM's are to be considered at each level, in the same fashion as in Section 3.1, where now the domains depend on the level k in the tree: I k 0 = ( -L + l k , -3 l k ) ∪ (3 l k , L -l k ), and J k 0 = ( -l k , l k ).

For symmetric kernels, this amount of precomputation can be lowered to one EIM by level. The key to a multilevel method stands in the ability to derive recursion formulae to propagate information from the level κ to the level 0 in the upward pass, and from the level 0 to the level κ in the downward pass.

Black-Box FMM, at least in its original version, is proposed with a constant-order of interpolation, which enables to propagate the information through the levels of the tree without introducing additional errors. Considering a variable-order Chebyshev interpolation is possible, but one must carefully analyse the introduced error, as done in [6] in the context of H2-matrices. In our method, the domains of approximation depicted in Figure 1 naturally enable a variable-order interpolation depending on the level of the tree, for which the introduced error is the same as the EIM approximation error in (9). As a consequence, the variable-order aspect of the approximation is directly obtained at the required accuracy, and the length of the approximation is automatically optimized at each level of the tree, which leads to important gains for inhomogeneous kernels.

## 3.3.1 Recursion for the upward pass: multipole to multipole (M2M)

The goal is to compute

$$W _ { m } ^ { I ^ { k } } \coloneqq \sum _ { \bar { y } _ { j } \in I ^ { k } } \sigma _ { j } K ( x _ { m } ^ { I ^ { k } _ { 0 } }, \bar { y } _ { j } - c _ { I ^ { k } } ),$$

for all box I k at all level k with 0 ≤ k ≤ κ and all 1 ≤ m ≤ d k by recursion. The initialization of the recursion consists in computing W I κ m at level κ using (14). Then, suppose W I k m is available at level k , there holds

$$W _ { m } ^ { I ^ { k - 1 } } = \sum _ { \vec { y } _ { j } \in I ^ { k - 1 } } \sigma _ { j } K ( x _ { m } ^ { I _ { 0 } ^ { k - 1 } }, \vec { y } _ { j } - c _ { I ^ { k - 1 } } ) = \sum _ { J ^ { k } \in \mathcal { C } ( I ^ { k - 1 } ) } \sum _ { \vec { y } _ { j } \in J ^ { k } } \sigma _ { j } K ( x _ { m } ^ { I _ { 0 } ^ { k - 1 } } + c _ { J ^ { k - 1 } - c _ { J ^ { k } }, \vec { y } _ { j } - c _ { J ^ { k } } ), \quad ( 1 5 )$$

where C ( I k -1 ) denotes the children boxes of I k -1 . To obtain a recursion on W I k -1 m , i.e. to compute W I k -1 m using some values of W J k m , for J k being boxes at the level k , we will use the EIM approximations already computed at level k to derive (9). One must check that the corresponding EIM approximations are evaluated on the right domains, namely I k 0 and J k 0 .

For all ¯ y j ∈ J k , -l k ≤ ¯ y j -c J k ≤ l k , i.e. ¯ y j -c J k ∈ J k 0 . Then, for all J k ∈ C ( I k -1 ), -l k ≤ c I k -1 -c J k ≤ l k (see for instance Figure 3). Moreover, since x I k -1 0 m ∈ I k -1 0 = ( -L + l k -1 , -3 l k -1 ) ∪ (3 l k -1 , L -l k -1 ), it is inferred that

$$- L + l _ { k } \leq x _ { m } ^ { J _ { 0 } ^ { k - 1 } } + c _ { I ^ { k - 1 } } - c _ { J ^ { k } } \leq - 5 l _ { k } \leq - 3 l _ { k } \text{ or } 3 l _ { k } \leq 5 l _ { k } \leq x _ { m } ^ { J _ { 0 } ^ { k - 1 } } + c _ { I ^ { k - 1 } } - c _ { J ^ { k } } \leq L - l _ { k },$$

i.e. x I k -1 0 m + c I k -1 -c J k ∈ I k 0 . Therefore, the EIM from level k can then be used to write the following approximation:

$$K ( x _ { m } ^ { J ^ { k - 1 } _ { 0 } } + c _ { I ^ { k - 1 } } - c _ { J ^ { k } }, \bar { y } _ { j } - c _ { J ^ { k } } ) \approx \sum _ { p, q = 1 } ^ { d ^ { k } } \Delta _ { p, q } ^ { I ^ { k } _ { 0 }, J ^ { k } _ { 0 } } K ( x _ { m } ^ { J ^ { k - 1 } _ { 0 } } + c _ { I ^ { k - 1 } } - c _ { J ^ { k } }, y _ { p } ^ { J ^ { k } _ { 0 } } ) K ( x _ { q } ^ { J ^ { k } _ { 0 } }, \bar { y } _ { j } - c _ { J ^ { k } } ), \quad ( 1 6 )$$

where the approximation error is of the same order as in the (9).

Then, injecting (16) into (15), it is inferred that

$$\int n, \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \end{cases} \\ W ^ { I ^ { k - 1 } _ { m } } _ { m } & \approx \sum _ { p, q = 1 } ^ { d ^ { k } } \Delta ^ { I ^ { k } _ { 0 }, J ^ { k } _ { 0 } } _ { p, q } \sum _ { J ^ { k } \in \mathcal { C } ( I ^ { k - 1 } ) } K ( x ^ { I ^ { k - 1 } _ { 0 } } _ { m } + c _ { I ^ { k - 1 } } - c _ { J ^ { k } }, y ^ { J ^ { k } _ { p } } _ { p } ) \sum _ { \tilde { \nu } _ { j } \in J ^ { k } } \sigma _ { j } K ( x ^ { I ^ { k - 1 } _ { 0 } } _ { m } + c _ { I ^ { k - 1 } } - c _ { J ^ { k } }, y ^ { J ^ { k } _ { p } } _ { p } ) W ^ { J ^ { k } } _ { q }, \\ \text{which would be an exception for the normal use in the same fashion as for the normal}$$

which provides a recursive formula for the upward pass. In the same fashion as for the monolevel version, we then compute

$$\hat { W } _ { l } ^ { I ^ { k } } \coloneqq \sum _ { m = 1 } ^ { d ^ { k } } \Delta _ { l, m } ^ { I ^ { k } _ { 0 }, J ^ { k } _ { 0 } } W _ { m } ^ { I ^ { k } }, \\.$$

for all box I k at all level k with 0 ≤ k ≤ κ and all 1 ≤ l ≤ d k .

## 3.3.2 Transfer pass: multipole to local (M2L)

We denote I ( I k ), the set of boxes that are the children of the boxes in the neighborhood of the parent of I k but are well-separated from I k at level k , and call it the interaction list of I k , see Figure 2.

The M2L step consists in computing

$$g _ { m ^ { \prime } } ^ { I ^ { k } } \coloneqq \sum _ { J ^ { k } \in \mathcal { I } ( I ^ { k } ) } \sum _ { l = 1 } ^ { d ^ { k } } K ( x _ { m ^ { \prime } } ^ { J ^ { k } _ { 0 } }, y _ { l } ^ { J ^ { k } _ { 0 } } + c _ { J ^ { k } } - c _ { I ^ { k } } ) \hat { W } _ { l } ^ { J ^ { k } },$$

for all box I k at all level k with 0 ≤ k ≤ κ and all 1 ≤ m ′ ≤ d k .

Figure 2: Representation of the neighborhood and the interaction list of a box I k and level k .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0210v3-figure-71465c66-90d5-4016-954d-d042a322b094.jpg)

## 3.3.3 Recursion for the downward pass: local to local (L2L)

The goal is to compute

$$l _ { m ^ { \prime } } ^ { I ^ { k } } \coloneqq \sum _ { J ^ { k } \in \mathcal { F } ( I ^ { k } ) } \sum _ { l = 1 } ^ { d ^ { k } } K ( x _ { m ^ { \prime } } ^ { J _ { 0 } ^ { k } }, y _ { l } ^ { J _ { 0 } ^ { k } } + c _ { J ^ { k } } - c _ { I ^ { k } } ) \hat { W } _ { l } ^ { J ^ { k } },$$

for all box I k at all level k with 0 ≤ k ≤ κ and all 1 ≤ m ′ ≤ d k by recursion. The initialization of the recursion consists in computing l I 2 m ′ at level 2 using (20). Then, suppose l I k m ′ is available, there holds

$$l ^ { I ^ { k + 1 } } _ { m ^ { \prime } } = g ^ { I ^ { k + 1 } } _ { m ^ { \prime } } + \sum _ { J ^ { k } \in \mathcal { F } ( \mathcal { P } ( I ^ { k + 1 } ) ) \ J ^ { k + 1 } \in \mathcal { C } ( J ^ { k } ) } \sum _ { l = 1 } ^ { d ^ { k + 1 } } K ( x ^ { J ^ { k + 1 } } _ { m ^ { \prime } }, y ^ { J ^ { k + 1 } } _ { l } + c _ { J ^ { k + 1 } } - c _ { I ^ { k + 1 } } ) \hat { W } ^ { J ^ { k + 1 } } _ { l }$$

$$l _ { m ^ { \prime } } ^ { I ^ { k + 1 } } & = g _ { m ^ { \prime } } ^ { I ^ { k + 1 } } + \sum _ { J ^ { k } \in \mathcal { F } ( \mathcal { P } ( I ^ { k + 1 } ) ) } \sum _ { J ^ { k + 1 } \in \mathcal { C } ( J ^ { k } ) } \sum _ { l = 1 } ^ { d ^ { k + 1 } } K ( x _ { m ^ { \prime } } ^ { J ^ { k + 1 } _ { 0 } }, y _ { l } ^ { J ^ { k + 1 } _ { 0 } } + c _ { J ^ { k + 1 } } - c _ { J ^ { k + 1 } } ) \hat { W } _ { l } ^ { J ^ { k + 1 } _ { 0 } } \\ & = g _ { m ^ { \prime } } ^ { I ^ { k + 1 } _ { 0 } } + \sum _ { J ^ { k } \in \mathcal { F } ( \mathcal { P } ( I ^ { k + 1 } ) ) } \sum _ { J ^ { k + 1 } \in \mathcal { C } ( J ^ { k } ) } \sum _ { l = 1 } ^ { d ^ { k + 1 } _ { 0 } } K ( x _ { m ^ { \prime } } ^ { J ^ { k + 1 } _ { 0 } } + c _ { J ^ { k + 1 } } - c _ { \mathcal { P } ( I ^ { k + 1 } ) }, y _ { l } ^ { J ^ { k + 1 } _ { 0 } } + c _ { J ^ { k + 1 } } - c _ { \mathcal { P } ( I ^ { k + 1 } ) } ) \hat { W } _ { l } ^ { J ^ { k + 1 } _ { 0 } }, \\ \text{where $\mathcal{P} ( I ^ { k + 1 } ) \text{ denotes the parent of $I^{k+1}$}.$$

where P ( I k +1 ) denotes the parent of I k +1 .

It is seen on Figures 3 and 4 that

$$- l _ { k + 1 } \leq c _ { I ^ { k + 1 } } - c _ { \mathcal { P } ( I ^ { k + 1 } ) } \leq l _ { k + 1 },$$

and that

$$- L + 3 l _ { k + 1 } \leq c _ { J ^ { k + 1 } } - c _ { \mathcal { P } ( I ^ { k + 1 } ) } \leq - 7 l _ { k + 1 } \text{ or } 7 l _ { k + 1 } \leq c _ { J ^ { k + 1 } } - c _ { \mathcal { P } ( I ^ { k + 1 } ) } \leq L - 3 l _ { k + 1 }.$$

Figure 3: Possible values for c I k +1 -c P ( I k +1 ) and c J k +1 -c J k , and smallest possible values for c J k +1 -c P ( I k +1 ) and c P ( I k +1 ) -c J k , with J k well-separated from P ( I k +1 ) at level k , and J k +1 child of J k .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0210v3-figure-9586308c-f2dd-452e-b0b7-c9c01714524d.jpg)

$$\text{Since $x_{m^{\prime}}_{0}^{k+1},y_{f}^{k+1}^{k+1} \in J_{0}^{k+1}$, it is then readily verified that $x_{m^{\prime}}_{0}^{k+1}+c_{f}_{f}+1}-c_{\mathcal{P}(f^{k+1}) \in J_{0}^{k}$ and}$$

$$- L + 2 l _ { k + 1 } = - L + l _ { k } \leq y _ { l } ^ { j ^ { k + 1 } _ { 0 } } + c _ { J ^ { k + 1 } } - c _ { \mathcal { P } ( l ^ { k + 1 } ) } \leq - 6 l _ { k + 1 } = - 3 l _ { k } \text{ or } 3 l _ { k } \leq y _ { l } ^ { j ^ { k + 1 } _ { 0 } } + c _ { J ^ { k + 1 } } - c _ { \mathcal { P } ( l ^ { k + 1 } ) } \leq L - l _ { k },$$

Figure 4: Largest possible value for c J k +1 -c P ( I k +1 ) and c P ( I k +1 ) -c J k , with J k well-separated from P ( I k +1 ) at level k , and J k +1 child of J k .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0210v3-figure-9930312a-d2bc-40cc-8c80-d61e6b5c3da4.jpg)

i.e. y J k +1 0 l + c J k +1 -c P ( I k +1 ) ∈ I k 0 . Therefore, the EIM from level k can then be used to write the following approximation

$$K ( x _ { m ^ { \prime } } ^ { J _ { 0 } ^ { k + 1 } } + c _ { I ^ { k + 1 } } - c _ { \mathcal { P } ( I ^ { k + 1 } ) }, y _ { l } ^ { J _ { 0 } ^ { k + 1 } } + c _ { J ^ { k + 1 } } - c _ { \mathcal { P } ( I ^ { k + 1 } ) } ) \approx \\ \sum _ { p ^ { \prime }, q ^ { \prime } = 1 } ^ { d ^ { k } } \Delta _ { p ^ { \prime }, q ^ { \prime } } ^ { J _ { 0 } ^ { k }, J _ { 0 } ^ { k } } K ( x _ { m ^ { \prime } } ^ { J _ { 0 } ^ { k + 1 } } + c _ { I ^ { k + 1 } } - c _ { \mathcal { P } ( I ^ { k + 1 } ) }, y _ { p ^ { \prime } } ^ { I _ { 0 } ^ { k } } ) K ( x _ { q ^ { \prime } } ^ { J _ { 0 } ^ { k } }, y _ { l } ^ { J _ { 0 } ^ { k + 1 } } + c _ { J ^ { k + 1 } } - c _ { \mathcal { P } ( I ^ { k + 1 } ) } ). \\ \sum _ { p ^ { \prime }, q ^ { \prime } = 1 } ^ { d ^ { k } } \cdots \cdots \cdots \cdots \cdots \cdots.$$

Then, injecting (22) into (21), it is inferred that

$$& \L ^ { \cdot } _ { m ^ { \prime } } \approx g ^ { l k + 1 } _ { m ^ { \prime } } + \sum _ { J ^ { k } \in \mathcal { F } ( \mathcal { P } ( I ^ { k + 1 } ) ) } \sum _ { J ^ { k + 1 } \in \mathcal { C } ( J ^ { k } ) } \sum _ { p ^ { \prime }, q ^ { \prime } = 1 } ^ { d ^ { k } } \Delta ^ { J ^ { k } _ { 0 }, I ^ { k } _ { 0 } } _ { p ^ { \prime }, q ^ { \prime } } \sum _ { l = 1 } ^ { d ^ { k + 1 } } K ( x ^ { J ^ { k + 1 } _ { 0 } } _ { m ^ { \prime } } + c _ { I ^ { k + 1 } } - c _ { P ( I ^ { k + 1 } ) }, y ^ { J ^ { k + 1 } _ { 0 } } _ { p ^ { \prime } } ) \\ & \quad K ( x ^ { J ^ { k } _ { 0 } } _ { q ^ { \prime } }, y ^ { J ^ { k + 1 } _ { 0 } } _ { l } + c _ { J ^ { k + 1 } } - c _ { P ( I ^ { k + 1 } ) } ) \hat { W } ^ { J ^ { k + 1 } _ { l } } _ { l }. \\ & \quad \nu ^ { k + 1 } _ { \nu ^ { k + 1 } }$$

Consider K x ( J k 0 q ′ , y J k +1 0 l + c J k +1 -c P ( I k +1 ) ) = K x ( J k 0 q ′ + c P ( I k +1 ) -c J k , y J k +1 0 l + c J k +1 -c J k ). Since J k ∈ F P ( ( I k +1 )) and J k is the parent of J k +1 , it is seen on Figures 3 and 4 that

$$- l _ { k + 1 } \leq c _ { J ^ { k + 1 } } - c _ { J ^ { k } } \leq l _ { k + 1 },$$

and that

$$- L + 4 l _ { k + 1 } \leq c _ { \mathcal { P } ( I ^ { k + 1 } ) } - c _ { J ^ { k } } \leq - 8 l _ { k + 1 } \text{ or } 8 l _ { k + 1 } \leq c _ { \mathcal { P } ( I ^ { k + 1 } ) } - c _ { J ^ { k } } \leq L - 4 l _ { k + 1 }.$$

Since x J k 0 q ′ , y J k +1 0 l ∈ J k +1 0 , J k ∈ F P ( ( I k +1 )), and J k is the parent of J k +1 , it is readily verified that y J k +1 0 l + c J k +1 -c J k ∈ J k 0 and

$$- L + l _ { k } \leq - L + 3 l _ { k + 1 } \leq x _ { q ^ { \prime } } ^ { J ^ { k } _ { 0 } } + c _ { P ( l ^ { k + 1 } ) } - c _ { J ^ { k } } \leq - 7 l _ { k + 1 } \leq - 3 l _ { k } \text{ or } 3 l _ { k } \leq x _ { q ^ { \prime } } ^ { J ^ { k } _ { 0 } } + c _ { P ( l ^ { k + 1 } ) } - c _ { J ^ { k } } \leq L - l _ { k },$$

i.e. x J k 0 q ′ + c P ( I k +1 ) -c J k ∈ I k 0 . Therefore, the EIM from level k can then be used again to write the following approximation

$$\dots \dots \circ \lambda r r \dots \dots \dots \\ K ( x _ { q ^ { \prime } } ^ { J _ { 0 } ^ { k } } + c _ { \mathcal { P } ( I ^ { k + 1 } ) } - c _ { J ^ { k } }, y _ { l } ^ { J _ { 0 } ^ { k + 1 } } + c _ { J ^ { k + 1 } } - c _ { J ^ { k } } ) & \approx \\ \sum _ { u, v = 1 } ^ { d ^ { k } } \Delta _ { u, v } ^ { I _ { 0 } ^ { k }, J _ { 0 } ^ { k } } K ( x _ { q ^ { \prime } } ^ { J _ { 0 } ^ { k } } + c _ { \mathcal { P } ( I ^ { k + 1 } ) } - c _ { J ^ { k } }, y _ { u ^ { \prime } } ^ { J _ { 0 } ^ { k } } ) K ( x _ { v ^ { \prime } } ^ { J _ { 0 } ^ { k } }, y _ { l } ^ { J _ { 0 } ^ { k + 1 } } + c _ { J ^ { k + 1 } } - c _ { J ^ { k } } ). \\ \lambda \dots,$$

Injection (24) into (23), making use of the formula (18) to express ˆ W J k +1 l , it is inferred that

$$l _ { m ^ { \prime } } ^ { l ^ { k + 1 } } & \approx g _ { m ^ { \prime } } ^ { l ^ { k + 1 } } + \sum _ { J ^ { k } \in \mathcal { F } ( \mathcal { P } ( l ^ { k + 1 } ) ) } \sum _ { J ^ { k + 1 } \in \mathcal { C } ( J ^ { k } ) } \sum _ { p ^ { \prime }, q ^ { \prime } = 1 } ^ { d ^ { k } } \Delta _ { p ^ { \prime }, q ^ { \prime } } ^ { J ^ { k } _ { 0 }, l ^ { k } _ { 0 } } \sum _ { l = 1 } ^ { d ^ { k + 1 } } K ( x _ { m ^ { \prime } } ^ { J ^ { k + 1 } _ { 0 } } + c _ { l ^ { k + 1 } } - c _ { \mathcal { P } ( l ^ { k + 1 } ) }, y _ { p ^ { \prime } } ^ { l ^ { k } _ { 0 } } ) \\ & \stackrel { d ^ { k } } { \stackrel { d ^ { k } } { \longrightarrow } } \cdots I ^ { k } \cdots \cdots I ^ { k } \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \dots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots i ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k +1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \dots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ {k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \ddots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots I ^ {k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k - 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdets \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots I ^ { k + 1 + 1 } \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots I ^ {k + 1 } \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cdots I ^ { k + 1 } \cd$$

d

d

∑

I

k

k

k

k

k

k

+1

k

+1

k

+1

∆

0

,J

J

J

I

J

I

,J

k

+1

u,v

0

∑

J

q

′

P

(

I

)

J

0

l

J

J

0

l,m

m

k

+1

k

k

+1

k

u

v

0

K x

(

+

c

-

c

, y

0

)

K x

(

0

, y

+

c

-

c

)

∆

0

W

.

u,v

=1

m

=1

Making use of K x ( I k 0 v , y J k +1 0 l + c J k +1 -c J k ) = K x ( I k 0 v + c J k -c J k +1 , y J k +1 0 l ) and K x ( J k 0 q ′ + c P ( I k +1 ) c J k , y J k 0 u ) = K x ( J k 0 q ′ , y J k 0 u + c J k -c P ( I k +1 ) ), and reorganizing some terms, there holds

$$-$$

$$l _ { m ^ { \prime } } ^ { I ^ { k + 1 } } \approx g _ { m ^ { \prime } } ^ { I ^ { k + 1 } } + \sum _ { p ^ { \prime }, q ^ { \prime } = 1 } ^ { d ^ { k } } \Delta _ { p ^ { \prime }, q ^ { \prime } } ^ { J _ { 0 } ^ { k }, I _ { 0 } ^ { k } } K ( x _ { m ^ { \prime } } ^ { J _ { 0 } ^ { k + 1 } } + c _ { I ^ { k + 1 } } - c _ { \mathcal { P } ( I ^ { k + 1 } ) }, y _ { p ^ { \prime } } ^ { I _ { 0 } ^ { k } } )$$

$$\sum _ { j ^ { k } \in F ( P ( l ^ { k + 1 } ) ) \, u = 1 } ^ { d ^ { k } } K ( x _ { q ^ { \prime } } ^ { j ^ { k } }, y _ { u ^ { \prime } } ^ { k ^ { \prime } } + c _ { j ^ { k } } - c _ { P ( l ^ { k + 1 } ) } ) \sum _ { v = 1 } ^ { d ^ { k } } \Delta _ { u, v ^ { \prime } } ^ { t ^ { k } _ { 0 }, j ^ { k } _ { 0 } } ^ { d ^ { k } _ { 0 } } \sum _ { l, m = 1 } ^ { d ^ { k + 1 } } \Delta _ { l, m } ^ { t ^ { k + 1 } _ { 0 }, j ^ { k + 1 } _ { 0 } } \sum _ { j ^ { k + 1 } \in C ( j ^ { k } ) } K ( x _ { v ^ { \prime } } ^ { j ^ { k } _ { 0 } } + c _ { j ^ { k } } - c _ { j ^ { k + 1 }, j ^ { k } _ { 1 } } ) W _ { m } ^ { t ^ { k + 1 } },$$

︸

︷︷

︸

=

W

J k

︸

= ˆ

J k

v

using (17)

︷︷

using (18)

︸

︸

︸

=

l

P

(

q

′

︷︷

using (20)

(26)

which provides a recursive formula for the upward pass. In the same fashion as for the monolevel version, we then compute

$$\hat { l } ^ { I ^ { \kappa } } _ { l ^ { \prime } } \coloneqq \sum _ { m ^ { \prime } = 1 } ^ { d ^ { \kappa } } \Delta ^ { J ^ { \kappa } _ { 0 }, I ^ { \kappa } _ { 0 } } _ { l ^ { \prime }, m ^ { \prime } } l ^ { I ^ { \kappa } } _ { l ^ { \prime } },$$

for all box I κ at level κ and all 1 ≤ l ′ ≤ d κ .

## 3.3.4 Algorithm and complexity

Using the recursion formulae derived above, the multilevel procedure is detailed in Algorithm 3.

## Algorithm 3 Multilevel EIFMM

- 1. Compute W I κ m = ∑ ¯ y j ∈ I κ σ K x j ( I κ 0 m , y ¯ j -c I κ ), for all box I κ at level κ and all 1 ≤ m ≤ d κ
- 2. Compute W I k m = ∑ d k +1 p,q =1 ∆ I k +1 0 ,J k +1 0 p,q ∑ J k +1 ∈C ( I k ) K x ( I k 0 m + c I k -c J k +1 , y J k +1 0 p ) W J k +1 q , for all box I k at all levels κ -1 ≥ k ≥ 0 and all 1 ≤ m ≤ d k +1
- 3. Compute ˆ W I k l = ∑ d k m =1 ∆ I k 0 ,J k 0 l,m W I k m , for all box I k at all level k with 0 ≤ k ≤ κ and all 1 ≤ l ≤ d k
- 4. Compute g I k m ′ = ∑ J k ∈I ( I k ) ∑ d k l =1 K x ( J k 0 m ′ , y J k 0 l + c J k -c I k ) ˆ W J k l , for all box I k at all level k with 0 ≤ k ≤ κ and all 1 ≤ m ′ ≤ d k .
- 5. Let l I 0 m ′ = g I 0 m ′ and compute l I k m ′ = g I k m ′ + ∑ d k -1 p ,q ′ ′ =1 ∆ J k -1 0 ,I k -1 0 p ,q ′ ′ K x ( J k 0 m ′ + c I k -c P ( I k ) , y I k -1 0 p ′ ) l P ( I k ) q ′ for all box I k at all levels 1 ≤ k ≤ κ and all 1 ≤ m ′ ≤ d k -1 κ κ
- 7. Compute f F , FMM (¯ ) = x i ∑ d κ l ′ =1 K x (¯ i -c I κ , y I κ 0 l ′ ) ˆ l I κ l ′ , for all box I κ and all x i ∈ I κ at level κ
- 6. Compute ˆ l I κ l ′ := ∑ d κ m ′ =1 ∆ J 0 ,I 0 l ′ ,m ′ l I κ l ′ for all box I κ at level κ and all 1 ≤ l ′ ≤ d κ

The complexity of the multilevel procedure is given in Table 2. For step 4 and the computation of the near interactions to be in the same order in N , one must impose (2 D κ ) ∝ N 2 (2 D ) -κ , leading to a overall complexity of the sum of order N . The choice that minimizes the overall cost is κ ∝ log( N ) and a constant number of points per box at the level κ . As for any multilevel FMM algorithm, the complexity of step 3 in Algorithm 2 has been reduced by replacing the sum over J ∈ L ( I ) (which is of complexity linear in N ) by a sum over J ∈ I ( I ) (which is of complexity independent in N ), making use of the fact that the set of boxes well-separated from I k at level k contains the set of boxes well-separated from the parent of I k at level k -1. By precomputing some terms, independent to the location of the particles, the complexity of the steps 2, 4 and 5 can be reduced.

## 4 Optimizing the overall complexity

## 4.1 Precomputations for the M2M and L2L steps

Recall that when using an EIM approximation, there holds ∆ = B -t Γ -1 , where the matrices B and Γ are triangular and constructed in a precomputation step. In practice, the condition number of the

I k

+1 )

W

u

Table 2: Complexity of each step of Algorithm 3 (for the most expensive level) and of the near interactions term ; n I denotes the maximum number of elements in the interaction list of a box (4 in 1D, 40 in 2D and 316 in 3D).

| Step number   | Complexity       |
|---------------|------------------|
| 1             | dN               |
| 2             | d 3 (2 D ) κ     |
| 3             | d 2 (2 D ) κ     |
| 4             | d 2 n I (2 D ) κ |
| 5             | d 3 (2 D ) κ     |
| 6             | d 2 (2 D ) κ     |
| 7             | dN               |
| f N           | N 2 (2 D ) - κ   |

matrix Γ worsens as the EIM approximation gets more accurate. It is well known in linear algebra that when a matrix Γ is ill-conditioned, a direct resolution of the linear system Γ x = b for some b ∈ R d provides a result much less polluted by numerical errors than first inverting the matrix Γ and then computing the matrix vector product Γ -1 b . Since the matrices Γ and B are triangular, computing ∆ = b y by first solving Γ x = b and then solving B y t = x is of the same complexity d 2 as the direct matrix-vector product.

Consider step 2 in Algorithm 3 at level k :

$$W _ { m } ^ { I ^ { k } } = \sum _ { p, q = 1 } ^ { d ^ { k + 1 } } \Delta _ { p, q } ^ { I _ { 0 } ^ { k + 1 }, J _ { 0 } ^ { k + 1 } } \sum _ { J ^ { k + 1 } \in \mathcal { C } ( I ^ { k } ) } K ( x _ { m } ^ { I _ { 0 } ^ { k } } + c _ { I ^ { k } } - c _ { J ^ { k + 1 } }, y _ { p } ^ { J _ { 0 } ^ { k + 1 } } ) W _ { q } ^ { J ^ { k + 1 } }.$$

Wesee that W I k m = ∑ J k +1 ∈C ( I k ) K I k ,J k +1 ∆ I k +1 0 ,J k +1 0 W J k +1 , where K I k ,J k +1 m,p = K x ( I k 0 m + c I k -c J k +1 , y J k +1 0 p ), 1 ≤ m ≤ d k , 1 ≤ p ≤ d k +1 . The matrix K I k ,J k +1 M M 2 := K I k ,J k +1 ∆ I k +1 0 ,J k +1 0 can be precomputed as ( (∆ I k +1 0 ,J k +1 0 ) ( t K I k ,J k +1 ) t ) t , where each matrix-vector product for each column of ( K I k ,J k +1 ) t is computed solving linear systems involving B and Γ matrices to preserve numerical accuracy, as explained in the previous paragraph. Notice that K I k ,J k +1 M M 2 only depends on the relative position of I k and J k +1 . Since J k +1 is a child of I k , there are actually only 2 D different operators K I k ,J k +1 M M 2 , where D denotes the dimension of the space containing the particles. We denote them K i M M 2 , 1 ≤ i ≤ 2 D , and step 2 in Algorithm 3 at level k is reduced to

$$W _ { m } ^ { I ^ { k } } = \sum _ { i = 1 } ^ { 2 ^ { D } } K _ { M 2 M } ^ { i } W _ { q } ^ { J ^ { k + 1 } _ { i } },$$

where J k +1 i denotes the i th child of I k . This step is now of complexity d 2 2 D (2 D k ) at level k , and the leading order of step 2 in Algorithm 3 is d 2 2 D (2 D κ ) -1 = d 2 (2 D κ ) .

The same optimization is done for step 5 in Algorithm 3, reducing the complexity of this step to the order d 2 (2 D κ ) .

## 4.2 Optimization of the M2L step

It is well known that the M2L step is the most expensive one. Numerous efforts have been made to reduce the execution time of this step without degrading the accuracy of the overall result, see in particular [25]. We use the optimization denoted SArcmp in [25]. It consists of a bivariate ACA applied to the collection of all the M2L operators at each level, followed by a truncated SVD. The bivariate ACA and truncated SVD is then applied for each M2L operator at each level to increase further the compression. The procedure is detailed for symmetric kernels for completeness of the presentation.

Consider a level k and the collection of all the M2L operators at this level. Consider a level k , a box I k of the tree at this level k , and a box J k in the interaction list of I k . The vector δ := c J k -c I k only depends on the relative position of the boxes I k and J k . Therefore, there are only n I = 7 D -3 D different M2L matrices K δ i,j := K x ( J k 0 i , y J k 0 j + ) per level (we recall that δ D denotes the dimension of the space containing the particles), see (19). We suppose that EIM algorithms have been carriedout for a certain error bound glyph[epsilon1] , leading to the choice of d k interpolation points. The M2L matrices are organized in one large matrix as follows K fat := [ K K ...K δ 1 δ 2 δ n I ], where K fat ∈ R d k × n d I k . Two rectangular matrices U, V are computed by a bivariate ACA such that | UV | F ≤ | glyph[epsilon1] K fat | F , where | · | F denotes the Frobenius norm. Then, two QR decompositions are computed, such that U = Q R U U and ( V ) t = Q R V V , and a SVD is computed on R U ( R V ) . t Denote s i , 1 ≤ i ≤ d k , the singular values of R U ( R V ) t in decreasing order, and ˆ , U ˆ V the unitary matrices such that R U ( R V ) t = ˆ diag( U s V ) ˆ . The two following filters are applied: (i) keep the first j singular values such that s j s 0 ≤ glyph[epsilon1] ≤ s j -1 s 0 , (ii) keep the first j singular values such that ∑ j i =1 s i ∑ d k i =1 s i ≥ 1 -glyph[epsilon1] ≥ ∑ j -1 i =1 s i ∑ d k i =1 s i . We denote r k the number of kept singular values. Denote ˆ U r k and ˆ V r k the first r k columns of respectively ˆ U and ˆ . V Following [13], we can show that the M2L step can be reduced to size r k × r k (instead of size d k × d k ) matrix-vector products involving the matrices C δ := ˆ U T r k K U δ ˆ r k , plus low complexity pre- and post-processing steps. Then, the same compression is applied locally for each compressed M2L operator C δ . This time, the bivariate ACA is applied on each C δ , producing the rectangular matrices U and V . The QR decompositions, SVD and filters are applied in the same fashion, selecting the singular values s k . The rectangular matrices ˆ U s k and ˆ V s k are constructed, and a diagonal matrix S s k , whose diagonal entries are the square roots of the selected singular values, is constructed. Then, we construct ˜ U s k = Q U U ˆ s k S s k ∈ R r k × s k and ˜ V s k = SV ˆ s k Q t V ∈ R s k × r k , so that C δ is approximated by ˜ U s k ˜ V s k . Thus, size r k × r k matrix-vector products involving the matrices C δ are replaced by two r k × s k matrix-vector products involving the matrices ˜ U s k and ˜ V s k .

## 5 Numerical experiments in 3D

## 5.1 Some elements on the implementation

The EIFMM has been implemented in C++, using the open source library scalfmm, see [9]. From [29], scalfmm is one of the fastest opensource library for FMM computations. We will mainly compare our implementation of the EIFMM to the Chebyshev interpolation-based implementation in scalfmm, which follows the Black-Box FMM [13]. In our code, the EIM precomputation and the SArcmp precomputation steps are performed in C++ using BLAS and LAPACK, and stored in binary files. These precomputations are then read by the modified scalfmm library, using also BLAS for the basic linear algebra operations.

## 5.2 Test-cases

We consider sets of points included in Ω = ( -0 5 . , 0 5) . 3 , randomly taken following (i) a uniform law in Ω, (ii) a uniform law on the sphere of radius 0 5 and (iii) a modification of the set (ii) in an ellipsoid, . see Figure 5.

When comparing runtimes, we have to make sure that the Black-Box FMM is used in a optimal way. There are two parameters to choose: the number of interpolation points per dimension and the accuracy glyph[epsilon1] c of the compression. First, we chose 3 interpolation points per dimension and compute the relative error with glyph[epsilon1] c = 10 -15 . Then, if this relative error is smaller than the desired one, we search by dichotomy the largest value of the glyph[epsilon1] c so that the relative error is smaller or equal to the desired one. Otherwise, we increase the number of interpolation points per dimension by one and repeat the process. That way, the runtimes comparison is fair in the sense that Block-Box FMM is not artificially slowed down by a too large number of interpolation points together with a too large value for glyph[epsilon1] c . Notice that this trick is necessary since we do not know in advance how the Chebyshev interpolation will perform in practice for a new kernel, while for the EIFMM, the error estimate automatically selects the right number of interpolation points.

Figure 5: Sets of points considered in the numerical applications, from left to right: random (5000 points), sphere (10100 points), and ellipsoid (1944 points). The values for the number of points have been chosen for representation reasons, higher values will be considered in the numerical experiments.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0210v3-figure-463e99d8-fd2d-4f15-8f10-155af4885477.jpg)

## 5.3 Performance for the Laplace kernel

In this section, we compare runtimes of the FMM summation using the Black-Box FMM and the EIFMM. The computations are carried out in sequential, on a laptop with 4 GB of RAM and a Intel(R) Core(TM) i5 CPU M 560 @ 2.67GHz. For the Black-Box FMM, we use the so-called symmetric M2L compression, based on local SVD for each M2L operator. The symmetric M2L compression exploit the symmetry of the Chebyshev nodes to express the M2L step as matrix-matrix products, and even if the complexity of this step is not improved, runtimes are drastically improved due to better reuse of memory cache by the processor. Proposed in [25], the symmetric compression is currently, to the author's knowledge, the most efficient compression strategy for the Black-Box FMM. Since the approximation (9) used in EIFMM has four summations instead of two for Black-Box FMM, the steps 3 and 6 in Algorithm 3 for EIFMM are not present in Black-Box FMM.

In Figure 6 are represented execution time with respect to the relative error on the FMM summation for the Laplace kernel, for the symmetric Black-Box FMM, the EIFMM and the fastest algorithm available in scalfmm. This latter is dedicated to the Laplace kernel, explaining why it performs better than the other two. The test-case corresponds to 10 6 random points taken in the unit cube, organized in a 6-level tree, see Figure 5, left picture. The execution times are better for Black-Box FMM at low accuracy, and equivalent between EIFMM and Black-Box FMM at high accuracy. Notice that since the Laplace kernel is homogeneous, the difficulty of the approximation problem does not change from one level to another in the tree. As a consequence, the variable-order capability of our algorithm is not exploited for this kernel. In what follows, inhomogeneous kernels are considered and gains on runtimes are observed. Precomputation times will be investigated in the following section as well.

## 5.4 Performance for some inhomogeneous kernels

The advantages of the EIFMM mainly apply to inhomogeneous kernels, for which the error estimation allows the selection of an optimal number of interpolation points at each level of the tree by the EIM.

## 5.4.1 1-million-point test-cases

The three test-cases presented in Figure 5 are considered: 10 6 random points in the unit cube with a 6level tree (denoted 'cube'), 992 046 points at the surface of a sphere and an ellipsoid, with respectively , 8-level and 9-level trees (denoted respectively 'sphere' and 'ellipsoid').

Figures 7, 8 and 9 present execution times of the long-range FMM accelerated terms with respect to the relative error on the FMM summation, respectively for the kernels cos 20 r r , e -r 2 and √ r 2 +1, for the three test-cases (cube, sphere and ellipsoid). The reported runtimes do not contain the times required to build the tree, nor the EIM, M2L compression, M2M and L2L precomputations. For the EIFMM, we use the SArcmp compression for the M2L operators, as described in Section 4.2. Due to smaller interpolation formulae, EIFMM enables memory savings with respect to Black-Box FMM. In the figures, missing data correspond to cases requiring more memory than available with our computer,

Figure 6: Execution time with respect to the relative error on the FMM summation for the Laplace kernel; black: symmetric Black-Box FMM, blue: EIFMM, red: fastest spherical algorithm available in scalfmm.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0210v3-figure-0665b4b8-c7b0-4910-9ecc-a8f5802b9ae9.jpg)

for which EIFMM can be used, but not Blakc-Box FMM. For the three kernels, and especially the last two, the EIFMM is much faster than the Black-Box FMM, because the error estimation and the greedy procedure in EIFMM allows the selection of an optimal number of interpolation points at each level of the tree. For instance, for the kernel e -r 2 on the ellipsoid and a relative error of 10 -6 , the Black-Box FMM needs 6 interpolation points per dimension (i.e. 216 interpolation points), whereas the EIFMM needs respectively 167, 74, 36, 21, 13, 10, 9, and 4 points at levels 2 to 9 of the tree. As we go deeper in the tree, the interpolation problems become easier, and the EIM selects fewer points to produce an interpolation of the same quality. √

Table 3 provides precomputations and executions times for the kernels cos(20 r ) r , e -r 2 and r 2 +1 on the cube test-case. Notice that in order not to artificially disadvantage Black-Box FMM, we report values with EIFMM set to a better relative accuracy and we neglect the precomputation times of BlackBox FMM (which are short anyway). In the reported cases, EIFMM becomes advantageous (i.e. the EIFMM precomputations times are compensated by faster long-range FMM terms computation time) if one has to compute more than 2 to 18 FMM summations.

√

| kernel                             | cos(20 r ) r   | cos(20 r ) r   | e - r 2    | e - r 2    | r 2 +1     | r 2 +1     |
|------------------------------------|----------------|----------------|------------|------------|------------|------------|
| algorithm                          | EIFMM          | BBFMM          | EIFMM      | BBFMM      | EIFMM      | BBFMM      |
| relative accuracy                  | 2 × 10 - 4     | 4 × 10 - 4     | 2 × 10 - 6 | 3 × 10 - 6 | 1 × 10 - 5 | 2 × 10 - 5 |
| EIM precomp time (s)               | 73             | -              | 30         | -          | 10         | -          |
| compression time (s)               | 177            | -              | 33         | -          | 1.1        | -          |
| tree construction time (s)         | 3              | 3              | 3          | 3          | 3          | 3          |
| long-range FMM terms comp time (s) | 52             | 69             | 5.4        | 37         | 3.5        | 4.3        |
| short-range terms comp time(s)     | 32             | 32             | 30         | 30         | 16         | 16         |
| minimum EIFMM calls                | 18             | 18             | 2          | 2          | 14         | 14         |

Table 3: Precomputations and executions times for the kernels cos(20 r ) r , e -r 2 and √ r 2 +1 on the cube test-case. The last line indicates the minimum FMM summations to compute for EIFMM to be faster than Black-Box FMM (BBFMM), taking all precomputations times into account.

## 5.4.2 100-million-point test-case

We consider 10 8 random points in the unit cube (Figure 5, left) in a 9-level tree, with the e -r 2 kernel. The computations are carried-out in sequential on a Intel(R) Xeon(R) CPU E5-4620 v2 @ 2.60GHz

Figure 7: Execution time of the long-range FMM accelerated terms with respect to the relative error on the FMM summation for the kernel cos 20 r r , black: symmetric Black-Box FMM, blue: EIFMM. From top to bottom and left to right: cube, sphere and ellipsoid test-cases.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0210v3-figure-e6f8b68a-59fd-4e63-a022-efa9569aaf9c.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0210v3-figure-37eca037-3645-41a6-bcd6-6d59c084e38f.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0210v3-figure-45d7ac0f-ffc7-4469-8fb7-5908bf37c1f7.jpg)

## and 128 GB of RAM.

Table 4 presents the relative accuracy and the execution times of all the steps for two simulations using EIFMM and Black-Box FMM. For a better accuracy, the EIFMM is faster than the BlackBox FMM, including the precomputation times (increasing the accuracy with Black-Box FMM would require more than 128 GB of RAM). In this case, EIFMM is interesting even for a single FMM computation. Notice that the precomputations for EIFMM are negligible for such large cases: the precomputation complexity is linear with respect to the number of levels in the tree (therefore logarithmic with respect to N ). Moreover, for inhomogeneous kernels, the approximation problem for an additional level is easier than the previous approximation problems, leading to even faster additional EIM and M2L compression precomputation times.

## 6 Conclusion and outlook

This work introduces a new Fast Multipole Method (FMM) using a low-rank approximation of the kernel based on the Empirical Interpolation Method (EIM), called the Empirical Interpolation Fast Multipole Method (EIFMM). The proposed multilevel algorithm is implemented in scalfmm, a FMM library written in C++. The important feature of the EIFMM is a built-in error estimation of the interpolation error made by the low-rank approximation of the far-field behavior of the kernel. As a consequence, the algorithm automatically selects the optimal number of interpolation points required to ensure a given accuracy for the result, leading to important execution time reduction for

Figure 8: Execution time of the long-range FMM accelerated terms with respect to the relative error on the FMM summation for the kernel e -r 2 , black: symmetric Black-Box FMM, blue: EIFMM. From top to bottom and left to right: cube, sphere and ellipsoid test-cases.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0210v3-figure-db145506-9740-406b-b9a7-a3841dc3e38b.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0210v3-figure-78c5a077-7ea3-4727-be7a-b30719d7f491.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0210v3-figure-0901a6b4-b3e3-480f-a53d-0d3f85a7a448.jpg)

inhomogeneous kernels, for which the difficulty of the approximation varies from one level of the tree to another.

## Acknowledgement

The author wish to thank Centre National d'Etudes Spatiales for financial support through the TOSCA committee, Guillaume Sylvand (Airbus Group Innovations, INRIA) for drawing his attention to fast multipole methods, thoroughly reading the manuscript and for many fruitfull discussions, Tony Leli` evre (CERMICS) and Alexandre Ern (CERMICS) for fruitful discussions, and Mathieu Ruffat (Diginext) for his C++ expertise.

## References

- [1] B. Alpert, G. Beylkin, R. Coifman, and V. Rokhlin. Wavelet-like bases for the fast solutions of second-kind integral equations. SIAM J. Sci. Comput. , 14(1):159-184, 1993.
- [2] M. Barrault, Y. Maday, N. C. Nguyen, and A. T. Patera. An 'empirical interpolation' method: application to efficient reduced-basis discretization of partial differential equations. Comptes Rendus Mathematique , 339(9):667 - 672, 2004.

Figure 9: Execution time of the long-range FMM accelerated terms with respect to the relative error on the FMM summation for the kernel √ r 2 +1, black: symmetric Black-Box FMM, blue: EIFMM. From top to bottom and left to right: cube, sphere and ellipsoid test-cases.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0210v3-figure-78caa22f-7855-48f7-a11b-e7510d462bd8.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0210v3-figure-6bcf60c6-6357-4d5e-af0e-e18b3d33407d.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0210v3-figure-3decbe7c-c25b-4cc6-aa1e-8002ef5d3db5.jpg)

- [3] M. Bebendorf. Approximation of boundary element matrices. Numerische Mathematik , 86(4):565589, 2000.
- [4] M. Bebendorf and S. Rjasanow. Adaptive low-rank approximation of collocation matrices. Computing , 70(1):1-24, 2003.
- [5] M. Bebendorf and R. Venn. Constructing nested bases approximations from the entries of nonlocal operators. Numerische Mathematik , 121(4):609-635, 2012.
- [6] S. B¨rm, M. L¨hndorf, and J.M. Melenk. Approximation of integral operators by variable-order o o interpolation. Numerische Mathematik , 99(4):605-643, 2005.
- [7] F. Casenave, A. Ern, and T. Leli`vre. A nonintrusive reduced basis method applied to aeroacoustic e simulations. Advances in Computational Mathematics , pages 1-26, 2014.
- [8] H. Cheng, Z. Gimbutas, P.G. Martinsson, and V. Rokhlin. On the compression of low-rank matrices. SIAM J. Sci. Comput. , 26(4):1389-1404, 2005.
- [9] O. Coulaud, B. Bramas, and C. Piacibello. Scalfmm, C++ Fast Multipole Method Library for HPC. http://scalfmm-public.gforge.inria.fr/doc/ .
- [10] A. Dutt, M. Gu, and V. Rokhlin. Fast algorithms for polynomial interpolation, integration, and differentiation. SIAM Journal on Numerical Analysis , 33(5):1689-1711, 1996.

Table 4: Accuracy and execution times for the 100 million points test-case. P2M, M2M, M2L, L2L and L2P corresponds respectively to the steps (1), (2-3), (4), (5-6) and (7) in Algorithm 3.

| algorithm                  | EIFMM       | Black-Box FMM   |
|----------------------------|-------------|-----------------|
| relative accuracy          | 7 × 10 - 5  | 6 × 10 - 4      |
| EIM precomp time           | 22.8 s      | -               |
| compression time           | 1.3 s       | -               |
| tree construction time (s) | 9 min 10 s  | 9 min 10 s      |
| P2M (1)                    | 1 min 13 s  | 2 min 17 s      |
| M2M (2-3)                  | 31 s        | 49 s            |
| M2L (4)                    | 10 min 39 s | 53 min 01 s     |
| L2L (5-6)                  | 41 s        | 59 s            |
| L2P (7)                    | 1 min 42 s  | 3 min 55 s      |
| total long-range FMM terms | 14 min 46 s | 1 h 01 min 01 s |
| short-range terms          | 19 min 36 s | 19 min 36 s     |
| total                      | 43 min 56 s | 1 h 29 min 47 s |

- [11] A. Dutt and V. Rokhlin. Fast fourier transforms for nonequispaced data. SIAM Journal on Scientific Computing , 14(6):1368-1393, 1993.
- [12] J. L. Eftang and B. Stamm. Parameter multi-domain 'hp' empirical interpolation. International Journal for Numerical Methods in Engineering , 90:412-428, 2012.
- [13] W. Fong and E. Darve. The black-box fast multipole method. Journal of Computational Physics , 228(23):8712 - 8725, 2009.
- [14] K. Giebermann. Multilevel approximation of boundary integral operators. Computing , 67(3):183207, 2001.
- [15] Z. Gimbutas, L. Greengard, and M. Minion. Coulomb interactions on planar structures: Inverting the square root of the laplacian. SIAM Journal on Scientific Computing , 22(6):2093-2108, 2001.
- [16] S.A. Goreinov, E.E. Tyrtyshnikov, and N.L. Zamarashkin. A theory of pseudoskeleton approximations. Linear Algebra and its Applications , 261(1-3):1 - 21, 1997.
- [17] L. Greengard and V. Rokhlin. A fast algorithm for particle simulations. Journal of Computational Physics , 135(2):280-292, 1997.
- [18] W. Hackbusch. A sparse matrix arithmetic based on H-matrices. part I: Introduction to Hmatrices. Computing , 62(2):89-108, 1999.
- [19] W. Hackbusch and S. B¨ orm. Data-sparse approximation by adaptive H2-matrices. Computing , 69(1):1-35, 2002.
- [20] W. Hackbusch and S. B¨ orm. H2-matrix approximation of integral operators by interpolation. Applied Numerical Mathematics , 43(1-2):129 - 143, 2002. 19th Dundee Biennial Conference on Numerical Analysis.
- [21] W. Hackbusch and Z.P. Nowak. On the fast matrix multiplication in the boundary element method by panel clustering. Numerische Mathematik , 54(4):463-491, 1989.
- [22] Y. Maday and O. Mula. A generalized empirical interpolation method: Application of reduced basis techniques to data assimilation. 4:221-235, 2013.
- [23] Y. Maday, N. C. Nguyen, A. T. Patera, and S. Pau. A general multipurpose interpolation procedure: the magic points. Communications On Pure And Applied Analysis , 8(1):383-404, 2008.

| [24]   | P.G. Martinsson and V. Rokhlin. An accelerated kernel-independent fast multipole method in one dimension. SIAM J. Sci. Comput , 26:1389-1404, 2005.                             |
|--------|---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| [25]   | M. Messner, B. Bramas, D. Coulaud, and E. Darve. Optimized M2L Kernels for the Chebyshev Interpolation based Fast Multipole Method. CoRR , abs/1210.7292, 2012.                 |
| [26]   | B. Peherstorfer, D. Butnaru, K. Willcox, and H.-J. Bungartz. Localized discrete empirical inter- polation method. SIAM Journal on Scientific Computing , 36(1):A168-A192, 2014. |
| [27]   | E. Tyrtyshnikov. Mosaic-skeleton approximations. CALCOLO , 33(1-2):47-57, 1996.                                                                                                 |
| [28]   | L. Ying, G. Biros, and D. Zorin. A kernel-independent adaptive fast multipole algorithm in two and three dimensions. J. Comput. Phys. , 196(2):591-626, 2004.                   |
| [29]   | R. Yokota. Fast multipole methods benchmark. https://sites.google.com/site/rioyokota/research/fmm                                                                               |
Fabien Casenave
A Fast Summation Method for translation invariant kernels
We derive a Fast Multipole Method (FMM) where a low-rank approximation of the
kernel is obtained using the Empirical Interpolation Method (EIM). Contrary to
classical interpolation-based FMM, where the interpolation points and basis are
fixed beforehand, the EIM is a nonlinear approximation method which constructs
interpolation points and basis which are adapted to the kernel under
consideration. The basis functions are obtained using evaluations of the kernel
itself. We restrict ourselves to translation-invariant kernels, for which a
modified version of the EIM approximation can be used in a multilevel FMM
context; we call the obtained algorithm Empirical Interpolation Fast Multipole
Method (EIFMM). An important feature of the EIFMM is a built-in error
estimation of the interpolation error made by the low-rank approximation of the
far-field behavior of the kernel: the algorithm selects the optimal number of
interpolation points required to ensure a given accuracy for the result,
leading to important gains for inhomogeneous kernels.