PBBFMM3D: a parallel black-box algorithm for kernel matrix-vector multiplication

TL;DR

This paper presents PBBFMM3D, a parallel black-box algorithm that efficiently computes kernel matrix-vector products with linear complexity, suitable for large-scale scientific computations on multi-core systems.

Contribution

It introduces a novel parallel black-box method for kernel matrix-vector multiplication with proven $O(N)$ complexity and practical implementation for shared-memory architectures.

Findings

Achieves up to 19x speedup on 32 cores

Demonstrates effectiveness in geostatistics PCA applications

Provides a scalable solution for large kernel matrices

Abstract

Kernel matrix-vector product is ubiquitous in many science and engineering applications. However, a naive method requires $O(N^2)$ operations, which becomes prohibitive for large-scale problems. We introduce a parallel method that provably requires $O(N)$ operations to reduce the computation cost. The distinct feature of our method is that it requires only the ability to evaluate the kernel function, offering a black-box interface to users. Our parallel approach targets multi-core shared-memory machines and is implemented using OpenMP. Numerical results demonstrate up to $19\times$ speedup on 32 cores. We also present a real-world application in geostatistics, where our parallel method was used to deliver fast principle component analysis of covariance matrices.

Tables

Table 1. Table 1: Parallel running time of PBBFMM3D . The kernel function is 1 / r 1 𝑟 1/r , and we used Chebyshev interpolation with p = 4 𝑝 4 p=4 . We set the number of tree levels to be 5, 6 and 7 for the three increasing problem sizes, and the errors are 2.10 ​ e- ​ 5 2.10 e- 5 2.10\text{e-}{5} , 2.08 ​ e- ​ 5 2.08 e- 5 2.08\text{e-}{5} , and 2.10 ​ e- ​ 5 2.10 e- 5 2.10\text{e-}{5} , respectively (independent of the number of cores used).

$N$	Time (seconds)
	1 core	2 cores	4 cores	8 cores	16 cores	32 cores
$8^2$e+4	$5.74$e+0	$3.32$e+0	$1.83$e+0	$9.93$e-1	$5.84$e-1	$4.44$e-1
$8^3$e+4	$4.72$e+1	$2.58$e+1	$1.40$e+1	$7.87$e+0	$4.40$e+0	$3.49$e+0
$8^4$e+4	$4.04$e+2	$2.15$e+2	$1.14$e+2	$6.37$e+1	$3.72$e+1	$2.13$e+1

Equations

$${\varphi }_i=\sum _{j=1}^N\mathcal{K}({\mathbf{x}}_i,{\mathbf{y}}_j){\sigma }_j,i=1,\dots ,N$$

$$\mathbf{\varphi }=\mathbf{K}\mathbf{\sigma },$$

$$\mathcal{S}\otimes \mathcal{S}\otimes \mathcal{S}=\left\{{\mathbf{x}}_i^{\ast }=(x_{i_1}^{\ast },x_{i_2}^{\ast },x_{i_3}^{\ast }),i_1,i_2,i_3=1,2,\dots ,p\right\},$$

$${\ell }_p(x_i^{\ast },x)={\Pi }_{1\le k\le p,k\ne i}\frac{x-x_k^{\ast }}{x_i^{\ast }-x_k^{\ast }},i=1,2,\dots ,p.$$

$$L_p({\mathbf{x}}_i^{\ast },\mathbf{x})={\ell }_p(x_{i_1}^{\ast },x_1){\ell }_p(x_{i_2}^{\ast },x_2){\ell }_p(x_{i_3}^{\ast },x_3),i=1,2,\dots ,p^3,$$

$$\mathcal{K}(\mathbf{x},\mathbf{y})\approx$$

$$\approx$$

$$\mathcal{S}=\left\{\cos \left(\frac{2k-1}{2p}\pi \right),k=1,2,\dots ,p\right\}$$

$$\mathcal{S}=\left\{\frac{k-1}{p-1},k=1,2,\dots ,p\right\}.$$

$${\mathbf{M}}_i^{\alpha }=\sum _{{\mathbf{y}}_j\in \alpha }{L}_p({\mathbf{y}}_i^{\ast },{\mathbf{y}}_j){\sigma }_j,i=1,2,\dots ,p^3,$$

$${\mathbf{M}}_i^{\alpha }=\sum _{\beta \in \mathcal{C}(\alpha )}\sum _{j=1}^{p^3}{L}_p({\mathbf{y}}_i^{\ast },{\mathbf{y}}_j^{\ast }){\mathbf{M}}_j^{\beta },i=1,2,\dots ,p^3,$$

$${\mathbf{F}}_i^{\alpha }=\sum _{\beta \in \mathcal{I}(\alpha )}\sum _{j=1}^{p^3}\mathcal{K}({\mathbf{x}}_i^{\ast },{\mathbf{y}}_j^{\ast }){\mathbf{M}}_j^{\beta },i=1,2,\dots ,p^3,$$

$${\mathbf{F}}_i^{\alpha }+=\sum _{j=1}^{p^3}{\mathbf{F}}_j^{\beta }{L}_p({\mathbf{x}}_j^{\ast },{\mathbf{x}}_i^{\ast }),i=1,2,\dots ,p^3,$$

$${\varphi }_i=\sum _{j=1}^{p^3}{\mathbf{F}}_j^{\alpha }{L}_p({\mathbf{x}}_j^{\ast },{\mathbf{x}}_i),$$

$${\varphi }_i+=\sum _{\beta \in \mathcal{N}(\alpha )}\sum _{{\mathbf{y}}_j\in \beta }\mathcal{K}({\mathbf{x}}_i,{\mathbf{y}}_j){\sigma }_j,$$

$$\mathcal{K}(\alpha \mathbf{x},\alpha \mathbf{y})={\alpha }^m\mathcal{K}(\mathbf{x},\mathbf{y}),$$

$$\mathcal{K}(\mathbf{x},\mathbf{y})=\mathcal{K}(\mathbf{y},\mathbf{x}),$$

$$\mathcal{K}(\mathbf{x},\mathbf{y})=-\mathcal{K}(\mathbf{y},\mathbf{x}).$$

$$\frac{\Vert \hat{K}\mathbf{\sigma }-K\mathbf{\sigma }{\Vert }_2}{\Vert K\mathbf{\sigma }{\Vert }_2},$$

Full text

1 Introduction

We consider the problem of computing kernel matrix-vector products, where the kernel function $\mathcal{K}({\bf x},{\bf y})$ is non-oscillatory, translation invariant, and sufficiently smooth (NOTIS) ($\mathcal{K}$ can be singular when ${\bf x}={\bf y}$). For example, a NOTIS function can be $1/\|{\bf x}-{\bf y}\|$ or exp($-\|{\bf x}-{\bf y}\|$).

The problem can be formulated mathematically as evaluating

[TABLE]

where $\{{\bf x}_{i}\}_{i=1}^{N}$ and $\{{\bf y}_{i}\}_{i=1}^{N}$ are the target and the source data points in a cubical domain, respectively, and $\sigma_{j}$ is the weight associated with ${\bf y}_{j}$. In many applications, the two sets of points $\{{\bf x}_{i}\}_{i=1}^{N}$ and $\{{\bf y}_{i}\}_{i=1}^{N}$ may overlap. Algebraically, Eq. 1 can be written as the following matrix-vector multiplication:

[TABLE]

where ${\bm{\sigma}}=[\sigma_{1},\ldots,\sigma_{N}]^{T}$ and $\bm{\phi}=[\phi_{1},\ldots,\phi_{N}]^{T}$ are two vectors, and $\bm{K}$ is an $N$-by-$N$ matrix with $\bm{K}_{ij}=\mathcal{K}({\bf x}_{i},{\bf y}_{j})$. Such type of computation arises in many science and engineering fields, such as kernel methods in statistical learning and machine learning [1, 2], data assimilation methods in geosciences [3, 4], particle simulations and boundary integral/element methods in computational physics [5, 6, 7], dislocation dynamics simulations in material science [8, 9], etc. To compute $\bm{\phi}$ in Eq. 2, a naive direct evaluation requires $\mathcal{O}(N^{2})$ operations, which is prohibitive when $N$ is large.

1.1 Related work

One special but important instance of Eq. 1 is when the data points $\{{\bf x}_{i}\}_{i=1}^{N}$ and $\{{\bf y}_{i}\}_{i=1}^{N}$ lie on a regular grid. In this case, Eq. 1 can be evaluated exactly (up to round-off errors) using the fast Fourier transform (FFT), which requires $\mathcal{O}(N\log(N))$ operations. Although the FFT can be extended to handle non-uniform data distributions [10, 11], but its efficiency decreases for highly irregular distributions in three dimensions (3D).

The fast multipole method (FMM) is a general framework that has been successfully applied to non-uniform data distributions. In the FMM, we partition the problem domain hierarchically into subdomains of different scales and exploit the multi-scale decomposition in the evaluation of Eq. 1. In particular, we evaluate exactly the calculation associated with adjacent subdomains at the finest scale; and we evaluate approximately the calculation associated with non-adjacent subdomains at every scale. Overall, the procedure requires $\mathcal{O}(N)$ operations with a provable accuracy. While the FMM has been derived for specific kernels that appear frequently in computational physics [5, 12, 13, 14, 15], the derivation may be tedious, difficult or even impossible for an arbitrary NOTIS kernel function.

Therefore, black-box algorithms have been developed, which require only the evaluation of a given kernel function. Such methods can be classified into two groups. The first consists of methods that approximate the kernel function (away from the origin) with polynomials, such as Legendre polynomials or Chebyshev polynomials [16, 17, 18, 9]. The other group consists of methods that compute the so-called equivalent densities or the so-called skeletons for every subdomain to efficiently represent the contained source points and their weights [19, 20, 21, 22]. Theoretically, this approach is justified by the potential theory for kernel functions that are fundamental solutions of non-oscillatory elliptic partial differential equations.

1.2 Contributions

In this paper, we present a parallel implementation of the black-box method in [18, 9] for evaluating Eq. 1 using $\mathcal{O}(N)$ memory and computation. The key idea is using Lagrange interpolation to construct approximations of the kernel function $\mathcal{K}({\bf x},{\bf y})$ when ${\bf x}$ and ${\bf y}$ are distant. Unlike fast algorithms that have been developed for specific kernel functions, our method applies to a wide range of functions. Successful stories include applications of our method in dislocation dynamics simulations [9] and aquifer characterization [23]. In those two applications, no fast algorithm exists for the two kernel functions—the Green’s function in anisotropic elasticity and the isotropic exponential function. In particular, our method requires only a black-box routine to evaluate the kernel function, and thus can be integrated easily with other codes. For example, our package was used in an iso-geometric boundary element method to obtain a solver of $\mathcal{O}(N)$ complexity [7]. Other examples of using (an earlier version of) our code are in the elastic formulation of the displacement discontinuity method for the simulation of micro-seismicity [24, 25, 26]. Extension of our algorithm for solving/factorizing kernel matrices has been explored in [27, 28].

To evaluate Eq. 1, our method follows the general FMM machinery as follows. First, the problem domain containing $\{{\bf x}_{i}\}_{i=1}^{N}$ and $\{{\bf y}_{i}\}_{i=1}^{N}$ is partitioned in a hierarchical fashion (see Fig. 1(a)), and the partitioning is associated with a tree data structure where every tree node represents a subdomain in the hierarchy. Second, a post-order traversal of the tree is performed, and the “multipole coefficients” associated with every tree node is computed as a compact representation of the source points and corresponding weights that the subdomain contains. This step is often called an upward pass. Third, at every level of the hierarchy, the “local coefficients” of every tree node is computed using the “multipole coefficients” of its interaction list (two nodes are in each other’s interaction list if their parents are adjacent but they are not; see Fig. 1(b)). Fourth, a pre-order traversal of the tree is performed to “accumulate” the local coefficients of all nodes to those at the leaf level. This step is often called a downward pass. Finally, the contribution from adjacent subdomains are evaluated exactly for all leaf nodes. Technically speaking, multipole and local coefficients are terminologies used in the original FMM [5]. Here we adhere to the same terms for their counterparts in our algorithm, which we will rigorously define in Section 2.3.

Our parallel algorithm is based on three observations of the above procedure: (1) the computation of multipole coefficients and local coefficients during the upward pass and the downward pass is embarrassingly parallel for tree nodes at the same level; (2) the computation of local coefficients based on multipole coefficients is embarrassingly parallel for all nodes; and (3) the exact evaluation of the contribution from adjacent subdomains is embarrassingly parallel for all leaf nodes. While our parallel algorithm does not exploit task-level parallelism [29], which usually requires a well-designed runtime system to implement, our results show satisfactory parallel speedups.

To summarize, we present parallel black-box FMM in 3D (PBBFMM3D) to compute the kernel matrix-vector multiplication in Eq. 1 for NOTIS kernels. The algorithm provably requires $\mathcal{O}(N)$ memory and operations. Our parallel implementation is based on OpenMP111https://www.openmp.org/ and targets multi-core shared-memory machines. Nowadays, on-node parallelism is becoming increasingly important because of the wide adoption of multi-core architectures. Our method also serves as the basis for a multi-node distributed-memory solution, which can be implemented on top of our code using the so-called local essential tree and space-filling curve [30, 31]. Extensions to many-core architectures (e.g., GPUs) was considered in [32, 33], which used blocking schemes to reduce memory movement and to improve the arithmetic intensity (flop-to-word ratio). Our code is publicly available at

https://github.com/ruoxi-wang/PBBFMM3D

1.3 Outline

The rest of this paper is organized as follows. Section 2 describes the algorithm focusing on the data dependency and the parallelism. Section 3 describes the software architecture and the user interface. Section 4 presents the accuracy and the running time of PBBFMM3D. Section 5 draws the conclusion.

2 Parallel black-box algorithm

In this section, we present the black-box algorithm for evaluating Eq. 1. In particular, we focus on the data dependency and illustrate our parallel strategy. Although our code deals with 3D problems, we use 1D and 2D examples here for pictorial illustration.

2.1 Hierarchical domain decomposition

As with other multilevel methods, our approach is based on a hierarchical decomposition of the problem domain to achieve linear complexity. Specifically, a cubical domain is divided into eight subdomains through binary partitioning along every coordinate. Then, every subdomain is divided recursively until the number of data points in every subdomain is less than a prescribed constant. This hierarchical domain decomposition is naturally associated with an octree tree data structure, where the root stands for the entire domain, and the other nodes stand for subdomains at different levels in the hierarchy; see Fig. 1(a). For the rest of the paper, we use the terms subdomain and tree node interchangeably.

Note the hierarchical decomposition is generally non-uniform ($\mathcal{T}$ is adaptive). But we choose to use a uniform-tree data structure in PBBFMM3D for the ease of parallel implementation, which is reasonably efficient as long as the point distribution is not extremely irregular. With a uniform tree, some leaf nodes may end up having few data points. In PBBFMM3D, empty nodes are skipped in the algorithm, and the overhead of using a uniform tree compared with an adaptive one is from processing nodes that have only a few points.

Given the tree $\mathcal{T}$, the parent $\mathcal{P}(\alpha)$ and the children $\mathcal{C}(\alpha)$ of a node $\alpha$ are naturally defined. We also define the neighbors $\mathcal{N}(\alpha)$ and the interaction list $\mathcal{I}(\alpha)$ of a node $\alpha$ as below; see Fig. 1(b) for an pictorial illustration.

Definition 2.1** (Neighbors and interaction list).**

Given a hierarchical tree structure,

• •

$\mathcal{N}(\alpha)$: the adjacent nodes (subdomains) of $\alpha$ at the same level in the hierarchy including $\alpha$ itself. The number of neighbors is generally $3^{d}$, where $d$ is the spatial dimension.

• •

$\mathcal{I}(\alpha)=\mathcal{C}(\mathcal{N}(\mathcal{P}(\alpha)))/\mathcal{N}(\alpha)$: the non-adjacent nodes (subdomains) at the same level in the hierarchy whose parents are neighbors of $\mathcal{P}(\alpha)$. The number of nodes in the interaction list is generally $6^{d}-3^{d}$, where $d$ is the spatial dimension.

2.2 Separation of variables (low-rank approximation)

The key idea of our algorithm is to approximate the kernel function $\mathcal{K}({\bf x},{\bf y})$ through polynomial interpolation when the target point ${\bf x}$ and the source point ${\bf y}$ are distant. In this section, we focus on the situation where ${\bf x}$ and ${\bf y}$ are inside two non-adjacent subdomains, respectively. Technically speaking, ${\bf x}$ and ${\bf y}$ are well-separated. For the following discussion, we need a set of $p$ interpolation nodes $\mathcal{S}=\{x^{*}_{1},x^{*}_{2},\ldots,x^{*}_{p}\}$ on the real line. Then we can form interpolation nodes in 3D with tensor products:

[TABLE]

where $i$ and $(i_{1},i_{2},i_{3})$ are 1D-index and 3D-index of the $p\times p\times p$ grid, respectively, e.g., $i_{1}=i\,\%\,p$, $i_{2}=i\mod p$, and $i_{3}=i\mod p^{2}$.

Definition 2.2** (Lagrange basis polynomials).**

The $p$-th order Lagrange basis polynomials in 1D are

[TABLE]

The $p$-th order Lagrange basis polynomials in 3D are tensor products of the $p$-th order Lagrange basis polynomials in 1D:

[TABLE]

where ${\bf x}=(x_{1},x_{2},x_{3})\in\mathbb{R}^{3}$, and ${\bf x}^{*}_{i}=(x^{*}_{i_{1}},x^{*}_{i_{2}},x^{*}_{i_{3}})\in\mathcal{S}\otimes\mathcal{S}\otimes\mathcal{S}$.

Definition 2.3** ($p$-th order polynomial interpolant of $\mathcal{K}({\bf x},{\bf y})$).**

Denote $\hat{\mathcal{K}}({\bf x},{\bf y})$ as the $p$-th order polynomial interpolant of $\mathcal{K}({\bf x},{\bf y})$, which is obtained through Lagrange interpolation on ${\bf x}$ and ${\bf y}$, respectively (${\bf x}$ and ${\bf y}$ are well-separated):

[TABLE]

where ${\bf x}^{*}_{i},{\bf y}^{*}_{j}\in\mathcal{S}\otimes\mathcal{S}\otimes\mathcal{S}$.

Observe that constructing the approximation $\hat{\mathcal{K}}({\bf x},{\bf y})$ requires only evaluations of the kernel function $\mathcal{K}({\bf x}^{*}_{i},{\bf y}^{*}_{j})$.

The PBBFMM3D code offers two options for interpolation nodes, namely, Chebyshev nodes

[TABLE]

and equally spaced nodes

[TABLE]

With the Chebyshev nodes, the approximation in Definition 2.3 is nearly optimal among polynomials of the same order. More importantly, the error decays as $\mathcal{O}(\rho^{-p})$ if the kernel function $K({\bf x},{\bf y})$ is analytic and bounded in the “Bernstein ellipse” of foci 1 and -1 with semimajor and semiminor axis lengths summing to $\rho$ [34]. With equally spaced nodes, the matrix $\mathcal{K}({\bf x}^{*}_{i},{\bf y}^{*}_{j})$ in Definition 2.3 is a block-Toeplitz-Toeplitz-block matrix, which has a reduced memory footprint and can be applied in $\mathcal{O}(p^{3}\log(p))$ time using the FFT [9]. Note the Lagrange interpolation of high degree over equally spaced nodes does not always converge, even for smooth functions, which is known as Runge’s phenomenon. In practice, we find low-order approximations sufficiently accurate in many applications. The scheme of Lagrange interpolation on uniform nodes can be stabilized by fitting a polynomial of degree $d<2\sqrt{p}$ using least-squares.

2.3 Black-box FMM algorithm

To evaluate Eq. 1, our algorithm has the following four stages.

Upward pass. A post-order traversal of $\mathcal{T}$ is performed to compute the “multipole coefficients” of every subdomain, which encodes information of the source points and their weights contained in the subdomain. For every leaf node $\alpha$ in $\mathcal{T}$, the particle-to-moment (P2M) translation is executed:

[TABLE]

where ${\bf y}_{j}$ and ${\bf y}^{*}_{i}$ are the source points and interpolation nodes in $\alpha$, respectively. For every non-leaf node $\alpha$ in $\mathcal{T}$, the moment-to-moment (M2M) translation is executed:

[TABLE]

where ${\bf y}^{*}_{i}$ and ${\bf y}^{*}_{j}$ are the interpolation nodes in $\alpha$ and $\beta$, respectively. 2. 2.

Far-field interaction. The “local coefficients” of every node is computed using the “multipole coefficients” of its interaction list. For every node $\alpha$ in $\mathcal{T}$, the moment-to-local (M2L) translation is executed:

[TABLE]

where ${\bf x}^{*}_{i}$ and ${\bf y}^{*}_{j}$ are the interpolation nodes in $\alpha$ and $\beta$, respectively. 3. 3.

Downward pass. A pre-order traversal of $\mathcal{T}$ is performed to “accumulate” the “local coefficients” at leaf nodes. This is effectively a “transpose” of the upward pass. For every node $\alpha$ in $\mathcal{T}$, the local-to-local (L2L) translation is executed:

[TABLE]

where $\beta=\mathcal{P}(\alpha)$, ${\bf x}^{*}_{i}$ and ${\bf x}^{*}_{j}$ are the interpolation nodes in $\alpha$ and $\beta$, respectively. For every leaf node $\alpha$ in $\mathcal{T}$, the local-to-particle (L2P) translation is also executed:

[TABLE]

where ${\bf x}_{i}$ and ${\bf x}^{*}_{j}$ are the target points and interpolation nodes in $\alpha$, respectively. 4. 4.

Near-field interaction. The contribution from neighbors is evaluated exactly. For every leaf node $\alpha$, the particle-to-particle (P2P) translation is executed:

[TABLE]

where ${\bf x}_{i}$ and ${\bf y}_{j}$ are the target and the source points in $\alpha$ and $\beta$, respectively.

To summarize, the contribution from neighbors is calculated exactly; the contribution from the interaction list is approximated using the polynomial interpolant in Definition 2.3; and the contribution from the remaining leaf nodes are approximated through coarser levels in the tree. A diagram illustrating the above process is shown in Fig. 2. Although our algorithm is presented using scalar operations for ease of illustration, it is implemented using BLAS2 and BLAS3 subroutines.

Below, we present a theorem for the computational cost and memory footprint of the above algorithm. Our primary focus is showing the complexity in terms of the number of points $N$, and thus we present the computational cost of one M2L translation as $\mathcal{O}(p^{6})$, the cost of a naive matrix-vector product. For the same reason, we use an $\mathcal{O}(p^{6})$ estimate for the storage of every M2L translation operator. Given the M2L translation is a bottleneck in the FMM, we implement acceleration techniques in PBBFMM3D, which is discussed in Section 2.5.

Theorem 2.1** (Computational cost and memory footprint).**

The computational cost and memory footprint of the black-box FMM algorithm are both $\mathcal{O}(N)$.

Suppose every leaf node has at most $n_{0}$ points (typically $64\sim 128$). The number of leaf nodes and the number of tree nodes are both $\mathcal{O}(N/n_{0})$. In the above algorithm, every node requires constant amount of work (independent of $N$) at every stage:

Upward pass. P2M: $\mathcal{O}(p^{3}\,n_{0})$ work for every leaf node. M2M: $\mathcal{O}(8\,p^{4})$ work for every non-leaf node, where every non-leaf node has at most 8 children in 3D. Note the M2M is a 3D tensor-vector multiplication. 2. 2.

Far-field interaction. M2L: $\mathcal{O}(189\,p^{6})$ work for every node, where the interaction list has at most $6^{3}-3^{3}=189$ nodes in 3D. 3. 3.

Downward pass. L2L: $\mathcal{O}(p^{4})$ work for every node. L2P: $\mathcal{O}(p^{3}\,n_{0})$ work for every leaf node. Note the L2L is a 3D tensor(transpose of the tensor in M2M)-vector multiplication. 4. 4.

Near-field interaction. P2P: $\mathcal{O}(27\,n_{0}^{2})$ work for every leaf node, where every node generally has at most $3^{3}=27$ neighbors in 3D.

Therefore, the computational complexity of the entire algorithm is $\mathcal{O}(N)$.

Regarding the memory footprint, every tree node stores $p^{3}$ “multipole coefficients” and $p^{3}$ “local coefficients”, which sums up to $\mathcal{O}(Np^{3}/n_{0})$. Since the kernel function is translational invariant, we need to precompute only $7^{3}-3^{3}=316$ M2L translation operators at every level for the $\log(N)$-level tree structure, which is $\mathcal{O}(316p^{6}\log(N))$ memory in total. Therefore, the memory footprint is $\mathcal{O}(2p^{3}N+316p^{6}\log(N))$.

2.4 Parallel algorithm

In this section, we analyze the parallelism in each of the four stages in the FMM algorithm, and we have implemented a parallel algorithm using the OpenMP API for shared-memory machines.

Upward pass. As stated earlier, the upward pass is a post-order traversal of $\mathcal{T}$, so a parallel post-order tree traversal using OpenMP tasks is implemented in PBBFMM3D. 2. 2.

Far-field interaction. The M2L translations are independent for all nodes. In addition, the translation typically requires the same amount of work for every node. Therefore, the OpenMP “parallel for” directive is employed on the loop over all nodes for M2L translations in PBBFMM3D. 3. 3.

Downward pass. As stated earlier, the downward pass is a pre-order traversal of $\mathcal{T}$, so a parallel pre-order tree traversal using OpenMP tasks is implemented in PBBFMM3D. 4. 4.

Near-field interaction. The P2P translations are independent for all leaf nodes. However, the translations between pairs of neighbors are work-heterogeneous due to the non-uniform distribution of the target and the source points. In PBBFMM3D, the OpenMP “parallel for” directive is employed on the loop over all leaf nodes for P2P translations, and different scheduling policies can used based on any prior knowledge of the point distribution.

Algorithm 1 shows the pseudocode of the above parallel algorithm. A common alternative to the parallel tree traversal using OpenMP task is a level-by-level traversal with the OpenMP “parallel for” directive on the loop over all nodes at the same level. Although the P2M and the L2P translations are work-heterogeneous, the efficiency of the alternative approach may still be reasonable since the upward pass and the downward pass are usually not the bottleneck.

Here our parallel algorithm focuses on parallelizing each stage of the FMM algorithm, and haven’t exploit the concurrency across different stages. In principle, the near-field interaction stage does not depend on the others except for updating the results. So the work-heterogeneous P2P translations can be prioritized. Furthermore, notice the M2L translation of a node can happen as soon as its “multipole coefficients” have been computed. However, implementing these ideas efficiently typically requires a task-based runtime system, for which we refer interested readers to [29, 33].

2.5 Acceleration techniques

Since the far-field interaction and the near-field interaction usually dominate the entire computation, we introduce the acceleration techniques used in PBBFMM3D.

Homogeneous kernel

A kernel function is homogeneous if

[TABLE]

for $\alpha\not=0$ and $m$ is typically an integer. For example, $\mathcal{K}({\bf x},{\bf y})=1/\|{\bf x}-{\bf y}\|$ is a homogeneous kernel function of degree -1. Since the interpolation grids are fixed relative to the problem domain, the M2L translation operators of different levels differ only by a scaling constant. Hence, we store these operators for only the leaf level.

Symmetry and skew-symmetry

A kernel function is symmetric if

[TABLE]

and is skew-symmetric if

[TABLE]

In many applications, the source $\{{\bf x}_{i}\}_{i=1}^{N}$ and the target $\{{\bf y}_{i}\}_{i=1}^{N}$ are the same set of points. Therefore, the kernel matrix $\bm{K}$ becomes symmetric or skew-symmetric if the kernel function is such. This implies only one of the two P2P/M2L translation operators needs be stored between a pair of tree nodes that are either neighbors or in each other’s interaction list.

Fast M2L translation

Recall the definition of an M2L translation operator $\mathcal{K}({\bf x}^{*},{\bf y}^{*})$ in Definition 2.3, where ${\bf x}^{*}$ and ${\bf y}^{*}$ are interpolation grids in a pair of tree nodes that are in each other’s interaction list. In PBBFMM3D, there are two options for the interpolation grids: Chebyshev nodes or equally spaced nodes. With Chebyshev nodes, an SVD-based compression of the M2L operator is employed following the approach in [18] since the translation operator is observed to be numerically low rank. With equally spaced nodes, the M2L operator is a block-Toeplitz-Toeplitz-block matrix [9], where the $p^{3}$-by-$p^{3}$ matrix has only $(2p-1)^{3}$ unique entries and can be applied to a vector in $\mathcal{O}(p^{3}\log(p))$ time using the FFT.

Multiple right-hand-sides

In some applications, Eq. 1 needs to be evaluated with multiple weight vectors associated with the same set of source points. In PBBFMM3D, all weight vectors are grouped into a matrix as the input of the FMM algorithm, which allows using cache-friendly BLAS3 operations in our algorithm.

3 Software description.

In this section, we briefly discuss the software architecture of PBBFMM3D focusing on the black-box feature of the algorithm and the C++ and Python interfaces. More details can be found in the documentation at https://github.com/ruoxi-wang/PBBFMM3D. The code is written in C++ with the OpenMP API, and requires some standard linear algebra libraries including the BLAS222http://www.netlib.org/blas/, the LAPACK333http://www.netlib.org/lapack/ and the FFTW3 library444http://www.fftw.org/. The Boost Python Libraries555https://www.boost.org/doc/libs/1_70_0/libs/python/doc/html/index.html is required for using the Python interface.

The PBBFMM3D has the following three main classes as shown in Fig. 3.

• •

Class H2_3D_Tree sets parameters and creates the hierarchical partitioning of the problem domain. The parameters include (1) “Domain size”: side length of a cubical problem domain, (2) “Tree level”: the number of levels in the hierarchical partitioning, (3) “Interpolation type”: Chebyshev nodes or equally spaced nodes, (4) “Interpolation order”: the number of interpolation nodes used in Definition 2.3, (5) “SVD truncation error”: error from the compression of the M2L translation operators, which is by default the prescribed accuracy of the entire computation. The two key member functions are the following. Function PrecomputeM2L() precomputes the M2L operators, and BuildFMMHierarchy() creates the hierarchical partitioning and builds the corresponding data structure.

• •

Class H2_3D_Compute stores the information regarding the source and the target points including (1) “Target”: position of target points $\{{\bf x}_{i}\}_{i=1}^{N}$, (2) “Source”: position of source points $\{{\bf y}_{i}\}_{i=1}^{N}$, (3) “Weight”: weights $\{\sigma_{i}\}_{i=1}^{N}$ associated with the source points, and (4) “Number of weights”: number of weights associated with every source point. The five key member functions include FMMDistribute(), which assigns the source and the target points to leaf cells in the tree, and four functions correspond to the four translation stages described in Section 2.

• •

Class kernel_NAME defines the kernel function. The member function SetKernelProperty() sets the homogeneous and the symmetric properties of the kernel function. The member function EvaluateKernel() takes two data points ${\bf x}$ and ${\bf y}$ and returns the value of $\mathcal{K}({\bf x},{\bf y})$.

C++ & Python interfaces

Listing LABEL:lst:cpp and LABEL:lst:python show the core lines of a basic example using the C++ interface and the Python interface, respectively. The example evaluates Eq. 1 with the standard Gaussian kernel. The first line creates an object of Class kernel_Gaussian, which implements the standard Gaussian function. The class inherits from class H2_3D_Tree and takes input parameters. The second line creates the hierarchical partitioning of the problem domain. The last line does the computation and stores results in the output variable.

Customized kernel

While some commonly used kernel functions are already implemented in PBBFMM3D, defining a new kernel function is straightforward as shown in Listing LABEL:lst:myKernel. It requires implementing only two methods as follows. The EvaluateKernel() method takes a pair of source and target points and returns the function value, and the SetKernelProperty() method tells the homogeneous degree of the kernel and whether it is symmetric (see Section 2.5).

4 Numerical Results

In this section, we present numerical experiments with PBBFMM3D to show the accuracy, the sequential running time and the parallel scalability. A real-world application in geostatistics is also presented, where PBBFMM3D was used to speed up the calculations with a covariance matrix. We focus on two particular kernel functions here: ${1}/{\|{\bf x}-{\bf y}\|}$ and $\text{exp}(-\|{\bf x}-{\bf y}\|)$, where ${1}/{\|{\bf x}-{\bf y}\|}$, the Green’s function for the Laplace equation in 3D, is frequently used in computational physics; and $\text{exp}(-\|{\bf x}-{\bf y}\|)$ is a popular choice as the covariance function for a Gaussian process.

All experiments were performed on a Linux server with 192 GB of RAM and the Intel Xeon Platinum 8280 (“Cascade Lake”) with 56 cores on two sockets (28 cores/socket). The code was compiled with GCC 8.3.0, which implements version 4.5 of the OpenMP standard, and the code was linked with the Intel MKL library666https://software.intel.com/en-us/mkl version 2020.1.217 and the FFTW3 library version 3.3.8.

Parameters and notations

• •

$N$: the number of source/target points, which are randomly generated using the uniform distribution in the unit cube.

• •

$p$: the order of interpolation, which determines the accuracy of our algorithm. Note the number of interpolation nodes is $p^{3}$ in Definition 2.3 as a tensor product of $p$ nodes in each dimension.

• •

unif: using equally spaced nodes in Definition 2.3.

• •

cheb: using Chebyshev nodes in Definition 2.3.

• •

tree_level: the number of levels of the hierarchical partitioning of the problem domain.

• •

SVD accuracy: the compression accuracy of the M2L translation operator, which is chosen to be the same as the prescribed accuracy.

• •

$r$: notation of the Euclidean distance between a pair of source point ${\bf x}$ and a target point ${\bf y}$, i.e., $r=\|{\bf x}-{\bf y}\|_{2}$.

4.1 Accuracy & sequential running time

In this section, we focus on two parameters $N$ and $p$. First, we fix $N=10^{4}$ and increase $p$ to show the accuracy of PBBFMM3D. (We refer to [9] for the precomputation time and the memory footprint.) Then, we fix $p=4$ and increase $N$ to show the sequential running time.

Fig. 5 shows the (relative) error of evaluating Eq. 1 as a function of $p$, where the kernel functions are $1/r$ and $\text{exp}({-r})$, respectively. The error is defined as

[TABLE]

where $\hat{K}$ stands for the approximation constructed implicitly in PBBFMM3D. Since both kernel functions are analytic (away from the origin), the error of our interpolation using either Chebyshev nodes or uniform nodes decays exponentially.

Fig. 5 shows the scaling of the sequential running time with respect to the number of points $N$. We see that the time increases linearly as $N$ increases, as opposed to the quadratic increase of evaluating Eq. 1 naively.

4.2 Parallel scalability

In this section, we present the parallel running time of PBBFMM3D on up to 32 cores. We focus on the kernel function $1/r$ and use Chebyshev nodes in Definition 2.3 with $p=4$. Table 1 reports the parallel running time. To show the parallel scalability, we chose the number of cores to be a power of 2. But that is not required in PBBFMM3D. The baseline (1 core) is the serial version used in [18, 9]. As the table shows, we obtained approximately $19\times$ speedup on 32 cores for $N=8^{4}\times 10^{4}$.

Fig. 6(a) shows the strong scalability of PBBFMM3D, i.e., the running time using a sequence of increasing number of cores for a fixed problem size. We provided a breakdown of the running time into the four stages: upward pass, far-field interaction, downward pass, and near-field interaction. We see that the running time for all stages nearly halved as the number of cores doubled.

Fig. 6(b) shows the weak scalability of PBBFMM3D, i.e., the running time for a fixed problem size per core. So we increased the number of particles proportionally to the number of cores. As the figure shows, the time spent on each stage in the FMM only increased by a small amount when the problem size increased by $8\times$ (the number of cores also increased by $8\times$). The ideal runtime would stay unchanged due to the linear complexity of PBBFMM3D. Our implementation achieved an efficiency of 73% to 86% (serial time/parallel time).

4.3 Application in Gaussian Processes

Gaussian random field (GRF) theory [35, 36] has been widely used in interpolation and estimation of spatially correlated unknowns. For example, GRF methods can be used for estimating the permeability of the underground soil and rock, which is of critical interest to hydrogeologists and petroleum engineers [37, 38]. In GRF methods, a covariance matrix is required as the prior information of the underlying unknown field. However, practical applications typically require large numbers of unknowns, so dimension reduction techniques such as the principle component analysis are required.

To obtain the top-$k$ principle components, we employ PBBFMM3D to compute the truncated eigenvalue decomposition of the covariance matrix. In particular, we use a randomized method (Algorithm 5.3 in [39]) to calculate the top-$k$ eigenvalues and their associated eigenvectors, and the randomized algorithm requires evaluating Eq. 1 $\mathcal{O}(k)$ times (same target and source points but with $\mathcal{O}(k)$ different weights). While the original method requires $\mathcal{O}(kN^{2})$ operations, we accelerate the method with PBBFMM3D and arrive at the optimal $\mathcal{O}(kN)$ complexity, where $N$ is the number of data points.

In the next experiment, we randomly generated data points in the unit cube and employed the kernel function $e^{-r}$, one-dimensional Matérn kernel with smoothness 1/2. The goal was to compute the top-100 eigen-pairs of the covariance matrix, where 120 matrix-vector products (packed into one matrix-matrix product) were evaluated twice in the randomized algorithm. In the original randomized method, matrix-matrix products are computed by calling the dgemm subroutine in the Intel MKL library (with multi-threading on 32 cores). In the PBBFMM3D-accelerated method, the number of Chebyshev nodes was $p=4$, and the number of tree levels were 3, 4 and 5 for the three increasing problem sizes. We measured the error $\|\Lambda_{FMM}-\Lambda_{ext}\|_{2}/\|\Lambda_{ext}\|_{2}$, where $\Lambda_{ext}$ and $\Lambda_{FMM}$ are the eigenvalues computed using the original method and the PBBFMM3D-accelerated method, respectively. The errors for $N=10^{4}$ and $N=8\times 10^{4}$ are $1.6\times 10^{-4}$ and $2.8\times 10^{-5}$, respectively.

Fig. 7 shows the running time of the original method and our accelerated method. In Geostatistics applications, the original randomized method with exact matrix-matrix products is the common practice, and therefore we set that as our baseline. This allows us to show the best speedup and scaling one can achieve from leveraging PBBFMM3D. In the original method, we need to evaluate all entries of the covariance matrix and compute matrix-matrix products, both of which scale as $\mathcal{O}(N^{2})$. It is worth noting that when $N=64\times 10^{4}$, forming the entire covariance matrix requires 3.2 TB of memory.

In the PBBFMM3D-accelerated method, we need to evaluate only $\mathcal{O}(N)$ entries of the covariance matrix. As Fig. 7 shows, our method scales linearly with respect to the problem size $N$, and thus the speedup over the original method becomes more pronounced when $N$ increases.

In the last experiment, we apply PBBFMM3D to data points lying on an unstructured grid, as shown in Figure 8. Same as above, we compute the top-50 eigenvectors with the randomized algorithm accelerated by PBBFMM3D. For this experiment, we employ the kernel function $e^{-r}$, and we used $p=5$ Chebyshev nodes in PBBFMM3D.

5 Conclusions

We have introduced PBBFMM3D, a black-box algorithm/software for evaluating kernel matrix-vector multiplication on shared-memory machines. The target kernels are non-oscillatory translation invariant functions that are sufficiently smooth away from the origin. The user needs only provide a function routine that returns the kernel value given a source point and a target point, if the kernel function is not already implemented. (The user can also specify the degree of homogeneity and whether the kernel is symmetric or skew-symmetric to active specific optimizations.) Our algorithm requires $\mathcal{O}(N)$ memory and work, where $N$ is the number of data points.

A parallel algorithm is presented in this paper and implemented using OpenMP for shared-memory machines. We have presented parallel scalability results on up to 32 cores and achieved at most $19\times$ speedup. We have also presented an application in geostatistics, where we accelerated the computation of the truncated eigen-decomposition of covariance matrices.