A Study of Clustering Techniques and Hierarchical Matrix Formats for Kernel Ridge Regression
Elizaveta Rebrova Gustavo Chavez Yang Liu Pieter Ghysels Xiaoye Sherry Li

	





A Study of Clustering Techniques and Hierarchical Matrix Formats

 for Kernel Ridge Regression
    Elizaveta Rebrova[Department of Mathematics,
University of Michigan, 530 Church St, Ann Arbor, MI 48109,
U.S.A. erebrova@umich.edu] Gustavo Chávez[Computational Research Division,
Lawrence Berkeley National Laboratory, 1 Cyclotron Road, Berkeley,
CA 94720, U.S.A. {gichavez,liuyangzhuan,pghysels,xsli}@lbl.gov]  Yang Liu[2] Pieter
Ghysels[2] Xiaoye Sherry Li[2]
    
=============================================================================================================================================================================================================================================================================================================================================================================





We present memory-efficient and scalable algorithms for kernel
methods used in machine learning. 
Using hierarchical matrix approximations for the kernel matrix the memory requirements,
the number of floating point operations, and the execution time are
drastically reduced compared to standard dense linear algebra
routines. We consider both the general ℋ matrix hierarchical
format as well as Hierarchically Semi-Separable (HSS) matrices. Furthermore, we
investigate the impact of several preprocessing and clustering
techniques on the hierarchical matrix compression. Effective
clustering of the input leads to a ten-fold increase in efficiency
of the compression. The algorithms are implemented using the	
STRUMPACK solver library. These results confirm that —
with correct tuning of the hyperparameters — classification using
kernel ridge regression with the compressed matrix does not lose
prediction accuracy compared to the exact — not compressed —
kernel matrix and that our approach can be extended to
𝒪(1M) datasets, for which computation with the full
kernel matrix becomes prohibitively expensive. We present numerical 
experiments in a distributed memory environment up to 1,024 processors
of the NERSC's Cori supercomputer using well-known datasets 
to the machine learning community that range from dimension 8 up to 784.







§ INTRODUCTION

	
Kernel methods play an important role in a variety of applications in
scientific computing and machine learning (see, for example,
<cit.>). The idea is to implicitly map a set of data to
a high-dimensional feature space via a kernel function, which allows
performing a more sensitive training procedure. Given n data points
x_1, …, x_n ∈^d and a kernel function
: ^d ×^d →, the corresponding kernel matrix
K ∈^n× n is defined as K_ij = (x_i, x_j). Solving
linear systems with kernel matrices is an algebraic procedure required
by many kernel methods. One of the simplest examples is kernel ridge
regression in which one solves a system with the matrix
K + λ I, with λ a regularization parameter and I the
identity matrix. The limitation of this approach is in the lack of
scalability. The number of data-points, n, is typically very large,
and a direct solve would require 𝒪(n^3) operations, even
requiring 𝒪(n^2) complexity just to construct the full
kernel matrix.

Acceleration of kernel methods has been studied extensively in
scientific computing research, e. g. <cit.>, mostly
using low-rank matrix approximations for K. This is, however, based
on an assumption which is not valid in general. Consider for example one of the 
kernel matrices, defined by	the Gaussian radial basis function:
    

    K_ij :=  exp( - 1/2x_i - x_j^2/h^2).

    
Note that for h → 0, K approaches the identity matrix, while for
h →∞ it is nearly a rank one matrix (all elements
K_ij∼ 1). Intermediate values of h interpolate between these
“easy” cases.  The value of h is determined by the expected accuracy of
the machine learning algorithm on a certain dataset (e.g. by 	cross-validation).
Therefore, we cannot simply assume the “easy” structure of the
matrix.  However, the off-diagonal part of the kernel matrix typically
has a fast singular value decay, which means kernel matrices are good
candidates for hierarchical low-rank solvers like
STRUMPACK <cit.>.
	
For a discussion on the existence of efficient far-field compression of the off-diagonal blocks we refer the reader to <cit.> and the references therein.
To empirically confirm the off-diagonal low-rank property, we examine one dataset,
GAS1K, whose kernel matrix K has dimension N=1000. In
Figure <ref>, we plot the singular values of the off-diagonal block K(1,2) of size 500, with h values varying from small to large. Figure <ref> plots the singular values of the entire kernel matrix. We use both the natural ordering of the rows/columns
as well as a reordering of the rows/columns based on a recursive
two-means (2MN) clustering algorithm applied to the input data (see
Section <ref>).  For the same K(1,2) block,
Table <ref> lists the number of singular values
larger than 0.01; we call this the effective rank.  As can be seen
from both Figure <ref> and Table <ref>,
The 2MN preprocessing leads to much faster decay for h ∼ 1,
significantly improves the potential benefits of a solvers that
exploits the off-diagonal low-rank property. Reordering of the input
data is the main subject of this paper.	
	



    


 §.§ Main Contributions


We propose to use hierarchical matrix approximations for the
kernel. We use both Hierarchically Semi-Separable (HSS) as well as
ℋ matrices. The HSS algorithms are implemented as part of
the STRUMPACK library. STRUMPACK — STRUctured Matrix PACKage — is
a fast linear solver and preconditioner for sparse and dense systems.
One salient feature of STRUMPACK is that the HSS construction uses
adaptive randomized sampling for numerical rank detection, which
requires a black-box matrix times vector multiplication routine as
well as access to selected elements from the kernel
matrix (<ref>). There is no need to explicitly store the whole
matrix K; we call this the partially matrix-free interface.
This feature is particularly attractive for kernel methods, because
forming the complete K may consume too much memory.

Although STRUMPACK works algebraically with any input matrix,
different row and column orderings affect the numerical ranks of the
off-diagonal blocks, and hence the performance.  Intuitively, for an
off-diagonal block K(I,J) to have a small rank, the data points in
the clusters I and J need to have as little interaction as
possible. This suggests a preprocessing step using clustering
algorithms, namely, (a) find groups of points with large inter-group
distances and small intra-group distances, and (b) permute the rows
and columns of the matrix K such that the points of each group have
consecutive indices.
	
Another performance-critical aspect is how the random sampling is
performed. If we use traditional matrix-matrix multiplication to
perform the sampling K R, where R consists of a number of random
vectors, the entire solution time will be dominated by this
operation. Here, we can exploit the special structure of matrix K
and use a faster structured sampling method with an ℋ
matrix.
	
To summarize, the main contributions of this work are:


	
  * We present an algorithm allowing fast approximate kernel matrix computations with linear scalability of the factorization and solution phases.
	
  * For the preprocessing step, we explore various row and column orderings to improve the efficiency of HSS approximations in kernel methods.
	
  * We use a fast sampling method based on the ℋ matrix approximation of K, removing the bottleneck of the sampling phase for the HSS construction.
	
  * We use the auto-tuning framework OpenTuner for the tuning of the hyperparameters in classification using kernel ridge regression.
	
  * We report scalability experiments up to 1,024 cores, with real-world datasets from the UCI <cit.> up to N=4.5M data points for training, and dimension 784.

	
We experimented with a number of clustering techniques (including
versions of agglomerative and hierarchical clusterings, as well as
divisive 2-means, kd-tree and PCA-tree clusterings) and achieved
improvements up to 10 × in terms of memory usage, with much
lower ranks, of the compressed matrix compared to the naive
application of STRUMPACK (without clustering) and up to 4 ×
compared to k-d tree based reordering.  We believe that the clustering
techniques might be useful for more general classes of data related
matrices, in order to reveal the implicit hierarchical off-diagonal
low-rank structure, and will make them amenable to fast and scalable
algorithms. Asymptotically quasi-optimal memory consumption is key for the kernel ridge regression to be able to process large datasets.
	


 §.§ Previous Work


A number of methods have been developed in the intersection of scientific
computing and machine learning, with the goal to accelerate kernel methods.
Here we give some highlights of the important methods, although we do not
have an exhaustive list covering the entire field.

The best rank r approximation of a matrix would be given by direct
low-rank matrix factorization methods, based on singular value decomposition.
However, even with multiple improvements (in particular, due to the use of
randomization <cit.>, see also a survey <cit.>),
the total complexity of such methods stays quadratic in the number of
samples 𝒪(n^2log r). They also require creating and storing the
entire kernel matrix K.

When the kernel matrix exhibits globally low rank, Nyström methods
are shown to be among the best methods (see, for example, <cit.>).
Unfortunately, not all kernel matrices can be well approximated by
low-rank matrices in a global sense.

Other approximation approaches include fast multipole method (FMM)
partitions that try to split the data points into different boxes and
quantify interactions between the points that are far apart
(such as <cit.>) – however, it seems that these methods work well
only for low dimensional data. 
	
The next idea, also crucial for our work, is to combine the low rank approximation with a clustering of the data points. Some of the previous algorithms using this idea are the following. Clustered Low-Rank Approximation (CLRA, <cit.> and its parallel version <cit.>) was created to process large graphs. It starts with the clustering of the adjacency matrix, and then computes a low-rank approximation of each cluster (i.e., diagonal block), e.g. using singular value decomposition. Memory Efficient Kernel Representation (MEKA, <cit.>)  also performs clustering, and then applies Nyström approximation to within-cluster blocks (to avoid computing all within block entries), and introduces a sampling approach to capture between-block information. Block Basis Factorization (BFF, <cit.>) improves upon MEKA in several ways, including computation of the low-rank basis vectors from a larger space, introducing a more sophisticated sampling procedure, and estimating near-optimal k, the number of clusters for the initial data splitting.

Our approach exploits both clustering and low rank property, but in a different way. Instead of initial splitting of data into k clusters, we construct a binary tree of embedded clusters. Then we use low rank property of the off-diagonal sub-blocks of the kernel matrix, corresponding to the inter-cluster links, as well as the hierarchical connection between the clusters (and, respectively, sub-blocks) for efficient approximation of the low rank bases.
    
The line of work that is closest to ours and also based on off-diagonal
low rank is presented in a series of
papers <cit.>,
where the block-diagonal-plus-low-rank hierarchical matrix format
is used to approximate the kernel matrix. The authors first developed an 
𝒪(d N log N) algorithm ASKIT to construct the approximate
representation for the kernel matrix <cit.>, and later
an 𝒪( N log N) algorithm INV-ASKIT to perform a
factorization of the approximate matrix <cit.>,
which can be used as a direct linear solver.
	
Our current work differs from the INV-ASKIT approach in several ways:

1) we use the ℋ and HSS matrix formats in the approximation, 2) we use ULV
factorization, rather than the Sherman-Morrison-Woodbury
formula used in INV-ASKIT, and 3) we compare a number of
clustering methods to help reduce the off-diagonal rank, while
INV-ASKIT only used the k-d tree ordering.
	



§ KERNEL RIDGE REGRESSION FOR CLASSIFICATION

    
Ridge regression is probably the most elementary algorithm that can be kernelized. Classical ridge regression is designed to find the linear hyperplane that approximates the data labels well, and at the same time does not have too large coefficients, namely

    _w∑_i = 1^n (y_i - w· x_i)^2 + λw_2^2,

where x_i are data points (rows of the data matrix X), y_i's are their labels, y = (y_1, …, y_n) and w is the normal vector to the target hyperplane. It can be proved (see, for example, <cit.>) that the optimal w is given by

    ŵ := X^T (X X^T + λ I)^-1 y.

The kernel trick introduces a way to replace the matrix X X^T by a certain kernel matrix K, effectively substituting the scalar products x_i· x_j = X_ij by the elements K(x_i, x_j) = K_ij, that represent the scalar product in some higher dimensional space. 
    
Specifically, our experimental results for the two-class classification using ridge regression with the Gaussian kernel (<ref>) are obtained by the following Algorithm <ref>:



The choice of parameters (h, λ) is based on a particular dataset and usually made by a cross-validation. The prediction accuracy is
computed as the fraction of correctly predicted test labels:
	
    Acc = |{y' ≡y̅'}|/m,

where y̅' ∈{± 1}^m are the true test labels (and y' ∈{± 1}^m are predicted test labels).

In Algorithm <ref>, the most time-consuming
computation is solving the linear system during the training stage
(Step 2).  One key observation is that the final prediction accuracy
does not require many digits of the solution weight vector w, since
it contributes only to the sign calculation in Step 4. Therefore, we
can use an inexact but faster linear solver, such as STRUMPACK, to
alleviate the performance bottleneck at Step 2. The goal of preprocessing Step 0 is to have  matrix K (constructed on Step 1) in an HSS format, at minimum rank.
    
Other algorithms are also suggested for the kernel matrix compression
(see <cit.>). In <cit.> the authors show the
results of the binary classification by the kernel ridge regression.
There does not seem to be one best solution suggested so far.

A possible future work is to fully compare the classification accuracy
and the speed with STRUMPACK to that obtained from the method
described in <cit.>.
    
Finally, although Algorithm <ref> is defined for two classes, it is easy to adapt it for the multi-class classification. The simplest way to do it is to make one-vs-all predictions. To distinguish between c > 2 classes, we would need to construct c binary classifiers, that differ from the Algorithm <ref> only in Step 4. Namely,

    y'(c)_i := |w^T(c) · K'(i)|,
 
which interprets now as the level of confidence that i-th test point belongs to the class c. Then the class of x'_i is defined as

    c(i) := _c y'(c)_i.





§ DATA SPARSE FORMATS

In this section, we briefly describe two hierarchical matrix formats
used in this work, the hierarchically semi-separable (HSS)
representation, in Section <ref>, and the hierarchical
ℋ matrix format, in Section <ref>. These matrix
representations use a hierarchical partitioning of the matrix into
smaller blocks, some of which can be compressed using low-rank
approximations.

The two hierarchical matrix formats play different roles in this work.
The HSS matrix has higher construction cost (if we use traditional
dense matrix multiplication in the sampling stage) but very low
factorization and solve cost, whereas the ℋ matrix has
lower construction and matrix-vector multiplication cost but higher
factorization cost. Our approach is to use the fast matrix-vector
multiplication capability of the ℋ matrix to speed up the
HSS matrix construction.
	


 §.§ Hierarchically Semi-Separable Matrix Representation



    
The HSS representation, as illustrated in Figure <ref>, uses a
block 2 × 2 partitioning of a matrix
A ∈ℝ^n × n, with a similar partitioning applied
recursively to the diagonal blocks. This recursive partitioning
defines a tree, as illustrated in Figure <ref>. With a
node i in the tree, an index set
I_i ⊂{1,…,n} is associated. At the last
level of the recursion, the diagonal blocks, i.e., A(I_i, I_i) are
stored as (small) dense matrices. All off-diagonal blocks A_ij are
compressed using a low-rank factorization U_i B_ij
V_j^T. Moreover, the column basis matrix U_i, for a node i with
children c_1 and c_2 in the hierarchy is defined as
U_i = [ U_c_1     0;     0 U_c_2 ]Ũ_i, and hence only the smaller matrix Ũ_i is
stored at node i. Only at the leaf nodes, where
U_i ≡Ũ_i, are the U_i stored explicitly. A similar
relation holds for the V_i basis matrices, and is referred to as the
nested basis property. The HSS data structure is implemented in the
STRUMPACK (STRUctured Matrix PACKage)
library <cit.>. STRUMPACK is a sparse direct
solver and preconditioner. It uses HSS compression for the sparse
triangular factors. However, the HSS kernels implemented in STRUMPACK
can also be used directly on dense matrices, for instance coming from
integral equations, the boundary element method, electromagnetic
scattering etc. For the construction of HSS matrices, STRUMPACK
implements a randomized algorithm from <cit.>. This algorithm
requires a matrix times (multiple) vector product for the random
sampling phase. A fast multiplication routine is crucial to get good
performance. STRUMPACK also implements a ULV
factorization <cit.> algorithm, and a
corresponding routine to solve a linear system with the factored HSS
matrix. If a fast sampling routine is available, both HSS compression
and factorization have 𝒪(rn) complexity, with r the
maximum HSS rank.
	


 §.§ ℋ Matrix Representation 


        

ℋ (and ℋ^2) matrices are yet another group of
hierarchical matrix formats for building fast linear solvers. Contrary to
HSS, where all off-diagonal blocks are low-rank compressed
(weak admissibility), ℋ formats only compress
the well-separated sub-blocks (strong admissibility). For example, the off-diagonal blocks in Figure <ref> are recursively partitioned into smaller ones and low-rank compressed if they are well-separated (i.e., admissible). Strong admissibility
leads to well-bounded numerical ranks for admissible blocks in linear
systems arising from high-dimensional applications (e.g., 3D acoustic 
and electromagnetic scattering problems and high-dimensional kernel matrices).
Therefore, the construction of an ℋ matrix for the 
kernel matrix can be performed in quasi-linear time and memory.
However, it is well-known that the inversion of ℋ matrices
requires significantly higher computation overhead when compared to
weak-admissibility solvers, such as those based on HSS.
Experiments showed that using an ℋ solver to solve the
linear system in Step 2 of Algorithm <ref> is much
slower than HSS due to the inversion bottleneck.
Therefore, instead of using the ℋ solver to solve the system
directly, we use it only to compress the kernel matrix and then 
to accelerate the HSS construction as described below.

Recall that during the HSS construction, STRUMPACK requires rapid multiplication of the kernel matrix and its transpose to random vectors in the sampling stage.
Instead of direct application of the matrices to vectors, one can leverage Fast Gauss Transform <cit.>, analytical/algebraic fast multipole methods to accelerate the construction.
To this end, we tailored and adapted an ℋ solver <cit.> to
kernel matrices. The low-rank representation of an admissible block is computed via a hybrid-ACA scheme that constructs a low-rank factorization of its submatrix that represents interactions between closely located points. The selection of the submatrix, however, is based on a trade-off between efficiency and
accuracy of the factorization. 	
	




§ DATASET CLUSTERING


This section discusses how appropriate preprocessing of the
input data can drastically improve the effectiveness of HSS
approximation. Section <ref> describes a
number of clustering algorithms to reduce the ranks of the HSS
off-diagonal blocks.



 §.§ Kernel Matrices and HSS Structure


Intuitively, a kernel matrix can be viewed as a similarity matrix
K = (K_ij)_i,j = 1, …, n, where

    K_ij =  similarity score between  x_i  and  x_j.

Such a matrix K, defined by a set of objects x_1, …, x_n is
always square and symmetric. These objects can be anything from two
integers, two real valued vectors, to particles, words etc., provided
that we know how to compare them. Kernel matrices are used to improve
algorithms to classify these objects, 
e.g. allowing more sophisticated boundaries between classes. The Gaussian
kernel (<ref>), studied throughout this paper, is but one example
kernel matrix; though probably the most widely used one. Let
x_1, …, x_n be data points in d dimensional space ^d,
which are said to be “similar” if they are close to each other in
Euclidean distance.

Broadly speaking, the preprocessing takes advantage of the fact that
the interaction between two well separated clusters of data points can
be approximated accurately when expressed in terms of the interaction
between a smaller number of representative points from each
cluster. This same idea is used very successfully in tree-codes <cit.>, the
fast multipole method (FMM) <cit.>, matrix skeletonization <cit.> and
interpolative decomposition <cit.>. Hence, splitting the input data in
clusters with large inter cluster distances leads to lower ranks for
the off-diagonal blocks of K, and hence less memory usage and faster
algorithms. This low rank property is illustrated in
Figure <ref> with and without preprocessing, and it
motivates the idea of using an HSS solver like STRUMPACK to speed up
kernel matrix computations.  Due to the exponential decay of the
Gaussian kernel, many elements in the kernel matrix, away from the
diagonal, are actually negligibly small. Furthermore, the low rank pattern is
proven in <cit.>.



 §.§ Preprocessing by Data Clustering


Reordering the input data x_1, …, x_n corresponds to
applying a permutation symmetrically to the rows and columns
of the kernel matrix. Preprocessing is applied in the
following steps:



  * Partition the data points in two clusters with large inter-group distance and small intra-group distance.
	
  * Reorder the input data such that points in the same group occupy consecutive indices, increasing the data-sparsity of the off-diagonal blocks in the corresponding kernel matrix.
	
  * Use these two index ranges to partition the corresponding node in the HSS tree, as illustrated in Figure <ref>. Repeat Steps 1 and 2 recursively for both partitions computed in Step 1.


Step 1 suggests that this preprocessing is a clustering
problem. The task of finding clusters of points is widely
studied, with numerous clustering algorithms described in the
literature, without there being one absolute best solution.
Different algorithms perform better for different applications
and clustering quality can often be traded for execution time
or memory usage (see, for example, <cit.>).

Our specific requirements for an optimal clustering come from
(a) the properties of the HSS data structure and (b) the general
goal to minimize memory usage of the hierarchical matrix data
structures. The latter requires a clustering algorithm that
does not construct the kernel matrix K, or an equivalent
n × n distance matrix, explicitly. It should also be
fast, preferably 𝒪(n) or 𝒪(n log n).
The HSS partitioning is defined by a hierarchical structure; a
possibly unbalanced and incomplete tree. However, in order to
be able to exploit enough parallelism and to reduce memory
usage, the tree should be deep enough, i.e., have a small
enough maximum leaf size.

Rather than the standard dissimilarity metrics measuring
clustering quality, the following performance metrics are used
in this work:


    
  * Memory (MB): the sum of the memory used by all the individual smaller matrices in the HSS structure: D_i, U_i, V_i, B_ij, B_ji (see
    Section <ref>).
    
  * Accuracy of classification (%): the percentage of correctly predicted labels in the test set (with the parameters h and λ chosen based on the validation set).
    
  * Time (s): the time required for compression into HSS form, for factorization of the HSS matrix, and for solution of the linear system.
    
  * Maximum rank: the largest rank encountered in any of the off-diagonal blocks of the HSS structure.




 §.§ Selected Preprocessing
Methods


Previous work on approximation of kernel matrices used
reorderings based on ball tree clustering, see for
instance <cit.>, or based on k-d tree
clustering, see <cit.>.

We have compared a variety of different clustering techniques
and their variations, both divisive and agglomerative, on
several real world datasets. Each of the clustering algorithms
used in the numerical results section are divisive, i.e., they
use a top-down, recursive split of the input points into two
separated clusters. The recursive splitting continues while
clusters are bigger than a certain leaf size, chosen to be
16 for HSS. This leaf size is the size of the diagonal
blocks in the HSS structure. The leaf size should not affect
the accuracy of the hierarchical matrix representation, but it
affects the memory usage.

Agglomerative methods, in contrast to the divisive strategy, although very good at reducing memory and ranks of the HSS structure, did not show competitive performance.
We experimented with a variety of hierarchical clustering methods, and typical disadvantages arising were either very unbalanced class sizes, or lack of parallelism (𝒪(n^2) scaling, requiring to construct and store the complete distance matrix).
        
For the experiments we consider four orderings:

No preprocessing (NP): This is the baseline to
compare with: the input is not reordered, no information about
mutual distances is used to permute the matrix. The HSS tree
is a complete binary tree, constructed by recursively
splitting index sets in two equal (± 1) parts.

Recursive two-means (2MN) This special case of the
well-known k-means clustering algorithm is applied
recursively to define the HSS tree. K-means is an iterative
algorithm which works as follows. Pick two points (at random)
to represent two clusters; for each point in the data set,
compute the distance to those two points and find which one is
closer, assign points to the closest cluster; compute the
center off each cluster and take that as the new representative
point of the cluster; repeat until no points change cluster.
Typically only a few iterations are required. However, the
procedure is relatively sensitive to the choice of initial
cluster representatives. Initially, we pick one point randomly 
and select the second one with a probability proportional to 
the distance from the first one.

K-d tree (KD) The data is split along the coordinate
dimension of maximum spread, at the mean value for that
coordinate. Splitting at the mean is sensitive to outliers,
and can lead to very unbalanced trees. Alternatively,
splitting at the median value always results in a well
balanced tree. However, using the mean leads to lower memory
usage and when the input data is normalized, the sensitivity
with respect to outliers is less pronounced. If the resulting
clusters are still too unbalanced, i.e., when
(100·size(cluster1) < size(cluster2)), we
fall back to splitting at the median. This clustering is
applied recursively, where at each step of the recursion, a
new direction of maximum spread is determined.

Principal component analysis (PCA) At each step of
the recursive clustering, the data is split according to the
mean value in the projection onto the first principal
component (i.e. direction of the maximum spread). We expect
this to be a better clustering than the simpler k-d tree
method, at a somewhat higher cost.




§ NUMERICAL RESULTS
  

Experiments were performed at NERSC's Cori supercomputer. Each Cori node has two sockets, each socket is a 16-core Intel Xeon Processor E5-2698 v3 (“Haswell”) processor at 2.3 GHz and 128 GB DDR4 memory.



 §.§ Datasets Description

	
We use real-world datasets coming from life sciences, physical sciences and artificial intelligence. The reported datasets are: SUSY, HEPMASS (high-energy physics, Monte Carlo simulated kinematic properties of the particles in the accelerator), COVTYPE (predicting forest cover type from cartographic variables), GAS (measurements from chemical sensors to distinguish between gases with  different concentration levels), PEN and LETTER (handwritten digits and letters recognition). All mentioned datasets are taken from the UCI repository <cit.>. Finally, to illustrate the performance on a dataset with large dimension, we used the MNIST dataset of handwritten digits (extended 8M dataset, including shifts and rotations of the classical MNIST digits dataset) <cit.>.
	
For the datasets having multiple classes we perform one-vs-all prediction. In particular, in MNIST and PEN we predict digit 5, in LETTER we predict letter A, in COVTYPE – type 3 (Ponderosa Pine), in GAS – gas number 5. Prediction accuracy might differ significantly if one would predict some other class.
	


 §.§ Preprocessing Comparison


        

We report in Table <ref> our main performance metrics – memory and accuracy. The table compares different preprocessing methods with seven datasets. Memory usage heavily depends on parameter h, which is illustrated on the GAS10K dataset in Fig. <ref>. Furthermore, in the context of hierarchical low-rank approximations, memory is proportional to performance, since the number of flops is proportional to the numerical ranks of the approximation.
	


Recursive two means (2MN) preprocessing shows best memory performance 
for all the h values.
With STRUMPACK tolerance set to be at most 0.1, the prediction accuracy
does not seem to depend on the preprocessing methods,
see Table <ref>.
Moreover, for the 10K datasets reported in Table <ref> this accuracy matches the accuracy we get using the full non-compressed kernel matrix in Algorithm <ref>. The main disadvantage of 2MN is higher variance of the resulting rank, and, to lower extent, memory. The numbers reported for 2MN are average over three runs of the algorithm. The instability can be avoided by choosing non-random start points at every step of two means clustering. All datasets were normalized to have zero mean and unit standard deviation columns.	The experiments with non-normalized datasets, and with datasets normalized to have maximum absolute value one have shown significantly lower accuracy (e.g. for MNIST2M dataset).
	




 §.§ Hyperparameter Tuning

    
The algorithm parameters h and λ are key to determine the
predictive capabilities of the matrix approximation. Consider the HSS
approximation of matrix K + λ I used in this work.  When
the parameter λ changes, we only need to update the diagonal entries
of the HSS matrix, and there is no need to perform HSS construction again.

However, a change to h requires to perform HSS reconstruction from scratch,
which is costly.
There are theoretical estimates that limit the search space for
h <cit.>, but it is problem dependent.
A fine grid search is too costly, see Figure <ref>.
In contrast, the black-box optimization techniques in the OpenTuner
package <cit.> in Figure <ref> required only
100 runs and converged to a tuning parameter with better prediction
accuracies than grid search. This technique drastically reduced the
computational requirements to select h and λ.
	

	


 §.§ Large-Scale Prediction


The appeal of using optimal algorithms is the ability to process large amount of data. Table <ref> shows the predictive capabilities of this method by using datasets that allow the training step to use data points in the order of millions of entries, at different dimensions. The reported prediction accuracy is tested against a test subset of the complete dataset, that is, labeled data that were not considered during training or hyperparameter tuning.
	






 §.§ Asymptotic Complexity

    
When hierarchical matrix approximations have constant ranks, such as in a broad class of elliptic partial differential equations, the HSS memory consumption are number of operations are strict 𝒪(N); however, recent theoretical results show that the numerical rank of kernel matrices in high dimension depends on the dimension of the dataset <cit.>. We experimentally confirm this rank growth in the additional memory requirements (Figure <ref>) and time for factorization (Figure <ref>).
    
The major benefit of the near-linear complexity in factorization and memory is that this method enables the use of kernel matrices for the large datasets. As an example, storing a 1M dense matrix requires 8,000GB, whereas the HSS construction used in this work just required 1.3 GB. A similar argument can be made for the factorization of such a matrix, with which a traditional Cholesky factorization of 𝒪(n^3) is intractable at this scale.
    

	


 §.§ Performance Details

    
Table <ref> shows the performance of the main algorithmic steps of our method. The first step is the ℋ construction to accelerate the otherwise 𝒪(n^2) HSS sampling. For instance, with naïve HSS sampling, the construction process of a kernel matrix required more than 2 hours for n=0.5M, now with the ℋ matrix sampling we can construct an HSS matrix for n=4.5M in about 10 minutes.
    
Currently, our ℋ matrix implementation is only a prototype 
code and is not optimized at all. Although it enabled us to experiment
with datasets as large as reported in the literature, it is only capable of
effectively using a subset of the processes that the HSS code can
use (factorization and solve). Nonetheless, the successful synergy between
the ℋ and HSS matrix formats motivates our future work to develop 
a robust and more scalable distributed memory ℋ matrix code. With that, we expect to achieve much faster ℋ construction time and sampling time.


	
Figure <ref> shows a strong scaling experiment up to 1,024 cores of the factorization phase of the kernel computation for the datasets used in Table <ref>. As mentioned in the previous section, a scalable and fast factorization for large datasets is critical for kernel matrix computations. We refer the reader to <cit.> for the distributed memory parallelization aspects of the ULV factorization used in this work. At large core count, the number of degrees of freedom per core decreases dramatically, while communication time starts to dominates, hence runtime starts to depart from the linear scaling line. The wall-clock time, however, is in the order of a handful of seconds even for the largest dataset considered in this work. Note that even though the SUSY dataset is larger than the MNIST dataset the overall factorization time is larger for MNIST, this is due to the fact that the dimension of MNIST is bigger than SUSY, and the dimension of the dataset has a direct effect on the necessary ranks of the approximation, and therefore in the required number of operations.







§ CONCLUSIONS


We showed that the HSS linear solvers, such as the one as implemented in STRUMPACK, are useful for a completely new area of machine learning applications. We proposed to use several relatively sophisticated ways to preprocess the kernel matrix (i.e. cluster points in the dataset), and showed that the preprocessing can significantly improve the compression rate. In all the literature we have seen the authors used either natural or k-d tree preprocessing and did not perform comparison of the various ordering techniques.
    
We presented performance data of an HSS-based complexity solver for kernel matrices scaling up to 1,024 cores of the NERSC Cori supercomputer. The construction of the HSS matrices used different preprocessing methods to minimize memory consumption for the solution of the classification task on high dimensional datasets at high prediction  accuracy.

Preliminary results show that an incomplete factorization as described in this work might be an effective preconditioner for the iterative solution of kernel matrices. We will report on the trade-offs and effectiveness of this strategy in future work.
	


§ ACKNOWLEDGEMENTS

We thank Wissam Sid Lakhdar for the fruitful discussions about hyperparameter
tuning and for facilitating his OpenTuner scripts.

This research was supported by the Exascale Computing Project (17-SC-20-SC),
a collaborative effort of the U.S. Department of Energy Office of Science
and the National Nuclear Security Administration.

The research of the first author was supported in part by an appointment
with the NSF Mathematical Sciences Summer Internship Program sponsored
by the National Science Foundation, Division of Mathematical
Sciences (DMS).



99
		
rand kernelEl Alaoui, A. and Mahoney, M.W., 2014. Fast randomized kernel methods with statistical guarantees. stat, 1050, p.2.

BachBach, F., 2013, June. Sharp analysis of low-rank kernel matrix approximations. In Conference on Learning Theory (pp. 185-209).

libsvmChang, C.C. and Lin, C.J., 2011. LIBSVM: a library for support vector machines. ACM transactions on intelligent systems and technology (TIST), 2(3), p.27.

Biros-3Yu, Chenhan D., March, W.B., Xiao, B. and Biros, G., 2016, May. INV-ASKIT: a parallel fast direct solver for kernel matrices. In Parallel and Distributed Processing Symposium, 2016 IEEE International (pp. 161-171). IEEE.

Biros-5Yu, Chenhan D., March, W.B., Xiao, B. and Biros, G., 2017, May. An Nlog N parallel fast direct solver for kernel matrices. In Parallel and Distributed Processing Symposium, 2017 IEEE International.

GMGittens, A. and Mahoney, M.W., 2016. Revisiting the Nyström method for improved large-scale machine learning. The Journal of Machine Learning Research, 17(1), pp.3977-4041.

hastieHastie, T., Tibshirani, R. and Friedman, J., 2009. Unsupervised learning. In The elements of statistical learning (pp. 485-585). Springer, New York, NY.
kernel methodsHofmann, T., Schölkopf, B. and Smola, A.J., 2008. Kernel methods in machine learning. The annals of statistics, pp.1171–1220.

uciLichman, M. 2013. UCI Machine Learning Repository http://archive.ics.uci.edu/ml. Irvine, CA: University of California, School of Information and Computer Science.

Biros-4March, W.B., Xiao, B., Tharakan, S., Chenhan, D.Y. and Biros, G., 2015, November. A kernel-independent FMM in general dimensions. In High Performance Computing, Networking, Storage and Analysis, 2015 SC-International Conference for (pp. 1-12). IEEE.

Biros-2March, W.B., Xiao, B. and Biros, G., 2015. ASKIT: Approximate skeletonization kernel-independent treecode in high dimensions. SIAM Journal on Scientific Computing, 37(2), pp.A1089-A1110.

2Martinsson, P.G., 2011. A fast randomized algorithm for computing a hierarchically semiseparable representation of a matrix. SIAM Journal on Matrix Analysis and Applications, 32(4), pp.1251–1274.

kdtreeOmohundro, S.M., 1989. Five balltree construction algorithms (pp. 1-22). Berkeley: International Computer Science Institute.

1Rouet, F.H., Li, X.S., Ghysels, P. and Napov, A., 2016. A distributed-memory package for dense hierarchically semi-separable matrix computations using randomization. ACM Transactions on Mathematical Software (TOMS), 42(4), p.27.

SarlosSarlos, T., 2006, October. Improved approximation algorithms for large matrices via random projections. In Foundations of Computer Science, 2006. FOCS'06. 47th Annual IEEE Symposium on (pp. 143-152). IEEE.

SDSavas, B. and Dhillon, I.S., 2011, April. Clustered low rank approximation of graphs in information science applications. In Proceedings of the 2011 SIAM International Conference on Data Mining (pp. 164-175). Society for Industrial and Applied Mathematics.

SiSi, S., Hsieh, C.J. and Dhillon, I., 2014, January. Memory efficient kernel approximation. In International Conference on Machine Learning (pp. 701–709).

3STRUMPACK website http://portal.nersc.gov/project/sparse/strumpack/

MahoneyWang, R., Li, Y., Mahoney, M.W. and Darve, E., 2015. Structured block basis factorization for scalable kernel matrix evaluation. preprint arXiv:1505.00398.

Mah_surveyMahoney, M.W., 2011. Randomized algorithms for matrices and data. Foundations and Trends® in Machine Learning, 3(2), pp.123-224.

SuiSui, X., Lee, T.H., Whang, J.J., Savas, B., Jain, S., Pingali, K. and Dhillon, I., 2012, September. Parallel clustered low-rank approximation of graphs and its application to link prediction. In International Workshop on Languages and Compilers for Parallel Computing (pp. 76-95). Springer, Berlin, Heidelberg.

WSWilliams, C.K. and Seeger, M., 2001. Using the Nyström method to speed up kernel machines. In Advances in neural information processing systems (pp. 682-688).

BirosYu, C.D., Huang, J., Austin, W., Xiao, B. and Biros, G., 2015, November. Performance optimization for the k-nearest neighbors kernel on x86 architectures. In Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis (p. 7). ACM.
Vancouver

FMM1Yang, C., Duraiswami, R., Gumerov, N.A. and Davis, L., 2003, October. Improved fast Gauss transform and efficient kernel density estimation. In Ninth International Conference on Computer Vision (p. 664). IEEE.

HH. Guo and Y. Liu and J. Hu and E. Michielssen, 2017. A butterfly-based direct integral-equation solver using hierarchical LU factorization for analyzing scattering from electrically large conducting objects. IEEE Transactions on Antennas and Propagation, 65(9), p.4742.

TreeJ. Barnes and P. Hut, 1986. A hierarchical O(NlogN) force-calculation algorithm A fast algorithm for particle simulations. Nature, 324, p.446.

FMML. Greegard and V. Rokhlin, 1987. A fast algorithm for particle simulations. Journal of Computational Physics, 73, p.325.

IDCheng, Hongwei, Zydrunas Gimbutas, Per-Gunnar Martinsson, and Vladimir Rokhlin, 2005. On the compression of low rank matrices. SIAM Journal on Scientific Computing, 26(4), p.1389.

SkelS.A.Goreinov, E.E.Tyrtyshnikov and N.L.Zamarashkin, 1997. A theory of pseudoskeleton approximations. Linear Algebra and its Applications, 261(1), p.1.

RBFWang, R., Li, Y., and Darve, E., 2017. On the numerical rank of radial basis function kernel matrices in high dimension. arXiv:1706.07883.

OPENTUNERJ. Ansel and S. Kamil and K. Veeramachaneni and J. Ragan-Kelley and J. Bosboom and U. M. O'Reilly and S. Amarasinghe, 2014. OpenTuner: An extensible framework for program autotuning. In 23rd International Conference on Parallel Architecture and Compilation Techniques (PACT) (p. 303-315). IEEE.

ghysels2017robustGhysels, P., Li, X. S., Gorman, C. and
Rouet, F. H., 2017.  A robust parallel preconditioner for
indefinite systems using hierarchical matrices and randomized
sampling, IPDPS 2017, pp.897–906. IEEE

chandrasekaran2006fastChandrasekaran, S., Gu, M. and Pals,
T, 2006.  A fast ULV decomposition solver for Hierarchically
Semiseparable representations, SIAM SIMAX, 28(3), pp.603–622.

hESTIMATEW. Silverman, 1986. Density Estimation for Statistics and Data Analysis, Chapman and Hall.

march2015William B. March, George Biros, 2017. Far-field compression for fast kernel summation methods in high dimensions. Applied and Computational Harmonic Analysis. volume 43, issue 1, pages 39-75.
	


    
