# Data-Driven Construction of Hierarchical Matrices with Nested Bases

Difeng Cai Department of Mathematics, Emory University, Atlanta, GA 30322
<EMAIL_ADDRESS>yxi26@emory.edu). The research of Difeng Cai and Yuanzhe Xi
is supported by NSF award OAC 2003720 and RTG Grant DMS-2038118. Hua Huang
School of Computational Science and Engineering, Georgia Institute of
Technology, Atlanta, GA 30332<EMAIL_ADDRESS>echow@cc.gatech.edu). The
research of Hua Huang and Edmond Chow is supported by NSF award OAC 2003683.
Edmond Chow† Yuanzhe Xi∗

###### Abstract

Hierarchical matrices provide a powerful representation for significantly
reducing the computational complexity associated with dense kernel matrices.
For general kernel functions, interpolation-based methods are widely used for
the efficient construction of hierarchical matrices. In this paper, we present
a fast hierarchical data reduction (HiDR) procedure with $O(n)$ complexity for
the memory-efficient construction of hierarchical matrices with nested bases
where $n$ is the number of data points. HiDR aims to reduce the given data in
a hierarchical way so as to obtain $O(1)$ representations for all nearfield
and farfield interactions. Based on HiDR, a linear complexity
$\mathcal{H}^{2}$ matrix construction algorithm is proposed. The use of data-
driven methods enables better efficiency than other general-purpose methods
and flexible computation without accessing the kernel function. Experiments
demonstrate significantly improved memory efficiency of the proposed data-
driven method compared to interpolation-based methods over a wide range of
kernels. Though the method is not optimized for any special kernel, benchmark
experiments for the Coulomb kernel show that the proposed general-purpose
algorithm offers competitive performance for hierarchical matrix construction
compared to several state-of-the-art algorithms for the Coulomb kernel.

##  1 Introduction

In various applications, the pairwise interaction between two objects is
characterized by a non-local kernel function. A system of $n$ objects then
gives rise to a $n$-by-$n$ dense kernel matrix. Such matrices arise frequently
in integral equations [22, 40, 41], astrophysics [4], statistics [26, 18],
machine learning [7, 23], etc. A computational challenge is that the naive
computational or storage cost associated with the dense kernel matrix is at
least $O(n^{2})$. For the Coulomb kernel and its variants, pioneering work
such as the Fast Multipole Method (FMM) [40, 41, 28] and the Barnes-Hut
algorithm [4] use multilevel approximation to successfully reduce the cost to
linear or quasilinear complexity. These techniques were later generalized into
the powerful algebraic framework of hierarchically low-rank matrices [33, 32,
9, 30] for efficiently approximating general dense kernel matrices with nearly
optimal computational cost. Two widely used classes of hierarchical matrices
are $\mathcal{H}$ matrices and $\mathcal{H}^{2}$ matrices. The $\mathcal{H}$
matrix yields an $O(n\log n)$ representation and the $\mathcal{H}^{2}$ matrix
yields an optimal $O(n)$ representation for approximating the kernel matrix
and computing the matrix-vector multiplication. For dense kernel matrices, the
efficient construction of these hierarchical representations associated with
general kernels remains a challenging problem. In this paper, we address this
issue for the $\mathcal{H}^{2}$ matrix representation. For general kernels,
interpolation-based methods are often adopted as a black-box tool used in
hierarchical matrix construction. However, as will be demonstrated in later
sections, the use of interpolation nodes (or points outside the given dataset)
may lead to loss of accuracy. We present a general-purpose data-driven
framework for hierarchical matrix construction, that resolves these issues and
offers improved efficiency.

Given a set of points and a kernel function, hierarchical matrix construction
starts with an adaptive partitioning of the data. The partitioning can be
encoded by a tree structure in which each node represents a subset generated
in the adaptive partitioning procedure. For example, the root node corresponds
to the entire set of $n$ points and its children nodes correspond to the
subsets generated from the first partitioning of the dataset. For hierarchical
matrix construction, directly compressing the submatrix corresponding to a
pair of nodes/subsets will be inefficient, since one of the subsets may
contain $O(n)$ points (see Figure 1(a)) and the total cost for all such pairs
will be at least $O(n^{2})$. We propose to first process the tree-structured
data to obtain a reduced representation in which each node only corresponds to
a small subset with $O(1)$ points while the farfield corresponds to $O(1)$
points as well. These $O(1)$ subsets are called _representor sets_ in [6]. See
Figure 1(b). Using representor sets, the multilevel compression of the kernel
matrix can be rapidly computed. To achieve optimal efficiency, we design a
hierarchical data reduction procedure with computational cost of $O(n)$.
Different from existing approaches, the procedure operates entirely on the
given data without accessing the kernel function and no algebraic compression
is performed. Hence the approach is termed _data-driven_. We show how to
incorporate hierarchical data reduction in hierarchical matrix construction to
obtain a fast algorithm for general kernels. Numerical experiments demonstrate
the competitive performance of the new method in terms of generality, memory
use, speed, and accuracy.

(a) Left: the farfield (shown in blue) of the orange box contains $O(n)$
points. Right: the $O(n)$ points are reduced to a subset of $O(1)$ points.
(For another orange box, the farfield is different and the reduced subset is
generally a different set of points.)

(b) The kernel matrix is shown, highlighting the interaction between points in
the orange region and points in its farfield. The arrows point to columns
corresponding to the points selected in the $O(1)$ subset. These columns form
an approximate basis for the farfield.

Figure 1: Farfield data reduction as a way to construct an approximate column
basis for the farfield block row in the kernel matrix

The rest of the manuscript is organized as follows. Section 2 reviews
hierarchical matrix structures and existing methods for constructing
hierarchical matrices. Section 3 introduces the hierarchical data reduction
algorithm. Based on this hierarchical data reduction, the complete
hierarchical matrix construction is presented in Section 4. Section 5 presents
numerical experiments to investigate different data reduction techniques.
Numerical results for the proposed data-driven hierarchical matrix
construction are given in Section 6. Concluding remarks are drawn in Section
7.

The notation used in this paper is listed below.

  * •

$|x-y|$ denotes the Euclidean distance between $x$ and $y$ in
$\mathbb{R}^{d}$;

  * •

diam($X$) denotes the diameter of set $X$, i.e., $\max\limits_{x,y\in
X}|x-y|$;

  * •

dist($X_{1},X_{2}$) denotes the distance between sets $X_{1}$ and $X_{2}$,
i.e. $\min\limits_{x\in X_{1},y\in X_{2}}|x-y|$;

  * •

card($X$) denotes the cardinality of set $X$.

  * •

$X^{*}$ denotes a representor set of $X$.

##  2 Review of Hierarchical Matrix Representations

Given a kernel function $\kappa(x,y)$ and a dataset
$X=\\{x_{1},\dots,x_{n}\\}$, the associated kernel matrix is defined by

$K=[\kappa(x_{i},x_{j})]_{i,j=1}^{n}.$

Dense kernel matrices are ubiquitous and arise in various applications, where
the kernel function measures the interaction between objects. In Coulombic
N-body simulations, $\kappa(x,y)=\frac{1}{|x-y|}$. In boundary integral
equations, $\kappa(x,y)$ can take very different forms, including $\Phi(x,y)$,
$\nabla_{v_{y}}\Phi(x,y)$, where $\Phi(x,y)$ denotes the fundamental solution
of the underlying differential operator and $v_{y}$ denotes the unit outer
normal at $y$ on the boundary. In certain structured matrix computations such
as those involving Cauchy and Cauchy-like matrices, $\kappa(x,y)$ is taken as
$\frac{1}{x-y}$ with $x,y\in\mathbb{C}$. In statistics and machine learning,
$\kappa(x,y)$ is often taken as the Gaussian kernel $e^{-|x-y|^{2}}$ or the
Laplace kernel $e^{-|x-y|}$. A major computational bottleneck in these
applications lies in the $O(n^{2})$ cost in storing the dense matrix $K$ and
performing operations such as matrix-vector multiplication.

To avoid the high cost in forming $K$ explicitly, hierarchically low-rank
matrix representations are used to approximate $K$. The hierarchical
representation is based on the fact that $K$ can be partitioned into blocks in
a hierarchical fashion and many blocks are numerically low rank. Low-rank
factors are computed for the blocks and are stored in the hierarchical
representation to replace original dense blocks. Two widely used hierarchical
representations are $\mathcal{H}$ and $\mathcal{H}^{2}$ [33, 9, 31, 30]. The
$\mathcal{H}$ matrix representation in general has a complexity of $O(n\log
n)$ in space while $\mathcal{H}^{2}$ has the optimal complexity of $O(n)$. The
matrix-vector multiplication can be computed in $O(n\log n)$ complexity for
$\mathcal{H}$ matrices and in $O(n)$ complexity for $\mathcal{H}^{2}$
matrices, which is much more efficient than directly multiplying $K$ by a
vector.

We review the mathematical description of hierarchical matrices in Section
2.1. Existing methods for constructing the hierarchical matrices are discussed
in Section 2.2.

###  2.1 Hierarchical matrix representations

In the following, we review the general algebraic framework of $\mathcal{H}$
and $\mathcal{H}^{2}$ matrices.

Given a dataset $X=\\{x_{i}\\}_{i=1}^{n}$ in $\mathbb{R}^{d}$, one builds a
tree structure by recursively partitioning $X$ spatially until no more than
$m=O(1)$ points are contained in each partitioned subset. The tree encodes the
subsets of $X$ generated by the adaptive partitions. Namely, the root node is
associated with $X$ and its children nodes are associated with subsets of $X$
created by the first partition. Inductively, each node is associated with a
non-empty subset of $X$. For node $i$, we denote by $X_{i}$ the subset
associated with that node. Hence $X_{\text{root}}=X$. An illustration is given
in Figure 2.

Figure 2: Adaptive partition of the dataset (left), label (at center) for each
subset (middle) and the associated partition tree (right).

The tree structure automatically yields a blockwise partition of the matrix
$K$. See Figure 3 for an example with data points lying inside an interval. We
define $K_{i,j}=[\kappa(x,y)]_{\begin{subarray}{c}x\in X_{i}\\\ y\in
X_{j}\end{subarray}}$. Block $K_{i,j}$ is approximated by a low-rank
factorization if $(i,j)$ satisfies an admissiblity condition that requires
$X_{i}$ and $X_{j}$ to be separated from each other to some extent. For
example, the pair $(i,j)$ is considered admissible if

$\text{diam}(X_{i})+\text{diam}(X_{j})\leq 2\tau|a_{i}-a_{j}|$ (1)

for some $\tau\in[0,0.7]$, where $a_{i}$ denotes the center of the box
associated with $X_{i}$ in the partition. The admissiblity condition in (1) is
used in [44, 15]. The parameter $\tau$, often called the _separation ratio_
[44], controls how separated the two subsets are. Smaller $\tau$ implies
better separation and the value $0.7$ can be replaced by other values less
than 1 (but not close to one). For any admissible $(i,j)$, the corresponding
submatrix $K_{i,j}$ is called an _admissible block_. An admissible block is
also referred as a _farfield_ block. Non-admissible blocks are referred as
_nearfield_ blocks. For a subset $X_{i}\subset X$, the union of all $X_{j}$
that are separated from $X_{i}$ in the sense of (1) is called the _farfield_
of $X_{i}$. Points that are not in the farfield of $X_{i}$ constitute the
_nearfield_ of $X_{i}$. In hierarchical matrices, each admissible block is
approximated by a low rank factorization,

$K_{i,j}\approx U_{i}B_{i,j}V_{j}^{T}$

where $U_{i}$, $V_{j}$ are column and row basis matrices, and $B_{i,j}$ is
called a _coupling matrix_. The maximum column size of all basis matrices
$U_{i}$ and $V_{j}$ is the approximation rank for $K$. A larger approximation
rank yields a more accurate approximation. For each node $i$, the _interaction
list_ of $i$ consists of nodes $j$ such that $X_{j}$ is well-separated from
$X_{i}$ but for the parent $p$ of $j$, $X_{p}$ is not well-separated from
$X_{i}$. The interaction list specifies the blockwise partition of the matrix
in which each admissible block has a low-rank approximation. For example, in
Figure 3, the interaction list of node $4$ is $\\{6,7\\}$; the interaction
list of node $10$ is $\\{8,12,13\\}$; The number of elements in any
interaction list is bounded by the so-called sparsity constant
$c_{\text{sp}}$, which is independent of $n$. The blockwise low-rank
representation yields the $\mathcal{H}$ matrix structure, which stores all
nearfield blocks and low-rank factors $U_{i},B_{i,j},V_{j}$. $\mathcal{H}$
matrices generally admit $O(n\log n)$ storage complexity. The more refined
$\mathcal{H}^{2}$ structure requires the basis matrices to be _nested_. The
nested bases property states that, if node $p$ has children
$c_{1},\dots,c_{k}$, then there exist _transfer matrices_ $R_{c_{i}}$,
$W_{c_{i}}$ such that

$U_{p}=\begin{bmatrix}U_{c_{1}}R_{c_{1}}\\\ \vdots\\\ U_{c_{k}}R_{c_{k}}\\\
\end{bmatrix},\quad V_{p}=\begin{bmatrix}V_{c_{1}}W_{c_{1}}\\\ \vdots\\\
V_{c_{k}}W_{c_{k}}\\\ \end{bmatrix}.$

See Figure 4 for an illustration. The sizes of the transfer matrices are equal
to $r$ if rank-$r$ factorizations are used for approximations to admissible
blocks. In practice, one uses $r=O(1)$ independent of $n$. Due to the nested
bases property, only basis matrices $U_{i}$, $V_{i}$ associated with leaf
nodes need to be stored in an $\mathcal{H}^{2}$ representation, in addition to
all transfer matrices (whose row and column sizes are $O(1)$) and nearfield
blocks. This results in $O(n)$ storage cost.

Figure 3: Hierarchical matrix structure for 1D problem (left to right): tree
(only top 4 levels are plotted), admissible blocks (colored) at levels 3,4,5,6
and all admissible blocks (colored) in the kernel matrix Figure 4: Admissible
block (10,12) and nested bases

Note that the above complexity estimates for storing hierarchical
representations assume that the low-rank factors have already been computed.
In practice, computing these hierarchical low-rank factors is usually the most
costly step, compared to applying the hierarchical representation to a vector.
Extensive research has focused on the efficient computation of the
hierarchical representation. We review several state-of-the-art methods in
Section 2.2.

###  2.2 General-purpose methods

For an $n$-by-$n$ kernel matrix $K$ associated with a general kernel function
$\kappa(x,y)$, the hierarchical matrix representation can be computed in
linear or quasilinear time with a variety of techniques, including
interpolation [33, 30], adaptive cross approximation (ACA) [5, 6], hybrid
cross approximation (HCA) [8], SMASH [15], etc. These methods work for general
kernels (for example, non-symmetric, non-translationally-invariant). For
special kernel functions, such as fundamental solutions of certain elliptic
partial differential equations (PDEs), layer potentials in integral equations,
Coulomb interactions in electrostatics, analytic methods such as the fast
multipole method [40, 41, 28, 44, 19] and its variants [1, 45] can be used to
construct a hierarchical matrix representation efficiently.

When constructing the hierarchical matrix with _nested bases_ , one needs to
compute basis matrices $U_{i}$, $V_{i}$ for submatrices $K_{X_{i}Y_{i}}$ that
account for the interaction between $X_{i}$ with $O(1)$ points and its entire
farfield $Y_{i}$ with $O(n)$ points (see Figure 1(a)). In order to achieve
$O(n)$ optimal complexity for building the hierarchical representation, each
basis matrix $U_{i}$ (as well as $V_{i}$) must be computed with $O(1)$
complexity. This requires that the matrix $K_{X_{i}Y_{i}}$, which is
$O(1)$-by-$O(n)$ in size, must not be formed. For general kernels, a widely
used technique to construct a column basis matrix $U_{i}$ in $O(1)$ complexity
is based on interpolation. For $K_{X_{i}Y_{i}}$, interpolating the kernel
function $\kappa(x,y)$ at $r$ nodes $Q=\\{q_{1},\dots,q_{r}\\}$:

$\kappa(x,y)\approx\sum_{i=1}^{r}\kappa(q_{i},y)L_{i}(x)$

yields a rank-$r$ approximation:

$K_{X_{i}Y_{i}}\approx U_{i}K_{QY_{i}},$ (2)

where $L_{i}$ is the Lagrange polynomial corresponding to node $q_{i}$ and
$U_{i}=[L_{k}(x)]_{\begin{subarray}{c}x\in X_{i}\\\ k=1:r\end{subarray}}$. Due
to its generality and efficiency in computing $U_{i}$, this interpolation
method is used in a number of general-purpose hierarchical matrix algorithms,
e.g. [33, 9, 8, 15].

Another way of finding column basis matrix for $K_{X_{i}Y_{i}}$ is through
subset selection [6, 25]. The column basis matrix is chosen as the submatrix
corresponding to a judiciously chosen $O(1)$ subset $Y_{i}^{*}$ from $Y_{i}$.
Similar to interpolation, subset selection can be applied to general kernel
functions. The reference [6] presents an efficient hierarchical scheme to
select representor sets for all nodes using two steps: top-down and bottom-up.
The cost of the algorithm in [6] is dominated by computing farfield
representor sets in the top-down step. Computing representor sets for all
nodes $i$ leads to $O(n\log n)$ complexity for a balanced tree with $O(\log
n)$ levels. The total complexity of the resulting hierarchical matrix
construction is $O(n\log n)$, instead of the optimal $O(n)$ complexity
achieved by interpolation-based methods. However, compared to using
interpolation nodes (which are generally _outside_ the given dataset),
selecting subsets directly from the dataset is more memory-efficient for low-
rank approximation. See Section 6 for a detailed discussion.

##  3 Fast Hierarchical Data Reduction (HiDR)

To facilitate the fast construction of hierarchically low rank
representations, we propose an efficient preprocessing scheme to reduce the
tree-structured data so that each node in the partition tree induces $O(1)$
cost in the subsequent hierarchical matrix construction process. Specifically,
let $X_{i}$ be the set of points corresponding to node $i$ and $Y_{i}$ be the
farfield of $X_{i}$. The data reduction aims to find representor sets
$X_{i}^{*}\subset X_{i}$ with $O(1)$ points and $Y_{i}^{*}\subset Y_{i}$ with
$O(1)$ points for each node $i$. Note that a naive data reduction for $Y_{i}$
with $O(n)$ points into a subset of evenly spaced points as shown in Figure 1
will lead to $O(n)$ computational complexity. The cost can be reduced to
$O(1)$ with a carefully designed hierarchical procedure presented in Section
3.1.

We present the hierarchical data reduction algorithm in Section 3.1 and verify
that it scales linearly with the size of the data in Section 3.2. Several
algorithms for performing data reduction are discussed in Section 3.3.

###  3.1 Linear complexity hierarchical data reduction

The fast HiDR consists of two traversals of the tree: bottom-up and top-down.
The $O(1)$ representor set $X_{i}^{*}$ for $X_{i}$ is computed in the bottom-
up pass and the $O(1)$ farfield representor set $Y_{i}^{*}$ for $Y_{i}$ is
computed in the top-down pass. A building block for the hierarchical scheme is
a DataReduct subroutine that takes the form:

${\tt DataReduct}(X,k)\to X^{*},$

where $X^{*}$ is a subset (representor set) of the input $X$ and $k$ is a
parameter that specifies the size of $X^{*}$. There are several options for
the subroutine DataReduct to obtain $X^{*}$ from $X$. A detailed discussion is
presented in Section 3.3. Given a set of points, the subroutine selects a
subset of evenly spaced points whose size is bounded by a prescribed constant
and scales linearly with the size of the input data. The HiDR is designed such
that the input dataset for DataReduct is always $O(1)$ in size.

Figure 5: Bottom-up pass of the tree: representor sets $X_{i}^{*}$ are first
computed for nodes $i$ at the deepest level and then upper levels. Nodes with
$X_{i}^{*}$ computed are wiped out. The three trees from left to right
correspond to the 1st, 3rd, 5th configurations in Figure 6.

In the bottom-up sweep, starting from leaf nodes $i$, each $X_{i}$ contains
$O(1)$ points and thus the data reduction from $X_{i}$ to $X_{i}^{*}$ induces
$O(1)$ cost only. After children nodes have been processed, we define for each
parent $p$ an intermediate set $S_{p}$ as the union of reduced sets
$X_{i}^{*}$ for all its children $i$. The representor set $X_{p}^{*}$ is
obtained by applying data reduction to the intermediate set $S_{p}$. Since
each parent has at most $C$ children ($C=2,4,8$ for binary tree, quadtree,
octree, respectively) and each $X_{i}^{*}$ is $O(1)$ in size, the intermediate
set $S_{p}$ is always $O(1)$ in size. Thus the cost of computing the
representor set for $p$ is always $O(1)$. Recursively, all non-root nodes can
be processed in $O(1)$ complexity in the bottom-up procedure (see Figure 5).
Figure 6 shows the data reduction to generate $X_{i}^{*}$ from the lowest
level to upper levels of the tree constructed in Figure 2 where each
$X_{i}^{*}$ contains 2 points.

Figure 6: Fast hierarchical data reduction: bottom-up pass. Each $X_{i}$ is
reduced to $X_{i}^{*}$ such that $X_{i}^{*}$ contains at most 2 points. The
1st, 3rd, 5th arrows correspond to the subroutine DataReduct for the boxes at
the bottom level of the tree, where points in red are output of DataReduct.
The 2nd, 4th arrows correspond to merging children boxes and going up in the
tree.

In the top-down sweep, starting from the top nodes $i$ in the tree with non-
empty interaction list, we define the intermediate set $T_{i}$ as the union of
$Y_{p}$ ($p$ denotes the parent of $i$) and $X_{j}^{*}$ for all nodes $j$ in
the interaction list of $i$. The representor set $Y_{i}^{*}$ for the farfield
$Y_{i}$ of $X_{i}$ is computed by applying data reduction to the intermediate
set $T_{i}$. Since the cardinality of the interaction list of $i$ is bounded
by the sparsity constant $c_{\text{sp}}$, which is $O(1)$, and each
$X_{j}^{*}$ is already $O(1)$ in size, we see that $T_{i}$ only contains
$O(1)$ points and consequently computing $Y_{i}^{*}$ has $O(1)$ cost. Once
parent nodes have been processed, we define the intermediate set $T_{i}$ for
each child $i$ as the union of all $X_{j}^{*}$ from its interaction list and
$Y_{p}^{*}$ from its parent $p$. From this definition, $T_{i}$ is also $O(1)$
in size. Similar to the above, the representor set $Y_{i}^{*}$ is then
obtained by applying data reduction to $T_{i}$. Recursively, each reduced
farfield representation $Y_{i}^{*}$ can be computed in $O(1)$ complexity via
the top-down procedure. Figure 7 shows the top-down data reduction for
farfield to generate $Y_{i}^{*}$ (circled), as $i$ goes from a node near the
root node to a leaf node. In Figure 7, each $Y_{i}^{*}$ contains at most 4
points only.

Figure 7: Fast hierarchical data reduction: top-down pass. For each
$X_{i}^{*}$ in an orange box, the blue region contains the entire farfield,
where each representor set $Y_{i}^{*}$ (circled) contains at most 4 points.

The HiDR algorithm is summarized in Algorithm 1. In practice, the cost of HiDR
is lower than the hierarchical low-rank compression. HiDR reduces the
computational cost of the subsequent hierarchical matrix construction.

Algorithm 1 Hierarchical data reduction (HiDR)

_Input:_ The adaptive partition tree $\mathcal{T}$ for dataset $X$, the
collection of subsets $X_{i}$ for all leaf nodes $i$, prescribed maximum size
$r_{1}$ of $X_{i}^{*}$ and maximum size $r_{2}$ of $Y_{i}^{*}$
_Output:_ Reduced representations $X_{i}^{*}$ and $Y_{i}^{*}$ for all nodes
$i$

1: for all $i\in\mathcal{T}$ do

2: $Y_{i}^{*}=\emptyset$, $S_{i}=\emptyset$

3: if $i$ is a leaf node then

4: $S_{i}=X_{i}$

5: end if

6: end for

7: for each level (from bottom to top) do

8: for all $i$ at this level do

9: if $i$ is a parent then

10: $S_{i}=\bigcup\limits_{c\in\text{ch}(i)}X_{c}^{*}$ with ch($i$) the set of
children of node $i$

11: end if

12: $X_{i}^{*}$ = DataReduct($S_{i}$,$r_{1}$)

13: end for

14: end for

15: for each level (from top to bottom) do

16: for all $i$ at this level do

17: if $i$ has non-empty farfield then

18: $T_{i}=Y_{p}^{*}\bigcup X_{j}^{*}$ over all $j$ in the interaction list of
$i$ ($p$ denotes the parent of $i$)

19: $Y_{i}^{*}$=DataReduct($T_{i}$,$r_{2}$)

20: end if

21: end for

22: end for

23: return $X_{i}^{*}$, $Y_{i}^{*}$ for all nodes $i\in\mathcal{T}$

###  3.2 Complexity analysis

In this section, we show that Algorithm 1 (HiDR) has linear complexity with
respect to the number of points in $X$.

######  Theorem 3.1.

For the given dataset $X$ that contains $n$ points, let $\mathcal{T}$ be a
partition tree for $X$, in which each leaf node corresponds to a subset of $X$
with $O(1)$ points and the sparsity constant $c_{\text{sp}}=O(1)$. Then the
complexity of Algorithm 1 is $O(n)$.

###### Proof.

We first analyze the complexity in the bottom-up pass. Since for each leaf
node $i$, $X_{i}$ contains $O(1)$ points, the DataReduct for $X_{i}$ has $(1)$
cost. We next prove by induction that the cost for each nonleaf node is also
$O(1)$. For a nonleaf node $i$, assume that for each child $c$, $X_{c}^{*}$
contains $O(1)$ points. Then the intermediate set $S_{i}$ in Line 10 of
Algorithm 1 contains $O(1)$ points because $i$ has $O(1)$ children and
$X_{c}^{*}$ contains $O(1)$ points for each child $c$. Consequently, the cost
of DataReduct applied to $S_{i}$ is $O(1)$. This shows that the cost to obtain
$X_{i}^{*}$ is also $O(1)$. From the induction, we conclude that the cost for
_each_ node is $O(1)$ and each reduced subset $X_{i}^{*}$ has $O(1)$ points.

Now we analyze the complexity for the top-down pass. First note that each
$Y_{i}^{*}$ has at most $O(1)$ points according to the construction.
Consequently, each $T_{i}$ in Line 18 of Algorithm 1 contains at most $O(1)$
points because the interaction list of $i$ contains at most
$c_{\text{sp}}=O(1)$ nodes and each $X_{j}^{*}$ has $O(1)$ points. This
implies that the DataReduct of $T_{i}$ to generate $Y_{i}^{*}$ has $O(1)$
complexity.

Overall, we see that for each node $i$, the associated total cost to compute
$X_{i}^{*}$ and $Y_{i}^{*}$ is $O(1)$. Since there are $O(n)$ nodes in the
tree, the total cost for all nodes is $O(n)$. This completes the proof of the
theorem. ∎

###  3.3 Data reduction methods

In this section, we provide several algorithms for performing data reduction.
The goal is to select a subset of $X_{i}$ and of $Y_{i}$ such that the
selected subsets preserve the geometry of $X_{i}$ and $Y_{i}$. For an input
set $X$, the subroutine takes the simple form: DataReduct($X$,$k$) $\to$
$X^{*}$ with $X^{*}\subset X$ the selected subset whose size is controlled by
the parameter $k=O(1)$. Note that DataReduct only depends on the input data
and is independent of any kernel function. It has been shown in [17] that the
choice of the subset is essential for the accuracy and robustness of the low
rank approximation. According to the results in [17, 16], a subset evenly
distributed over the containing set can offer an improved approximation
robustness and accuracy over one that is not. The methods below can be used to
generate such a subset efficiently. Many of them rely on a reference set with
good uniformity, such as a uniform tensor grid. In addition to a tensor grid,
it was shown recently that deep neural networks can be used to generate
distributions with good uniformity [10].

We briefly review some of the existing data reduction methods below. An
empirical comparison of these methods for the low-rank approximation of kernel
matrices is presented in Section 5.2.

##### Farthest point sampling.

Given $X$ and a target size $k$ for the reduced subset $S$ of $X$, farthest
point sampling (FPS) constructs $S$ in a sequential manner. $S$ is initialized
with one point only. Then FPS searches for a point in $X\backslash S$ that is
farthest from $S$ and adds the point to $S$. This procedure is repeated until
$S$ reaches size $k$. FPS generates evenly distributed subsets and has been
widely used in computational geometry [24, 39, 42]. FPS was recently proposed
for computing low rank approximations [16].

##### Volume-based data reduction.

Volume based data reductions choose a subset $S$ of $X$ via a reference grid
$Q$ with $O(1)$ points inside the computational domain. For example, $Q$ can
be chosen as a tensor grid (cf. [6]) inside the rectangular domain that
encloses the data. $S$ is chosen to be the collection of points in $X$ that
are closest to each point in $Q$.

##### Surface-based data reduction.

Following the same idea as the volume-based method, we can also use a
reference set $Q$ based on surfaces constructed from the given data $X$. In
the surface-based method, we define $Q$ as the union of points distributed on
surfaces near the boundary of the computational domain. For example, we can
construct ellipsoids centered at the center of a rectangular domain that
encloses $X$. The principal semi-axes of the ellipsoids are chosen to be equal
to $\gamma$ times the width of the rectangular domain in each dimension, where
$\gamma>0$ is a hyperparameter. In Section 5.2, we use three ellipsoids with
$\gamma=0.3,0.6,1.2$.

##### Anchor net method.

The anchor net method [17] is a newly proposed subset selection method based
on approximating the geometry of the given dataset with low discrepancy
subsets. For low-rank approximation, the anchor net method is shown to achieve
a good time-accuracy trade-off in practice and is particularly efficient for
high-dimensional data (cf. [17, 16]).

##  4 Data-driven Hierarchical Matrix Construction

In this section, we first show how to extract a low-rank factorization
instantly for an admissible block after the prepossessing procedure HiDR in
Section 4.1 and then present an algorithm with $O(n)$ complexity (Algorithm 2)
for constructing an $\mathcal{H}^{2}$ matrix representation in Section 4.2. We
analyze the computational complexity of Algorithm 2 in Section 4.3. The
proposed method enjoys the following features:

  1. a)

black-box general-purpose (kernel independent) construction of the
hierarchical low rank format;

  2. b)

optimal $O(n)$ complexity, where $n$ is the number of points in the given
dataset;

  3. c)

better efficiency for data from complex geometry compared to general-purpose
approaches as well as specialized kernel-dependent techniques.

Algorithm 2 Data-driven (DD) HiDR-based $\mathcal{H}^{2}$ matrix construction

_Input:_ Dataset $X$, kernel function $\kappa(x,y)$, approximation tolerance
$\epsilon$, separation ratio $\tau$, maximum number of points $q$ for a leaf
node
_Output:_ $\mathcal{H}^{2}$ matrix representation

1: Apply adaptive partitioning to $X$ to generate the partition tree
$\mathcal{T}$ with at most $q$ points for each leaf node and obtain subsets
$X_{i}$ for all _leaf_ nodes $i$

2: Determine approximation parameters $r_{1},r_{2}$ from $\epsilon$

3: Apply HiDR in Algorithm 1 with approximation parameters $r_{1},r_{2}$ to
obtain $O(1)$ representor sets $X_{i}^{*}$ and $Y_{i}^{*}$ for all nodes
$i\in\mathcal{T}$

4: for all $i\in\mathcal{T}$ do

5: Define $\bar{X}_{i}^{\text{(row)}}=\bar{X}_{i}^{\text{(col)}}=\emptyset$

6: if $i$ is a leaf node then

7: define $\bar{X}_{i}^{\text{(row)}}=\bar{X}_{i}^{\text{(col)}}=X_{i}$

8: end if

9: end for

10: for all non-root $i\in\mathcal{T}$ from bottom level to top level do

11: Apply getBasis($K_{\bar{X}_{i}^{\text{(row)}}Y_{i}^{*}}$) to obtain
$U_{i},\hat{X}_{i}^{\text{(row)}}$ and
getBasis($K_{Y_{i}^{*}\bar{X}_{i}^{\text{(col)}}}^{T}$) to obtain
$V_{i},\hat{X}_{i}^{\text{(col)}}$

12: Update
$\bar{X}_{p}^{\text{(row)}}=\bar{X}_{p}^{\text{(row)}}\cup\hat{X}_{i}^{\text{(row)}}$
and
$\bar{X}_{p}^{\text{(col)}}=\bar{X}_{p}^{\text{(col)}}\cup\hat{X}_{i}^{\text{(col)}}$,
where $p$ is the parent of $i$

13: end for

14: Define the $\mathcal{H}^{2}$ column and row basis matrices:
$U=\\{U_{i}\\}_{\text{leaf }i}$, $V=\\{V_{i}\\}_{\text{leaf }i}$, transfer
matrices $R=\\{R_{i}\\},W=\\{W_{i}\\}$:
$\displaystyle\begin{bmatrix}R_{c_{1}}\\\ \vdots\\\
R_{c_{k}}\end{bmatrix}=U_{i},\quad\begin{bmatrix}W_{c_{1}}\\\ \vdots\\\
W_{c_{k}}\end{bmatrix}=V_{i},\quad\text{if $i$ has children
$c_{1},\dots,c_{k}$},$ coupling matrices $B=\\{B_{i,j}\\}$:
$B_{i,j}=\begin{cases}K_{\hat{X}_{i}^{\text{(row)}}\hat{X}_{j}^{\text{(col)}}},&\text{if
$(i,j)$ is admissible}\\\ K_{X_{i}X_{j}}&\text{otherwise}\end{cases}$

15: return $\mathcal{H}^{2}$ representation: $U,V,R,W,B$

###  4.1 Approximating the entire farfield $K_{X_{i}Y_{i}}$

In this section, we show how to derive a low-rank approximation for the entire
farfield $K_{X_{i}Y_{i}}$ based on the representor sets $X_{i}^{*},Y_{i}^{*}$
returned by Algorithm 1.

When computing an approximate column basis for $K_{X_{i}Y_{i}}$, $X_{i}$
contains $O(1)$ points and its entire farfield $Y_{i}$ contains $O(n)$ points
for the leaf node $i$ (see Figure 1). In order for the entire algorithm to
have linear complexity, the column basis matrix of $K_{X_{i}Y_{i}}$ must be
computed in $O(1)$ complexity.

We first apply strong rank-revealing QR (SRRQR) factorization [29] to the
submatrix $K_{X_{i}Y_{i}^{*}}$:

$K_{X_{i}Y_{i}^{*}}=P\begin{bmatrix}I\\\
G\end{bmatrix}K_{\hat{X}_{i}Y_{i}^{*}},$ (3)

where $P$ is a permutation matrix, $||G||_{\max}$ is bounded by a prescribed
constant, and $\hat{X}_{i}$ is a subset of $X_{i}$ with $O(1)$ points. Then
the column basis matrix is chosen as $U_{i}=P[I;G^{T}]^{T}$ and the low-rank
approximation for the entire farfield reads

$K_{X_{i}Y_{i}}\approx P\begin{bmatrix}I\\\
G\end{bmatrix}K_{\hat{X}_{i}Y_{i}}.$ (4)

For notational convenience, we denote the procedure in (3)-(4) for computing
an approximate column basis for $K_{X_{i}Y_{i}}$ by

$\text{getBasis}(K_{X_{i}Y_{i}})=(U_{i},\hat{X}_{i}),$ (5)

where $U_{i}:=P[I;G^{T}]^{T}$ is the computed column basis and
$\hat{X}_{i}\subset X_{i}$. Note that the kernel matrix $K_{X_{i}Y_{i}}$ is
never formed because the input of “getBasis” is the kernel function and the
subsets $X_{i},Y_{i}^{*}$.

The cost to obtain $U_{i}$ and $\hat{X}_{i}$ from (5) is $O(1)$ as the matrix
$K_{X_{i}Y_{i}^{*}}$ is $O(1)$-by-$O(1)$. Also notice that

$\text{card}(\hat{X}_{i})\leq\text{rank}(K_{X_{i}Y_{i}^{*}})\leq\text{card}(Y_{i}^{*})=O(1).$
(6)

Thus, $\hat{X}_{i}$ always contains $O(1)$ points.

As we shall see in Section 6, the column basis $U_{i}$ derived from
$K_{X_{i}Y_{i}^{*}}$ can yield better accuracy than analytic methods such as
interpolation.

###  4.2 Computing hierarchical matrices with nested bases using HiDR

Hierarchical matrices with nested bases, e.g. $\mathcal{H}^{2}$ matrices, can
offer optimal $O(n)$ complexity in time and space when approximating an
$n$-by-$n$ kernel matrix. A black-box hierarchical matrix construction
proposed in [15] works for general kernel functions and allows for arbitrary
low-rank compression techniques. In this section, we show that the
hierarchical data reduction (HiDR) can be incorporated naturally into the
construction of $\mathcal{H}^{2}$ matrix representations via the general
framework proposed in SMASH [15]. SMASH employs a bottom-up procedure that
recursively applies rank-revealing factorization to the initial basis matrix
(with $O(1)$ entries) for each node in the tree. In [15], the initial basis
matrices are constructed via either interpolation or analytic expansion of the
kernel function. In this section, we leverage representor sets produced by
HiDR in Algorithm 1 to construct the initial basis matrices.

The full data-driven construction is presented in Algorithm 2. The algorithm
automatically determines the approximation parameters $r_{1},r_{2}$ for HiDR
in Algorithm 1 according to the approximation tolerance $\epsilon$ prescribed
by the user. The idea here is to apply the low-rank approximation in (4) to
the artificial kernel matrix $K_{Z_{1}Z_{2}}$ where
$Z_{1},Z_{2}\subset\mathbb{R}^{d}$ are well-separated subsets (in the sense of
(1)) of $O(1)$ random points. The parameters $r_{1},r_{2}$ are chosen
adaptively by increasing from $r_{1}=r_{2}=1$ to a point such that the
approximation error to $K_{Z_{1}Z_{2}}$ is smaller than $10^{-2}\epsilon$.
More sophisticated techniques like a posteriori error estimation (cf. [12, 13,
14, 11]) can also be studied to estimate the approximation error. After the
parameters are determined, Algorithm 1 first applies hierarchical data
reduction to $X$ associated with tree $\mathcal{T}$ to obtain representor
sets. Then the hierarchical matrix representation can be computed rapidly by
following the SMASH $\mathcal{H}^{2}$ construction and using
$K_{X_{i}Y_{i}^{*}}$ as the initial basis matrix for each leaf node $i$.

We perform numerical experiments in Section 6 to demonstrate that the new
data-driven method improves the matrix approximation accuracy of
interpolation-based SMASH algorithm [15]. Moreover, the cost of hierarchical
data reduction is smaller than the subsequent hierarchical matrix compression.

###  4.3 Complexity analysis

######  Theorem 4.1.

For the given dataset $X$ with $n$ points, let $\mathcal{T}$ be a partition
tree for $X$, in which each leaf node corresponds to a subset of $X$ with
$O(1)$ points and the sparsity constant $c_{\text{sp}}=O(1)$. Then the
complexity of Algorithm 2 is $O(n)$.

###### Proof.

Algorithm 2 follows the $\mathcal{H}^{2}$ construction in SMASH [15] with an
additional hierarchical data reduction (HiDR) in Line 3. According to Theorem
3.1, HiDR has $O(n)$ complexity. Thus to prove the $O(n)$ complexity of
Algorithm 2, it suffices to show that the cost per node is $O(1)$ in Lines
11–12.

We first show that $\bar{X}_{i}^{\text{(row)}}$ and
$\bar{X}_{i}^{\text{(col)}}$ contain at most $O(1)$ points. If $i$ is a leaf
node, then according to the definition in Line 7,
$\bar{X}_{i}^{\text{(row)}}=\bar{X}_{i}^{\text{(col)}}=X_{i}$ contain $O(1)$
points. If $i$ is a parent node, then after all children of $i$ have been
updated, $\bar{X}_{i}^{\text{(row)}}$ and $\bar{X}_{i}^{\text{(col)}}$ can be
written as

$\bar{X}_{i}^{\text{(row)}}=\bigcup\limits_{c\text{ is a child of
}i}\hat{X}_{c}^{\text{(row)}},\quad\bar{X}_{i}^{\text{(col)}}=\bigcup\limits_{c\text{
is a child of }i}\hat{X}_{c}^{\text{(col)}}.$

Since the number of children for every node is bounded from above by a
constant, and every subset $\hat{X}_{c}$ contains $O(1)$ points according to
(6), we see that $\bar{X}_{i}^{\text{(row)}}$ and $\bar{X}_{i}^{\text{(col)}}$
contain $O(1)$ points. This implies that the complexity in Line 12 is at most
$O(1)$ for all $i$.

Next we analyze the complexity in Line 11. Since $Y_{i}^{*}$ contains $O(1)$
points only, it follows that in Line 11, the input matrices
$K_{\bar{X}_{i}^{\text{(row)}}Y_{i}^{*}}$ and
$K_{Y_{i}^{*}\bar{X}_{i}^{\text{(col)}}}^{T}$ have $O(1)$ rows and columns.
Consequently, performing “getBasis” in Line 11 only takes $O(1)$ time.

Now we conclude that the total complexity in Lines 11–12 is $O(1)$. Therefore,
the total complexity for Algorithm 2 is $O(n)$. ∎

##  5 Data Reduction and Low Rank Approximation

In this section, we investigate different data reduction techniques for low
rank approximation. Section 5.1 presents an example to reveal a drawback of
methods that rely on points outside the given dataset (such as interpolation
nodes or random points) for computing the low-rank approximation to kernel
matrices. In Section 5.2, we compare the performance of the data reduction
methods of Section 3.3 for low-rank approximation.

###  5.1 Drawback of using points outside the given data for irregular
datasets

Interpolation nodes, random points, or in general points outside the given
dataset are commonly used in low-rank approximation to obtain approximate
column basis efficiently without forming the original kernel matrix. Since
these points are created artificially (not part of the given data), we call
these points _virtual points_. For given data $X$ and $Y$, virtual points are
constructed in rectangular domains $\Omega_{X}$ and $\Omega_{Y}$ that cover
$X$ and $Y$, respectively, and the kernel matrix associated with these virtual
points is used to compute the low-rank approximation.

One issue of of using virtual points is that it may lead to an incorrect
approximation to the kernel matrix. This is because the virtual points lie
outside the original data, and the kernel matrix involving these virtual
points may have a very different spectrum from that of the kernel matrix to be
approximated. Thus methods using virtual points may not be robust for
approximating general low-rank kernel matrices. To illustrate the issue, we
use the “nJ” dataset as illustrated in Figure 8(a), where $X$ contains 120
points in $\Omega_{X}=[-2.5,2.5]\times[0,3.75]$ and $Y$ contains 150 points in
$\Omega_{Y}=[0,5.13]\times[16.25,23.84]$. We consider the smooth kernel
function

$\kappa(x,y)=\sqrt{1+100|x-y+a|^{2}}$

with $a=[0,20]^{T}$. The same issue also arises for other kernels such as
Gaussians. The corresponding kernel matrix $K_{XY}$ (120-by-150) has rapidly
decaying singular values as shown in Figure 9 (dashed line). Consequently,
$K_{XY}$ can be approximated very well by a low-rank matrix.

Now consider $K_{X_{1}Y_{1}}$, with $X_{1}$ being 2000 random points selected
in $\Omega_{X}$ and $Y_{1}$ being 2000 random points selected in $\Omega_{Y}$.
The singular values of $K_{X_{1}Y_{1}}$ are an approximation to the continuous
singular values of the problem in (7). The computed singular values of
$K_{X_{1}Y_{1}}$ are also plotted in Figure 9. It can be seen that the
singular values of $K_{X_{1}Y_{1}}$ do not decay rapidly, compared to
$K_{XY}$.

Now consider $K_{X_{2}Y_{2}}$, where $X_{2}$ and $Y_{2}$ are $10\times 10$
Chebyshev points in $\Omega_{X}$ and $\Omega_{Y}$, respectively. Like
$K_{X_{1}Y_{1}}$, this matrix does not have singular values that decay as
rapidly as those of $K_{XY}$, and therefore an algebraic compression of
$K_{X_{2}Y_{2}}$ will not be an effective approximation for $K_{XY}$.

(a) Dataset $X$, $Y$

(b) 2000 random points in $\Omega_{X}$, $\Omega_{Y}$

(c) $10\times 10$ Chebyshev points in $\Omega_{X}$, $\Omega_{Y}$

Figure 8: Dataset $X,Y$ and different types of “virtual points” within
$\Omega_{X}$ and $\Omega_{Y}$. Figure 9: Distinct singular value patterns for
kernel matrices with original data $X\times Y$ and virtual points $X_{i}\times
Y_{i}$ ($i=1,2$) of Figure 8: singular values of $K_{XY}$ (120-by-150),
$K_{X_{1}Y_{1}}$ (2000-by-2000), $K_{X_{2}Y_{2}}$ (100-by-100). For
$K_{X_{1}Y_{1}}$, the largest 150 singular values are plotted.

##### Mathematical explanation.

Employing a continuous treatment of the matrix approximation problem ignores
the geometry of the discrete dataset. This can be problematic in general as
the continuous problem may have entirely different spectral properties
compared to the matrix. It is even possible that the kernel function is
_undefined_ at virtual points. For the model problem in Figure 8, the matrix
$K_{X_{1}Y_{1}}$ with virtual points $X_{1},Y_{1}$ is related to the following
integral operator:

$T:L^{2}(\Omega_{X}):\to
L^{2}(\Omega_{Y}),\quad(Tf)(x):=\int_{\Omega_{X}}\kappa(x,y)f(y)dy.$ (7)

The singular values of $K_{X_{1}Y_{1}}$ (with $X_{1}$ and $Y_{1}$ chosen as
described above) approximate the singular values (up to a scaling constant) of
the integral operator [2, 35, 43, 38, 3, 18]. These singular values do not
decay rapidly like those of $K_{XY}$. In essence, we see that virtual points
methods treat the matrix approximation as a _continuous_ problem and thus
ignore the geometry of the discrete data. When the continuous problem differs
substantially from the original _discrete_ problem (kernel matrix
approximation), the performance of methods that utilize virtual points can be
very unsatisfactory.

###  5.2 Comparison of data reduction methods

In this section, we perform experiments to compare the performance of the four
data reduction methods in Section 3.3: farthest point sampling (‘FPS’),
volume-based reduction (‘Volume’), surface-based reduction (‘Surface’), anchor
net method (‘AnchorNet’). These methods operate on the dataset and do not
require any kernel function.

##### Experiment setup

We consider low-rank approximation to the kernel matrix $K_{XY}$, where $X$
(198 points) and $Y$ (1577 points) are well-separated subsets from a dinosaur
manifold as shown in Figure 10. The diameter of $X$ is 58.21, and the distance
between $X$ and $Y$ is 21.275. We test three different kernel matrices
$K_{XY}$ corresponding to the kernel functions below:

$\frac{1}{|x-y|},\quad e^{-\frac{|x-y|^{2}}{900}},\quad|x-y|^{11}.$

To obtain the low-rank approximation, we first perform data reduction for $Y$,
and then build the factorization as described in Section 4.1 using (3) and
(4). The low-rank approximation error is measured by the relative matrix
approximation error in the 2-norm.

The error plots for three different kernels are shown in Figure 11. For each
plot, the horizontal axis denotes the number of points selected by the data
reduction method, namely the size of the subset $Y^{*}\subset Y$. Each curve
shows how the low-rank approximation error for a specific data reduction
method decays as we increase the size of $Y^{*}$. We see that ‘Volume’ and
‘AnchorNet’ offer the best performance and are almost indistinguishable from
each other in performance across all three kernels. ‘Surface’ achieves similar
performance for the first two kernels but is slightly worse than ‘Volume’ and
‘AnchorNet’ for the third kernel. The farthest point sampling ‘FPS’ performs
well but is not as accurate as the other three methods for the same number of
selected points for all kernels tested.

Figure 10: Section 5.2 dataset: $X$ and $Y$ well-separated

Figure 11: Comparison of data reduction methods for low-rank approximation to
$K_{XY}$ with dataset $X\times Y$ in Figure 10 and different kernels $k(x,y)$
on top of each plot

##  6 Numerical experiments

We present a series of numerical experiments in this section to illustrate the
performance of the proposed data-driven construction in Algorithm 2. The code
for the algorithm is available on GitHub111https://github.com/scalable-
matrix/H2Pack/tree/sample-pt-algo. The performance of the proposed data-driven
hierarchical matrix construction is shown in Section 6.1, inlcuding linear
scaling, generality for various kinds of kernels, and the efficiency of
hierarchical data reduction for varying kernel parameters. Comparison to the
state-of-the-art special-purpose methods for the Coulomb kernel is presented
in Section 6.2. Comparison to the widely used general-purpose method
(interpolation) for various kernels is presented in Section 6.3. For the data-
driven hierarchical construction, volume-based data reduction is used. For
experiments in Section 6.1 and Section 6.2, we use one compute node on the
Georgia Tech PACE-Hive cluster. This node has two sockets and 192 GB DDR4
memory. Each socket has an Intel Xeon Gold 6226 12-core processor.

The approximation error is measured by the relative matrix-vector product
error: $\frac{||Kz-\tilde{K}z||}{||z||}$, where $\tilde{K}$ denotes the
hierarchical approximation to the kernel matrix $K$ and $z$ is a standard
normal random vector. $||\cdot||$ denotes the 2-norm.

###  6.1 Data-driven construction: scaling, generality, once-for-all HiDR

This section has three objectives: (1) test the complexity of the proposed
data-driven approach in Section 6.1.1, including the hierarchical data
reduction (HiDR) and the resulting hierarchical matrix construction; (2)
illustrate the generality of the data-driven approach by testing different
kernels in Section 6.1.2; (3) apply HiDR once and use the representor sets to
construct hierarchical matrices for various types of kernels, including
Gaussian kernels with different bandwidths in Section 6.1.3.

####  6.1.1 Scaling test for different datasets

##### Datasets

Three datasets are used: cube, 3-sphere, Dino. The “cube” dataset contains
random samples from the uniform distribution in the unit cube $[0,1]^{3}$. The
“Dino” dataset is used in [15] and is illustrated in Figure 12. It consists of
points distributed on a dinosaur-shaped surface in three dimensions. The
“3-sphere” dataset (see Figure 12) consists of random points distributed on
the surface of three intersecting unit spheres whose centers form an
equilateral triangle with side length close to 1. Roughly the same number of
points is sampled from each sphere. Let $n$ denote the number of points in the
dataset. For the first three synthetic datasets, we test for $n$ from $10^{5}$
to $1.6\times 10^{7}$. For the Dino dataset, since the size of the original
data is fixed, we sample $n$ points randomly and vary $n$ from $10^{4}$ to
$1.5\times 10^{5}$.

Figure 12: 3-sphere (left) and Dino (right) datasets

The kernel function is chosen to be the Coulomb kernel $\frac{1}{|x-y|}$ and
the approximation error is $10^{-6}$ for each test. Figure 13 shows the
timings for hierarchical data reduction (‘HiDR’), $\mathcal{H}^{2}$ matrix
construction (‘build’), the resulting matrix-vector multiplication (‘matvec’)
with respect to $n$, respectively. All timings scale linearly with $n$ and the
hierarchical data reduction (‘HiDR’) has a much lower cost than the subsequent
hierarchical matrix construction (‘build’). The low cost of data reduction and
the kernel independence make the data-driven approach suitable for the case
when the kernel matrix changes frequently due to changes in the data or kernel
function.

Figure 13: Section 6.1.1 experiment: Timings of HiDR, $\mathcal{H}^{2}$ build,
matvec for $n$-by-$n$ Coulomb kernel matrices with three datasets in
$\mathbb{R}^{3}$

####  6.1.2 Scaling test for different kernels

In this section, we test the proposed data-driven algorithm for the kernel
functions in Table 1. We show that the general data-driven algorithm is
scalable for different types of kernel functions for the same approximation
accuracy. The “3-sphere” dataset is used.

For each kernel function, we measure the time cost of the proposed algorithm
as the size of data $n$ increases. The three types of costs - hierarchical
data reduction, hierarchical matrix construction, matrix-vector multiplication
- correspond to the three plots in Figure 14.

In Figure 14, each plot shows the timing for all four kernels. The relative
error for each case is $10^{-6}$. It is easily seen from Figure 14 that each
cost scales linearly with data size $n$. The algorithm is able to maintain
accuracy for different types of kernel functions.

Table 1: Kernel functions used in the experiments in Section 6.1. Here $\kappa_{1}(x,x)=0$. $\kappa_{1}(x,y)$ | $\kappa_{2}(x,y)$ | $\kappa_{3}(x,y)$ | $\kappa_{4}(x,y)$
---|---|---|---
$\frac{1}{|x-y|}$ | $\exp(-|x-y|^{2})$ | $\cos(x\cdot y)$ | $\exp{\left(-\frac{1}{1-0.1|x-y|^{2}}\right)}$

####  6.1.3 HiDR once for all

In this experiment - called “HiDR once for all” - we perform HiDR on the
dataset to obtain representor sets and then use the representor sets to
construct hierarchical matrix representation for the different kernel
functions in Table 1 and Gaussian kernels with different bandwidths. The key
here is that, for a fixed compression level, HiDR is only performed _once_ and
the same representor sets are used for _all_ kernels.

Figure 15 shows the approximation error of the hierarchical matrix with
respect to the average size of farfield representor sets $Y_{i}^{*}$.
Different error curves correspond to approximations to different kernels. We
see that by increasing the size of the farfield representor sets, the matrix
approximation error is reduced effectively for all kinds of kernel functions.
Since HiDR is only applied once, the precomputation cost is almost negligible
when amortized over multiple kernels. The accuracy as seen from Figure 15
justifies the data-driven construction with the efficient kernel-independent
HiDR. We see that the data-driven approach is particularly useful when
hierarchical matrices for different kernel functions or kernel parameters need
to be computed.

In applications like Gaussian processes, the bandwidth parameter for the
Gaussian kernel is unknown and is determined by an iterative algorithm.
Therefore, the bandwidth changes constantly, thus the kernel function. This
makes existing hierarchical matrix constructions inefficient because the
entire hierarchical algorithm needs to be run from scratch every time the
bandwidth changes. The proposed data-driven approach, however, performs data
reduction only _once_ and no matter what the bandwidth is, the hierarchical
matrix representation can be constructed rapidly based on the computed
representor sets. It can be seen from Figure 15(right) that for a wide range
of bandwidth values, the approximation error decays effectively as more points
are used in the representor sets. The nearly zero approximation error for the
bandwidth $L=0.01$ is due to the fact that the admissible block is almost a
zero matrix as $\exp(-|x-y|^{2}/0.01^{2})\approx 0$ when $x$ and $y$ are away
from each other.

Overall, it can be seen from the experiments that the data-driven approach
serves as a black-box tool for rapidly computing hierarchical matrices for
general kernel functions. It is especially efficient in the situation when
multiple kernel functions need to be approximated.

Figure 14: Section 6.1.2 scaling test for kernels in Table 1: CPU time for
HiDR (left), hierarchical matrix construction (middle), matrix-vector
multiplication (right)

Figure 15: Section 6.1.3: perform HiDR only once to obtain representor sets
and then use them to construct hierarchical matrices for multiple kernels:
matrix approximation error vs average size of farfield representor sets
$Y_{i}^{*}$. Left: kernels in Table 1; Right: Gaussian kernel
$\exp(-|x-y|^{2}/L^{2})$ with different bandwidth $L=10^{k}$ with
$k=-2,-1,0,1,2$

###  6.2 Comparison to special-purpose methods for the Coulomb kernel

In this section, we compare the new general-purpose data-driven (‘DD’)
construction to several optimized packages for the Coulomb kernel
$\kappa(x,y)=\frac{1}{|x-y|}$, for example,
FMM3D222https://fmm3d.readthedocs.io/en/latest/[21], PVFMM [36], and Proxy
Surface method (cf. [37, 20, 27]) implemented in H2Pack [34]. These methods
are specialized for the Coulomb kernel to offer better efficiency in practice
than the interpolation-based methods for constructing $\mathcal{H}^{2}$
matrices.

DD, FMM3D, PVFMM and H2Pack are compiled using Intel C/C++/Fortran compiler
v19.0.5 with optimization flags “-xHost -O3”. Intel MKL 19.0.5 is used in all
tested libraries to perform general matrix-vector and matrix-matrix
multiplications. DD, H2Pack, and FMM3D use one thread per CPU core and 24
cores on one computing node. PVFMM uses MVAPICH2 2.3.2 as the MPI backend and
uses one MPI process with 24 cores on one computing node.

We use the same datasets as in Section 6.1. For every method, the total time
is computed as

$\text{total time = precomputation + $\mathcal{H}^{2}$ construction + matrix-
vector multiplication}.$

FMM3D and proxy surface method do not have precomputation, while PVFMM and DD
require precomputation. For DD, the precomputation refers to HiDR.

Timings for precomputation, hierarchical build, matvec, are shown in Figure
16, Figure 17, Figure 18, respectively. The total time is shown in Figure 19.
The relative error (in 2-norm) for each test is $10^{-6}$.

From Figure 16, we see that the proposed hierarchical data-reduction (HiDR)
requires significantly lower precomputation cost compared to PVFMM. This is
due to the fact that the kernel matrix is _never_ accessed in HiDR and _no_
algebraic compression is computed. Moreover, we see from Figure 16 that the
advantage of the data-driven method becomes more obvious for irregular data
from a manifold, such as 3-sphere and Dino. It can be seen that PVFMM has
almost constant cost independent of the size of the dataset $n$. This is
because PVFMM treats matrix compression as a continuous problem (thus
independent of the size of data) and the precomputation involves solving
integral equations on spheres to facilitate the farfield compression.

From Figure 17, we see that FMM3D outperforms other methods in hierarchical
construction. This is because, unlike other methods, FMM3D does _not_ compute
and store a hierarchical matrix representation. Instead, it computes the
hierarchical representation when performing matrix-vector multiplication. The
other methods have similar performance, where no single method performs
significantly better than others across all datasets. It should be noted that,
for methods with precomputation, the hierarchical construction time can be
further reduced at the expense of more precomputation time. In principle, more
time spent in precomputation could yield faster hierarchical construction.

For matvec, Figure 18 shows that the data-driven method and proxy surface
method achieve similar performance that is in general better than FMM3D and
PVFMM. FMM3D is significantly slower than other methods due to the on-the-fly
hierarchical construction. For data sampled from a manifold, PVFMM is
outperformed by DD and Proxy Surface. The results in Figure 18 justify the
efficiency of the hierarchical matrix representation built from the fast HiDR
in Figure 16.

For the total computation time, as can be seen from Figure 19, we see that
Proxy Surface provides the best performance overall, followed by DD. The
advantage of DD is more evident for data from a manifold, e.g. 3-sphere, Dino.
In general, we see that the data-driven method leads to a lot more
computational savings when data is sampled from a low-dimensional manifold.

It should be emphasized that the Proxy Surface method is a specialized method
optimized for the Coulomb kernel to offer superior efficiency, while DD is a
general-purpose approach that can be applied to a variety of kernel functions
(cf. Table 1). Unlike the special-purpose methods, no analytic property of the
kernel function is used in the data-driven hierarchical construction. It is
nonetheless possible to design specialized data-driven algorithms for the
kernel function of interest to improve efficiency. Overall, we conclude from
the results in Figure 19 that the data-driven method, as a black-box tool for
hierarchical matrix computations, also offers excellent efficiency for special
kernels without utilizing any specific property of the kernel.

Figure 16: Section 6.2 experiment: Precomputation time of PVFMM and DD for
approximating $n$-by-$n$ Coulomb kernel matrices with datasets Cube, 3-sphere,
Dino

Figure 17: Section 6.2: Hierarchical construction time (after precomputation)
of FMM3D, PVFMM, Proxy Surface and DD for approximating $n$-by-$n$ Coulomb
kernel matrices with datasets Cube, 3-sphere, Dino

Figure 18: Section 6.2: Matrix-vector multiplication time of FMM3D, PVFMM,
Proxy Surface and DD for approximating $n$-by-$n$ Coulomb kernel matrices with
datasets Cube, 3-sphere, Dino

Figure 19: Section 6.2: Total construction time of FMM3D, PVFMM, Proxy Surface
and DD for approximating $n$-by-$n$ Coulomb kernel matrices with datasets
Cube, 3-sphere, Dino

###  6.3 Memory efficiency

In this section, we illustrate the memory efficiency of the proposed data-
driven approach by comparing it to interpolation-based hierarchical matrix
construction. We test the two general-purpose methods for the different kernel
functions listed in Table 1. The “3-sphere” dataset is used with $n=20000$
points.

In Figure 20, we plot the approximation error vs the memory use for each
method and each kernel function. The memory use is measured by the cost for
storing the hierarchical representation derived by the respective method. The
high memory use of interpolation-based construction is clearly seen from the
plots. In three dimensions, the number of interpolation nodes is $k^{3}$ if
$k$ interpolation nodes are used in each dimension. This number may exceed the
size of the admissible block to be approximated. We found that, in practice,
to achieve moderate to high approximation accuracy, a large number of
interpolation nodes is needed. The proposed data-driven method, on the other
hand, significantly reduces the memory needed to achieve a certain
approximation accuracy. Equivalently, we can also conclude that, for the same
approximation rank and memory requirement, the data-driven method is able to
provide a much more accurate hierarchical representation than the one derived
from interpolation. For large scale data, the memory efficiency of data-driven
construction would be even more prominent.

Figure 20: Section 6.3 experiment: error vs memory use of interpolation-based
and data-driven constructions for approximating the four kernel matrices
(Table 1) with the 3-sphere dataset

##  7 Conclusion

We proposed general-purpose data-driven hierarchical matrix construction
accelerated by a novel hierarchical data reduction (HiDR). The algorithm first
computes a reduced data representation following the tree structure and then
performs the hierarchical low-rank compression. Different from all existing
methods, HiDR entirely operates on the given dataset, _without_ accessing the
kernel function or kernel matrix. The complexity of the whole data-driven
construction is linear with respect to the data size. Compared to general-
purpose methods such as interpolation, the new data-driven framework requires
less memory for the same matrix approximation accuracy. For special kernels
like the Coulomb kernel, the general data-driven method, as a black-box
approach, demonstrates competitive performance when compared to specialized
methods optimized for the Coulomb kernel. The data-driven approach yields low
computational cost and is particularly efficient for data sampled from low-
dimensional manifolds. Future work includes extending the data-driven
framework to the efficient construction of hierarchical matrices for high
dimensional data in machine learning applications. One appealing feature is
that, when constructing hierarchical matrices for different kernel functions
associated with the same dataset, HiDR only needs to be performed once, which
significantly reduces the total computational cost. Additionally, the data
driven procedure can be optimized towards the special kernel function under
consideration. In the current presentation, we focus on providing a general
approach, which could be useful for general-purpose library for accelerating
kernel matrix computations with hierarchical matrix representations. Possible
improvements for special kernel functions will be investigated at a future
date.

## References

  * [1] C. R. Anderson.  An implementation of the fast multipole method without multipoles.  SIAM J. Sci. Statist. Comput., 13(4):923–947, 1992.
  * [2] K. E. Atkinson.  The numerical solution of the eigenvalue problem for compact integral operators.  Transactions of the American Mathematical Society, 129(3):458–465, 1967.
  * [3] K.E. Atkinson.  Convergence rates for approximate eigenvalues of compact integral operators.  SIAM Journal on Numerical Analysis, 12(2):213–222, 1975.
  * [4] J. Barnes and P. Hut.  A hierarchical O(N log N) force-calculation algorithm.  Nature, 324:446–449, December 1986.
  * [5] M. Bebendorf.  Approximation of boundary element matrices.  Numer. Math., 86(4):565–589, 2000.
  * [6] M. Bebendorf and R. Venn.  Constructing nested bases approximations from the entries of non-local operators.  Numer. Math., 121(4):609–635, 2012.
  * [7] Christopher M. Bishop.  Pattern Recognition and Machine Learning.  Springer, 2006.
  * [8] S. Börm and L. Grasedyck.  Hybrid cross approximation of integral operators.  Numer. Math., 101(2):221–249, 2005.
  * [9] S. Börm, L. Grasedyck, and W. Hackbusch.  Introduction to hierarchical matrices with applications.  Eng. Anal. Bound. Elem., 27(5):405–422, 2003.
  * [10] D. Cai.  Physics-Informed Distribution Transformers Via Molecular Dynamics and Deep Neural Networks.  preprint, https://ssrn.com/abstract=4028725.
  * [11] D. Cai and Z. Cai.  Hybrid a posteriori error estimators for conforming finite element approximations to stationary convection-diffusion-reaction equations.  arXiv preprint arXiv:2107.06341.
  * [12] D. Cai and Z. Cai.  A hybrid a posteriori error estimator for conforming finite element approximations.  Computer Methods in Applied Mechanics and Engineering, 339:320 – 340, 2018.
  * [13] D. Cai, Z. Cai, and S. Zhang.  Robust equilibrated a posteriori error estimator for higher order finite element approximations to diffusion problems.  Numerische Mathematik, 144(1):1–21, 2020.
  * [14] D. Cai, Z. Cai, and S. Zhang.  Robust equilibrated error estimator for diffusion problems: mixed finite elements in two dimensions.  Journal of Scientific Computing, 83(1):1–22, 2020.
  * [15] D. Cai, E. Chow, L. Erlandson, Y. Saad, and Y. Xi.  SMASH: Structured matrix approximation by separation and hierarchy.  Numerical Linear Algebra with Applications, 25(6):e2204, 2018.
  * [16] D. Cai, E. Chow, and Y. Xi.  Data-driven Linear Complexity Low-rank Approximation of General Kernel Matrices: A Geometric Approach.  preprint, 2022.
  * [17] D. Cai, J. G. Nagy, and Y. Xi.  Fast Deterministic Approximation of Symmetric Indefinite Kernel Matrices with High Dimensional Datasets.  SIAM J. Matrix Anal. Appl., accepted.
  * [18] D. Cai and P. S. Vassilevski.  Eigenvalue problems for exponential-type kernels.  Computational Methods in Applied Mathematics, 20(1):61–78, 2020\.
  * [19] D. Cai and J. Xia.  A stable matrix version of the fast multipole method: stabilization strategies and examples.  Electron. Trans. Numer. Anal., 54:581–609, 2021.
  * [20] Hongwei Cheng, Zydrunas Gimbutas, Per-Gunnar Martinsson, and Vladimir Rokhlin.  On the compression of low rank matrices.  SIAM Journal on Scientific Computing, 26(4):1389–1404, 2005.
  * [21] Hongwei Cheng, Leslie Greengard, and Vladimir Rokhlin.  A fast adaptive multipole algorithm in three dimensions.  Journal of computational physics, 155(2):468–498, 1999.
  * [22] D. Colton and R. Kress.  Integral equation methods in scattering theory.  Pure and applied mathematics. Wiley, 1983.
  * [23] Dennis Decoste and Bernhard Schölkopf.  Training invariant support vector machines.  Machine learning, 46(1-3):161–190, 2002.
  * [24] Yuval Eldar, Michael Lindenbaum, Moshe Porat, and Yehoshua Y Zeevi.  The farthest point strategy for progressive image sampling.  IEEE Transactions on Image Processing, 6(9):1305–1315, 1997.
  * [25] L. Erlandson, D. Cai, Y. Xi, and E. Chow.  Accelerating parallel hierarchical matrix-vector products via data-driven sampling.  In 2020 IEEE International Parallel and Distributed Processing Symposium (IPDPS), pages 749–758, 2020.
  * [26] Bengt Fornberg and Grady Wright.  Stable computation of multiquadric interpolants for all values of the shape parameter.  Computers & Mathematics with Applications, 48(5-6):853–867, 2004\.
  * [27] Adrianna Gillman, Patrick M. Young, and Per-Gunnar Martinsson.  A direct solver with o(n) complexity for integral equations on one-dimensional domains.  Frontiers of Mathematics in China, 7(2):217–247, 2012.
  * [28] L. Greengard and V. Rokhlin.  A fast algorithm for particle simulations.  J. Comput. Phys., 73:325–348, 1987.
  * [29] M. Gu and S. C. Eisenstat.  Efficient algorithms for computing a strong rank-revealing QR factorization.  SIAM J. Sci. Comput., 17(4):848–869, 1996.
  * [30] W. Hackbusch.  Hierarchical Matrices: Algorithms and Analysis.  Springer Series in Computational Mathematics. Springer Berlin Heidelberg, 2015.
  * [31] W. Hackbusch and S. Börm.  $\mathcal{H}^{2}$-matrix approximation of integral operators by interpolation.  Appl. Numer. Math., 43(1):129–143, 2002.
  * [32] W. Hackbusch, B. N. Khoromskij, and R. Kriemann.  Hierarchical matrices based on a weak admissibility criterion.  Computing, 73(3):207–243, 2004.
  * [33] W. Hackbusch, B.N. Khoromskij, and S.A. Sauter.  On $\mathcal{H}^{2}$-matrices.  In Hans-Joachim Bungartz, Ronald H. W. Hoppe, and Christoph Zenger, editors, Lectures on Applied Mathematics, pages 9–29. Springer, Berlin, 2000.
  * [34] H. Huang, X. Xing, and E. Chow.  H2Pack: High-performance H-2 matrix package for kernel matrices using the proxy point method.  ACM Transactions on Mathematical Software (TOMS), 47(1):1–29, 2020\.
  * [35] H. B. Keller.  On the accuracy of finite difference approximations to the eigenvalues of differential and integral operators.  Numerische Mathematik, 7(5):412–419, 1965.
  * [36] D. Malhotra and G. Biros.  Pvfmm: A parallel kernel independent fmm for particle and volume potentials.  Communications in Computational Physics, 18(3):808–830, 2015.
  * [37] P.G. Martinsson and V. Rokhlin.  A fast direct solver for boundary integral equations in two dimensions.  J. Comput. Phys., 205(1):1–23, 2005.
  * [38] J.E. Osborn.  Spectral approximation for compact operators.  Mathematics of computation, 29(131):712–725, 1975.
  * [39] Gabriel Peyré and Laurent D Cohen.  Geodesic remeshing using front propagation.  International Journal of Computer Vision, 69(1):145–156, 2006.
  * [40] V. Rokhlin.  Rapid solution of integral equations of classical potential theory.  Journal of Computational Physics, 60(2):187–207, 1985.
  * [41] V. Rokhlin.  Rapid solution of integral equations of scattering theory in two dimensions.  Journal of Computational Physics, 86(2):414–439, 1990.
  * [42] Thomas Schlömer, Daniel Heck, and Oliver Deussen.  Farthest-point optimized point sets with maximized minimum distance.  In Proceedings of the ACM SIGGRAPH Symposium on High Performance Graphics, pages 135–142, 2011.
  * [43] A. Spence.  Error bounds and estimates for eigenvalues of integral equations.  Numerische Mathematik, 29(2):133–147, 1978.
  * [44] X. Sun and N.P. Pitsianis.  A matrix version of the fast multipole method.  SIAM Rev., 43(2):289–300, 2001.
  * [45] 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.