Linear Time Kernel Matrix Approximation via Hyperspherical Harmonics

We present a new technique which leverages both the analytic properties of the underlying kernel (without being targeted to any specific kernel), and additional data-dependent compression to generate a low-rank approximation to the kernel matrix. The time and space cost of the approximation are linear in the size of the dataset. The main baseline technique with which we compare our results is the Nyström method for high-accuracy kernel-independent approximations. Our experiments show favorable comparison with the Nyström method for a host of kernels, dimensions, datasets, and error tolerances. In fact, we show a number of settings in which our method achieves a nearly-optimal tradeoff between rank and error, while still taking less time than the Nyström method.

Here we derive our analytic expansion and the subsequent additional compression step which comprise our new algorithm, which we refer to as a Harmonic Decomposition Factorization (HDF). We will begin by showing why an analytic expansion can lead to a low-rank approximation of a kernel matrix. Then we will show how a Chebyshev expansion of the kernel function can be expanded and rearranged into an expansion in hyperspherical harmonics. From there, our additional compression step examines the functions which form the coefficients in the harmonic expansion—when evaluated on the data points, these functions can be further reduced in number, thus reducing the overall rank of the factorization. Finally we present error and complexity analyses for the factorization.

We achieve a linear time low-rank approximation by analyzing the kernel’s functional form rather than algebraically operating on its associated matrix. Let $K$ be the kernel matrix whose entry in the $i$th row and $j$th column is $k(\|x_{i}-y_{j}\|)$ where $x_{i},y_{j}$ are $d$-dimensional points in datasets $\mathcal{X},\mathcal{Y}$ respectively, and $k$ is a piecewise smooth and continuous isotropic kernel. Our goal is to find $N\times r$ matrices $U,V$ with $r\ll N$ so that the $N\times N$ kernel matrix may be approximated by

We achieve this by finding functions $u_{l},v_{l}$ so that

Then we may define $U$ and $V$ as

This is the typical mechanism of analytic methods for kernel matrix approximation—$U$ and $V$ are formed by analyzing the kernel function and applying related functions to the data, rather than by examining the original kernel matrix. Thus our main goal is to find $u_{l},v_{l}$ so that (3) can be computed efficiently and is sufficiently accurate for as small $r$ as possible.

First we assume without loss of generality that ${0\leq\|x_{i}-y_{j}\|\leq 1}$ for all $x_{i},y_{j}$.²²2Mathematically, we’re absorbing the diameter of the point set into a lengthscale parameter in the kernel function. We begin by defining

and then forming a truncated Chebyshev expansion of $\bar{k}(r)$ for $-1\leq r\leq 1$

where $T_{i}$ is the $i$th Chebyshev polynomial of the first kind. The coefficients $a_{i}$ may be computed efficiently and accurately using e.g. the discrete cosine transform, see (Clenshaw & Curtis, 1960). This yields

for $0\leq\|x-y\|\leq 1$. Furthermore $a_{i}=0$ for odd $i$ since $\bar{k}$ is even by definition. To see this, note that ${a_{i}\coloneqq\int_{-1}^{1}\bar{k}(r)T_{i}(r)\mathrm{d}r}$, and when $i$ is odd the integrand is an even function times an odd function, hence the integral vanishes. Consequently, we henceforth assume that $p$ is even.

The coefficients in the Chebyshev polynomials are defined as $t_{i,j}$ so that

and, by definition of the Chebyshev polynomials,

Noting $\|x-y\|^{2}=\|x\|^{2}+\|y\|^{2}-2x^{T}y$ and plugging (6) into (5) yields

Applying the multinomial theorem and rearranging the sums (see Section A for details) results in

where we have absorbed $a_{i},t_{i,j}$, and the multinomial coefficients into constants $\mathcal{T}_{k_{1},k_{2},k_{3}}$.

Since $x^{T}y=\|x\|\|y\|\cos{\gamma}$ where $\gamma$ is the angle between the vectors $x,y$, we may write

Plugging this into (3.3) and performing additional rearranging/reindexing yields

where $\alpha=\frac{d}{2}-1$, $C_{j}^{\alpha}(\cos{\gamma})$ is the Gegenbauer polynomial of the $j$th order, and $\mathcal{T}^{{}^{\prime}}_{j,m,n}$ are constants which depend on the kernel but not on $x,y$.

We may finally put this in the form of (3) by using the hyperspherical harmonic addition theorem (Wen & Avery, 1985):

where $Z_{k}^{(\alpha)}$ is a normalization term, $\Upsilon_{k}^{h}$ are hyperspherical harmonics³³3We use a real basis of hyperspherical harmonics to avoid complex arithmetic. which are $d$-dimensional generalizations of the spherical harmonics (see (Avery, 1989) for their definition) and

Substituting (9) into (8) yields

This is of the form in (3), using

The number of functions in the approximate expansion (which is $r$ in (3)) will henceforth be referred to as the rank of the expansion. The cost of computing the low-rank expansion is the cost of populating the matrices $U$ and $V^{T}$, which is linear in the number of entries in the matrices. To achieve efficiency, we’d like the rank of the expansion to be as small as possible while still achieving the desired level of accuracy. Unfortunately, the rank of (3.3) is still large—proportional to $d^{p}$. In this section, we incorporate the datasets $\mathcal{X},\mathcal{Y}$ into a scheme that efficiently reduces the rank of (3.3).

First, we define

so that (3.3) may be written as

The function $r^{(k)}$ has rank $\frac{p-2k}{2}$,⁴⁴4$\mathcal{T}^{{}^{\prime\prime}}_{k,m,n}$ is nonzero only if $k,m,n$ have the same parity, see Section A. but if we can find an approximation $\tilde{r}^{(k)}\approx r^{(k)}$ with smaller rank and within our error tolerance (when evaluated on the data), then it may be substituted into (3.3) to reduce the overall rank.

We will find this low-rank approximation by efficiently computing a low-rank approximation of the associated $N\times N$ matrix defined by ${(\mathcal{R}^{(k)})_{ij}\coloneqq r^{(k)}(\|x_{i}\|,\|y_{j}\|).}$ Notably, $\mathcal{R}^{(k)}$ can be written as ${\mathcal{R}^{(k)}=X^{(k)}(Y^{(k)})^{T},}$ where

Letting $X^{(k)}=Q^{(k)}_{X}R^{(k)}_{X}$ and $Y^{(k)}=Q^{(k)}_{Y}R^{(k)}_{Y}$ be $QR$ factorizations we have

where $U^{(k)}\Sigma^{(k)}(V^{(k)})^{T}$ is an SVD of the matrix $R^{(k)}_{X}(R^{(k)}_{Y})^{T}$. Using an appropriate error tolerance $\tau$, we may find a lower rank representation of $\mathcal{R}^{(k)}$ by truncating all singular values of $\Sigma^{(k)}$ less than $\tau$. This results in

where

where the notation $U^{(k)}_{:,1:s_{k}}$ refers to the first $s_{k}$ columns of $U^{(k)}$, $\Sigma^{(k)}_{1:s_{k},1:s_{k}}$ means the first $s_{k}$ rows and columns of $\Sigma^{(k)}$, and $s_{k}$ is the index of the smallest singular value in $\Sigma^{(k)}$ greater than $\tau$. The selection of $\tau$ is discussed in Section 3.5.

Error is introduced into our approximation in two ways: the truncation of the Chebyshev expansion in (4) and the data-driven approximation to $\mathcal{R}^{(k)}$ in (14). Suppose we use a threshold of $\tau$ to truncate singuar values/vectors in (13), then the error $\varepsilon$ of the approximation for any pair of points will satisfy

For the first part, note that ${|C_{k}^{(\alpha)}(\cos{\gamma_{ij}})|\leq\binom{k+d-3}{k}}$ (see (Kim et al., 2012) Eq (1.1)). Further, since ${\|A\|_{\infty}\leq\sqrt{N}\|A\|_{2}}$ for any $N\times N$ matrix $A$, we have

For the second part, (Elliott et al., 1987) gives

Therefore, the total error satisfies

where ${\mathcal{C}_{p/2}\coloneqq\sum_{k=0}^{p/2}\binom{k+d-3}{k}/Z_{k}^{(\alpha)}}.$

The algorithm is summarized in Algorithm 1. Note that $p$ is an adaptively chosen parameter so that (4) is sufficiently accurate. The Chebyshev transform takes $\mathcal{O}(p)$ time via the discrete cosine transform. If an initial multiplication of Vandermonde matrices is performed outside the $k$ loop, the $QR$ factorizations can be performed in $\mathcal{O}(Np)$ time each (see, for example, (Brubeck et al., 2021)) and so contribute $\mathcal{O}(Np^{2})$ to the total cost. Forming $\bar{X}^{(k)},\bar{Y}^{(k)}$ takes $\mathcal{O}(Nps_{k})$ each, contributing $\mathcal{O}(Np^{2}s)$ to the total cost, where $s\coloneqq\max_{k}s_{k}$. The SVDs take $\mathcal{O}(p^{3})$ time each and so contribute $\mathcal{O}(p^{4})$ to the total cost.

The total number of harmonics needed by our algorithm will be bounded above by the number of harmonics of order less than or equal to ${p/2}$—in Section B we derive this to be ${\mathbf{H}_{p/2,d}=\binom{\frac{p}{2}+d-1}{d-1}+\binom{\frac{p}{2}+d-2}{d-1}}.$ Each harmonic costs $\mathcal{O}(d)$ time to evaluate per datapoint, hence the cost of computing the necessary harmonics will be $\mathcal{O}(N\mathbf{H}_{p/2,d}d)$. Once they are computed, the cost to populate the $U$ and $V$ matrices will be $\mathcal{O}(Nr)$ where $r$ is the final rank of the approximation. This rank is given by

Thus the total runtime complexity for the factorization is ${\mathcal{O}(N(r+\mathbf{H}_{p/2,d}d+p^{2}s)+p^{4})}$. After the factorization, matrix-vector multiplies may be performed in ${\mathcal{O}(Nr)}$ time.

We have described an additional data-dependent compression step based on matrices formed out of the polynomials within the analytic expansion. We currently make no attempt to perform further compression based on the observed rank of the harmonic components. This is logical when the number of points is high relative to the dimension—matrices formed from the harmonics will have little singular value decay owing to the orthogonality of the harmonic functions. However, if the number of dimensions is large relative to the the number of points, then the points are not space filling and the lack of additional compression based on the harmonics will result in poor efficiency, especially as the number of harmonics grows polynomially with the dimension (see Section C for further discussion).

Vandermonde matrices are notoriously ill-conditioned (Pan, 2016), and the accuracy of the $QR$ factorizations in the current presentation of our additional compression step can suffer when $p$ is large unless this is addressed. In our testing, this degradation only arises for very small length-scale parameters and/or very small error tolerances (notably, outside the range seen in our experiments). One step towards preventing potential conditioning issues could be to arrange the expansion so that Legendre polynomials of $\|x\|,\|y\|$ are used instead of powers of $\|x\|,\|y\|$ (in an analogous way to our replacement of powers of cosine with Gegenbauer polynomials of cosine, see Section A for details).

In many kernel matrix approximation algorithms, a non-negligible cost of generating the low-rank matrices is that matrix entries typically cost $\mathcal{O}(d)$ time to generate. For example, this is the case when a kernel evaluation is required between a point and a benchmark point, as in inducing point methods. On the one hand, our factorization benefits from the $\mathcal{O}(d)$ costs being confined⁵⁵5The Vandermonde-like matrices require only a single dimension-dependent computation of the norms of all datapoints. to the computation of the harmonics, each of which is reused across $s_{k}$ rows where $s_{k}$ is the rank in (14). On the other hand, this depends on an efficient implementation of hyperspherical harmonic computation—although we provide a moderately tuned implementation, this is the subject of further work as there are few available good pre-existing routines for this procedure.

To test the runtime of our algorithm for constructing the factorization, we generate datasets of points drawn from a normal distribution and subsequently scaled so that the maximum norm of a point in a dataset is $1$. Using the Cauchy kernel with lengthscale parameter $\sigma=1$, we vary the number of points and dimensions of this dataset and record the time required by the implementation to create the low-rank approximation. The relative error tolerance is fixed at $\varepsilon=0.005$. Figure 1 vizualizes the results of this experiment—across a variety of dimensions, our factorization requires linear time in the number of points.

To illustrate how efficiently we can tradeoff time and accuracy, we compare our factorization with the popular Nyström method where inducing points are chosen uniformly at random.⁶⁶6Times for the Nyström method include populating the interpolation matrix as well as factoring the inducing point matrix. We use a dataset of 10,000 points in three dimensions generated the same as before and using the same Cauchy kernel. This time, the error tolerance is varied and the factorization time and relative errors are recorded. Figure 2 shows that our factorization is consistently cheaper to produce than the Nyström factorization that achieves the same accuracy. This further translates to faster matrix vector products at the same accuracy. This is due largely to the HDF’s lower-rank approximation, which we investigate next.

To further explore the tradeoff between rank and relative error of our approximation, we conduct experiments where the error tolerance is varied and the rank of the necessary factorization along with the actual final relative error is collected. In this experiment, we use a dataset of 5000 points in five dimensions generated in the same way as above and use the Cauchy kernel $\left((k(r)=\frac{1}{1+(r/\sigma)^{2}}\right)$, Gaussian kernel ($k(r)=e^{-(r/\sigma)^{2}}$), Matérn kernel with $\nu=1.5$ (${k(r)=(1+\sqrt{3}(r/\sigma))e^{-\sqrt{3}r/\sigma}}$) and Matérn kernel with $\nu=2.5$ (${k(r)=(1+\sqrt{5}(r/\sigma)+(5/3)(r/\sigma)^{2})e^{-\sqrt{5}r/\sigma}}$). Figure 3 compares results between our method, the Nyström method, and the SVD (which provides the optimal rank/error tradeoff). In all cases our method has greater accuracy than the Nyström method for the same rank. Additionally, our method shows nearly optimal performance for lower rank approximations. This experiment also highlights the generality of our method—most existing analytic methods cannot be applied to this breadth of kernels and are, therefore, omitted from the experiment.

We also apply our factorization within the training process of kernel ridge regression and evaluate the predictions. Kernel ridge regression (Shawe-Taylor et al., 2004) is a kernel method which involves finding a weight vector $\mathbf{w}$ that solves the minimization problem

where $\mathbf{x}$ are the training points, $\mathbf{y}$ are the training labels, and $\lambda$ is a regularization parameter. The solution to this problem is

Predictions on test data are then performed using ${\hat{\mathbf{y}}=K(\mathbf{x}^{*},\mathbf{x})\mathbf{w}}$, where $\mathbf{x}^{*}$ are the test points and $\hat{\mathbf{y}}$ are our predictions of the test labels.

We use the California housing data set from the UCI Machine Learning Repository (Dua & Graff, 2017), which contains data drawn from the 1990 U.S. Census. The labels are the median house values, and the features are all other attributes, leading to a feature space of $d=8$ dimensions for a dataset of ${N=20,640}$ points. We scale features to be in $[0,1]$, and for each of five trials perform a random shuffle of the data followed by a $\frac{2}{3}$-$\frac{1}{3}$ train-test split.

We apply our factorization to the problem by generating low-rank approximations for the diagonal blocks of the kernel matrix in (16). The diagonal blocks are chosen to be associated with 30 clusters $C_{i}$ of points generated by $k$-means clustering (Lloyd, 1982), so that the $i$th block along the diagonal is given by ${B_{i}\coloneqq K_{C_{i},C_{i}}}$. A host of techniques exist for compressing the off-diagonal blocks (Ying et al., 2004; Ying, 2006; Ambikasaran et al., 2015; Ryan et al., 2021), but we perform no off-diagonal compression for this experiment to simplify the implementation. The associated linear system is then solved using the method of conjugate gradients (Greenbaum, 1997; Shewchuk, 1994), where the matrix-vector multiplies are accelerated by our compression. As a preconditioner we use a block diagonal matrix whose blocks are the inverses of the blocks of the true kernel matrix associated with the clusters.

For comparison, we also applied the Nyström method using the same ranks for the diagonal blocks. We record the relative errors of the diagonal block approximations, as well as the final mean squared errors (MSEs) when compared with the ground truth test labels. Figure 4 visualizes our findings. For all kernels, our method tends to find a more accurate approximation than Nyström for the same rank, as also demonstrated for the synthetic data in Figure 3. Further, owing to the more accurate approximations, our method yields MSEs closer to those found using dense operations.