Fast Low-Rank Kernel Matrix Factorization through Skeletonized Interpolation

TL;DR

This paper introduces a nearly optimal, efficient method for low-rank kernel matrix factorization using skeletonized interpolation, improving computational cost and accuracy for integral equation discretizations.

Contribution

It presents a novel skeletonized interpolation approach that constructs low-rank factorizations with near-optimal cost and stability, outperforming naive algorithms.

Findings

Achieves low-rank factorization at $\\mathcal{O}(nr)$ cost

Demonstrates asymptotic convergence and stability

Performs well in numerical experiments, near optimal rank

Abstract

Integral equations are commonly encountered when solving complex physical problems. Their discretization leads to a dense kernel matrix that is block or hierarchically low-rank. This paper proposes a new way to build a low-rank factorization of those low-rank blocks at a nearly optimal cost of $\mathcal{O}(nr)$ for a $n \times n$ block submatrix of rank r. This is done by first sampling the kernel function at new interpolation points, then selecting a subset of those using a CUR decomposition and finally using this reduced set of points as pivots for a RRLU-type factorization. We also explain how this implicitly builds an optimal interpolation basis for the Kernel under consideration. We show the asymptotic convergence of the algorithm, explain his stability and demonstrate on numerical examples that it performs very well in practice, allowing to obtain rank nearly equal to the optimal rank at a fraction of the cost of the naive algorithm.

Tables

Table 1. Table 1: Notations used in the paper

$𝒦$	The smooth kernel function
$𝒳,𝒴$	The spaces over which $𝒦$ is defined, i.e., $𝒳\times 𝒴$
$x$, $y$	Variables, $x\in 𝒳$, $y\in 𝒴$
$X$, $Y$	The mesh of points over which to approximate $𝒦$, i.e., $X\times Y$
$K$	The kernel matrix, $K=𝒦(X,Y)$, $K_{ij}=𝒦(x_i,y_j)$
$m$, $n$	$m=\left|X\right|$, $n=\left|Y\right|$
$\overline{X}$, $\overline{Y}$	The tensor grids of Chebyshev points
$\hat{X}$, $\hat{Y}$	The subsets of $\overline{X}$ and $\overline{Y}$ output by the algorithm used to build the low-rank approximation
$\overline{m}$, $\overline{n}$	The number of Chebyshev tensor nodes, $\overline{m}=\left|\overline{X}\right|$, $\overline{n}=\left|\overline{Y}\right|$
$r_0$	The “interpolation” rank of $𝒦$, i.e., $r_0=\min (\left|\overline{X}\right|,\left|\overline{Y}\right|)$
$r_1$	The Skeletonized Interpolation rank of $𝒦$, i.e., $r_1=\left|\hat{X}\right|=\left|\hat{Y}\right|$
$r$	The rank of the continuous SVD of $𝒦$
$S(x,\overline{X})$, $T(y,\overline{Y})$	Row vectors of the Lagrange basis functions, based on $\overline{X}$ and $\overline{Y}$ and evaluated at $x$ and $y$, respectively. Each column is one Lagrange basis function.
$\hat{S}(x,\hat{X})$, $\hat{T}(y,\hat{Y})$	Row vectors of Lagrange basis functions, based on $\hat{X}$ and $\hat{Y}$, built using the Skeletonized Interpolation and evaluated at $x$ and $y$, respectively. Each column is one function.
$w_k$, $w_l$	Chebyshev integration weights
$\mathrm{diag}\left({\overline{W}}_X\right)$, $\mathrm{diag}\left({\overline{W}}_Y\right)$	Diagonal matrices of integration weights when integration is done at nodes $\overline{X}$ and $\overline{Y}$

Equations

$$K_{ij}=\mathcal{K}(x_i,y_j)$$

$$\mathcal{K}(x,y)=\frac{1}{\Vert x-y{\Vert }_2}$$

$$a(x)u(x)+{\int }_{\tilde{\mathcal{Y}}}\mathcal{K}(x,y)u(y)\text{d}y=f(x)\forall x\in \tilde{\mathcal{X}}$$

$$a_i{u}_i+\sum _j{K}_{ij}{u}_j=f_i$$

$$K=US{V}^{\top }$$

$$\mathcal{K}(x,y)\approx \sum _i{\kappa }_i(x;y_0){\varphi }_i(y)$$

$$\mathcal{K}(x,y)\approx S(x,\tilde{X})\mathcal{K}(\tilde{X},\tilde{Y})T(y,\tilde{Y}{)}^{\top }=\sum _k\sum _l{S}_k(x)\mathcal{K}({\tilde{x}}_k,{\tilde{y}}_l)T_l(y)$$

$$\mathcal{K}(x,y)\approx \sum _s{\sigma }_s{u}_s(x)v_s(y)$$

$$\mathcal{K}(x,y)=\mathcal{K}(x,\tilde{Y})\mathcal{K}(\tilde{X},\tilde{Y}{)}^{-1}\mathcal{K}(\tilde{X},y)$$

$$X_{\text{piv}}\subset X,Y_{\text{piv}}\subset Y$$

$$\mathcal{K}(X,Y)\approx \mathcal{K}(X,Y_{\text{piv}})\mathcal{K}(X_{\text{piv}},Y_{\text{piv}}{)}^{-1}\mathcal{K}(X_{\text{piv}},Y)$$

$$K_w=\mathrm{diag}({\overline{W}}_X{)}^{1/2}\mathcal{K}(\overline{X},\overline{Y})\mathrm{diag}({\overline{W}}_Y{)}^{1/2}$$

$$K_w{P}_y=Q_y{R}_y$$

$$K_w^{\top }{P}_x=Q_x{R}_x$$

$$\mathcal{K}(X,Y)\approx \mathcal{K}(X,\hat{Y})\mathcal{K}(\hat{X},\hat{Y}{)}^{-1}\mathcal{K}(\hat{X},Y)$$

$$\mathcal{K}(X,Y)\approx S(X,\overline{X})\mathcal{K}(\overline{X},\overline{Y})T(Y,\overline{Y}{)}^{\top }$$

$$\mathcal{K}(x,y)=\sum _{s=1}^{\infty }{\sigma }_s{u}_s(x)v_s(y)$$

$$\mathcal{K}(X,Y)\approx \sum _{s=1}^r{u}_s(X){\sigma }_s{v}_s(Y)=\left[\begin{array}{ccc}{u}_1(X)& \cdots & u_r(X)\end{array}\right]\left[\begin{array}{ccc}{\sigma }_1& & \\ & \ddots & \\ & & {\sigma }_r\end{array}\right]\left[\begin{array}{c}{v}_1^{\top }(Y)\\ ⋮\\ v_r^{\top }(Y)\end{array}\right]$$

$$\mathcal{K}(X,Y)\approx S(X,\overline{X})\mathcal{K}(\overline{X},\overline{Y})T(Y,\overline{Y}{)}^{\top }$$

$$f(x)-S(x,\overline{X})f(\overline{X})=\frac{f^{(m)}(\xi )}{m!}\prod _{j=1}^m(x-{\overline{X}}_j)$$

$$\prod _{j=1}^m(x-{\overline{X}}_j)=x^m-r_{\overline{X}}(x)$$

$$\min_{\{\lambda_{s},v_{s}\}}\max_{\lambda}\Big{\|}\mathscr{K}(x,\lambda)-\sum_{s=1}^{m}\mathscr{K}(x,\lambda_{s})v_{s}(\lambda)\Big{\|}_{2},\qquad v_{s}(\lambda_{t})=\delta_{st}$$

$$\mathcal{K}(\hat{X},\hat{\Lambda }):=\left[\begin{array}{ccc}\mathcal{K}(\hat{X},{\lambda }_1)& \cdots & \mathcal{K}(\hat{X},{\lambda }_m)\end{array}\right]$$

$$\hat{S}(x,\hat{X})=\mathcal{K}(x,\hat{\Lambda })\mathcal{K}(\hat{X},\hat{\Lambda }{)}^{-1}$$

$$\hat{S}(x,\hat{X})\mathcal{K}(\hat{X},{\lambda }_s)=\mathcal{K}(x,{\lambda }_s)$$

$$\hat{S}(x,\hat{X})\mathcal{K}(\hat{X},\lambda )\approx \sum _{s=1}^m\hat{S}(x,\hat{X})\mathcal{K}(\hat{X},{\lambda }_s)v_s(\lambda )=\sum _{s=1}^m\mathcal{K}(x,{\lambda }_s)v_s(\lambda )\approx \mathcal{K}(x,\lambda )$$

$$\mathcal{K}(x,\lambda )=u(x)v(\lambda )$$

$$\hat{S}(x,\hat{X})=\mathcal{K}(x,{\lambda }_1)\mathcal{K}(x_1,{\lambda }_1{)}^{-1}$$

$$\hat{S}(x,\hat{X})\mathcal{K}(\hat{X},\lambda )=u(x)v({\lambda }_1)\frac{1}{u(x_1)v({\lambda }_1)}u(x_1)v(\lambda )=u(x)v(\lambda )$$

$$\mathscr{K}(X,Y)\approx\big{[}\mathscr{K}(X,\widehat{Y})\mathscr{K}(\widehat{X},\widehat{Y})^{-1}\big{]}\;\mathscr{K}(\widehat{X},\widehat{Y})\;\big{[}\mathscr{K}(\widehat{X},\widehat{Y})^{-1}\mathscr{K}(\widehat{X},Y)\big{]}$$

Full text

{AMS}

15-04, 15B99, 45-04, 45A05, 65F30, 65R20

1 Introduction

In this paper, we are interested in the low-rank approximation of kernel matrices, i.e., matrices $K_{ij}$ defined as

[TABLE]

for $x_{i}\in X=\{x_{1},\dots,x_{m}\}\subseteq\mathcal{X}$ and $y_{j}\in Y=\{y_{1},\dots,y_{n}\}\subseteq\mathcal{Y}$ and where $\mathscr{K}$ is a smooth function over $\mathcal{X}\times\mathcal{Y}$. A typical example is when

[TABLE]

and $X\subset\mathcal{X}$ and $Y\subset\mathcal{Y}$ are two well-separated sets of points.

This kind of matrices arises naturally when considering integral equations like

[TABLE]

where the discretization leads to a linear system of the form

[TABLE]

where $K$ is a dense matrix. While this linear system as a whole is usually not low-rank, one can select subsets of points $X\subset\mathcal{X}$ and $Y\subset\mathcal{Y}$ such that $\mathscr{K}$ is smooth over $\mathcal{X}\times\mathcal{Y}$ and hence $\mathscr{K}(X,Y)$ is low-rank. This corresponds to a submatrix of the complete $K$. Being able to efficiently compute a low-rank factorization of such submatrix would lead to significant computational savings. By “smooth” we usually refer to a function with infinitely many continuous derivatives over its domain. Such a function can be well approximated by its interpolant at Chebyshev nodes for instance.

Low-rank factorization means that we seek a factorization of $K=\mathscr{K}(X,Y)$ as

[TABLE]

where $U\in\mathbb{R}^{m\times r},V\in\mathbb{R}^{n\times r}$, $S\in\mathbb{R}^{r\times r}$, and $r$ is the rank. In that factorization, $U$ and $V$ don’t necessarily have to be orthogonal. One way to compute such a factorization is to first compute the matrix $K$ at a cost $\mathcal{O}\!\left(mn\right)$ and then to perform some rank-revealing factorization like SVD, rank-revealing QR or rank-revealing LU at a cost usually proportional to $\mathcal{O}\!\left(mnr\right)$. But, even though the resulting factorization has a storage cost of $\mathcal{O}\!\left((m+n)r\right)$, linear in the size of $X$ and $Y$, the cost would be proportional to $\mathcal{O}\!\left(mn\right)$, i.e., quadratic.

1.1 Notation

In the following, we will denote by $\mathscr{K}$ a function over $\mathcal{X}\times\mathcal{Y}$. $X$ and $Y$ are finite sequences of vectors such that $X\subset\mathcal{X}$ and $Y\subset\mathcal{Y}$ and $\mathscr{K}(X,Y)$ denotes the matrix $K_{ij}=\mathscr{K}(x_{i},y_{j})$. Small-case letters $x$ and $y$ denote arbitrary variables, while capital-case letters $\overline{X}$, $\widehat{X}$, $\widecheck{X}$, $\widetilde{X}$ denotes sequences of vectors. We denote matrices like $A(X,Y)$ when the rows refer to the set $X$ and the columns to the set $Y$. Table 1 summarizes all the symbols used in this paper.

1.2 Previous work

The problem of efficiently solving (1) has been extensively studied in the past. As indicated above, discretization often leads to a dense matrix $K_{ij}$. Hence, traditional techniques such as the LU factorization cannot be applied because of their $\mathcal{O}\!\left(n^{3}\right)$ time and even $\mathcal{O}\!\left(n^{2}\right)$ storage complexity. The now traditional method used to deal with such matrices is to use the fact that they usually present a (hierarchically) low-rank structure, meaning we can represent the matrix as a hierarchy of low-rank blocks. The Fast Multipole Method (FMM) [28, 12, 2] takes advantage of this fact to accelerate computations of matrix-vector products $Kv$ and one can then couple this with an iterative method. More recently, [10] proposed a kernel-independent FMM based on interpolation of the kernel function.

Other techniques compute explicit low-rank factorization of blocks of the kernel matrix through approximation of the kernel function. The Panel Clustering method [18] first computes a low-rank approximation of $\mathscr{K}(x,y)$ as

[TABLE]

by Taylor series and then uses it to build the low-rank factorization.

Bebendorf and Rjasanow proposed the Adaptive Cross Approximation [4], or ACA, as a technique to efficiently compute low-rank approximations of kernel matrices. ACA has the advantage of only requiring to evaluate rows or columns of the matrix and provides a simple yet very effective solution for smooth kernel matrix approximations. However, it can have convergence issues in some situations (see for instance [7]) if it cannot capture all necessary information to properly build the low-rank basis and lacks convergence guarantees.

In the realm of analytic approximations, [31] (and similarly [6], [7], [10] and [30] in the Fourier space) interpolate $\mathscr{K}(x,y)$ over $\mathcal{X}\times\mathcal{Y}$ using classical interpolation methods (for instance, polynomial interpolation at Chebyshev nodes in [10]), resulting in expressions like

[TABLE]

where $S$ and $T$ are Lagrange interpolation basis functions. Those expressions can be further recompressed by performing a rank-revealing factorization on the node matrix $\mathscr{K}(\widetilde{X},\widetilde{Y})$, for instance using SVD [10] or ACA [7]. Furthermore, [31] takes the SVD of a scaled $\mathscr{K}(\widetilde{X},\widetilde{Y})$ matrix to further recompress the approximation and obtain an explicit expression for $u_{r}$ and $v_{r}$ such that

[TABLE]

where $\{u_{s}\}_{s}$ and $\{v_{s}\}_{s}$ are sequences of orthonormal functions in the usual $L_{2}$ scalar product.

Bebendorf [3] builds a low-rank factorization of the form

[TABLE]

where the nodes $\widetilde{X}$ and $\widetilde{Y}$ are interpolation nodes of an interpolation of $\mathscr{K}(x,y)$ built iteratively. Similarly, in their second version of the Hybrid cross approximation algorithm, Börm and Grasedyck [7] propose applying ACA to the kernel matrix evaluated at interpolation nodes to obtain pivots $\widetilde{X}_{i}$, $\widetilde{Y}_{j}$, and implicitly build an approximation of the form given in Eq. 2. Both those algorithms resemble our approach in that they compute pivots $\widetilde{X},\widetilde{Y}$ in some way and then use Eq. 2 to build the low-rank approximation. In contrast, our algorithm uses weights, and has stronger accuracy guarantees. We highlight those differences in section 5.

Our method inserts itself amongst those low-rank kernel factorization techniques. However, with the notable exception of ACA, those methods often either rely on analytic expressions for the kernel function (and are then limited to some specific ones), or have suboptimal complexities, i.e., greater than $\mathcal{O}\!\left(nr\right)$. In addition, even though we use interpolation nodes, it is worth noting that our method differs from interpolation-based algorithm as we never explicitly build the $S(x,\widetilde{X})$ and $T(y,\widetilde{Y})$ matrices containing the basis functions. We merely rely on their existence.

$\mathcal{H}$-matrices [16, 17, 15] are one way to deal with kernel matrices arising from boundary integral equations that are Hierarchically Block Low-Rank. The compression criterion (i.e., which blocks are compressed as low-rank and which are not) leads to different methods, usually denoted as strongly-admissible (only compress well-separated boxes) or weakly-admissible (compress adjacent boxes as well). In the realm of strongly-admissible $\mathcal{H}$-matrices, the technique of Ho & Ying [22] as well as Tyrtyshnikov [29] are of particular interest for us. They use Skeletonization of the matrix to reduce storage and computation cost. In [22], they combine Skeletonization and Sparsification to keep compressing blocks of $\mathcal{H}$-matrices. [29] uses a somewhat non-traditional Skeletonization technique to also compress hierarchical kernel matrices.

Finally, extending the framework of low-rank compression, [9] uses tensor-train compression to re-write $\mathscr{K}(X,Y)$ as a tensor with one dimension per coordinate, i.e., $\mathscr{K}(x_{1},\dots,x_{d},y_{1},\dots,y_{d})$ and then compress it using the tensor-train model.

1.3 Contribution

1.3.1 Overview of the method

In this paper, we present a new algorithm that performs this low-rank factorization at a cost proportional to $\mathcal{O}\!\left(m+n\right)$. The main advantages of the method are as follows:

The complexity of our method is $\mathcal{O}\!\left(r(m+n)\right)$ (in terms of kernel function $\mathscr{K}$ evaluations) where $r$ is the target rank. 2. 2.

The method is robust and accurate, irrespective of the distribution of points $x$ and $y$. 3. 3.

We can prove both convergence and numerical stability of the resulting algorithm. 4. 4.

The method is very simple and relies on well-optimized BLAS3 (GEMM) and LAPACK (RRQR, LU) kernels.

Consider the problem of approximating $\mathscr{K}(x,y)$ over the mesh $X\times Y$ with $X\in\mathcal{X}$ and $Y\in\mathcal{Y}$. Given the matrix $\mathscr{K}(X,Y)$, one possibility to build a low-rank factorization is to do a rank-revealing LU. This would lead to the selection of

[TABLE]

called the “pivots”, and the low-rank factorization would then be given by

[TABLE]

In practice however, this method may become inefficient as it requires assembling the matrix $\mathscr{K}(X,Y)$ first.

In this paper, we propose and analyze a new method to select the “pivots” outside of the sets $X$ and $Y$. The key advantage is that this selection is independent from the sets $X$ and $Y$, hence the reduced complexity. Let us consider the case where $\mathcal{X},\mathcal{Y}=[-1,1]^{d}$. We will keep this assumption throughout this paper. Then, within $[-1,1]^{d}$, one can build tensor grids of Chebyshev points $\overline{X},\overline{Y}$ and associated integration weights $\overline{W}_{X},\overline{W}_{Y}$ and then consider the matrix

[TABLE]

Denote $r_{0}=\min(|\overline{X}|,|\overline{Y}|)$. Based on interpolation properties, we will show that this matrix is closely related to the continuous kernel $\mathscr{K}(x,y)$. In particular, they share a similar spectrum. Then, we select the sets $\widehat{X}\subset\overline{X}$, $\widehat{Y}\subset\overline{Y}$ by performing strong rank-revealing QRs [13] over, respectively, $K_{w}^{\top}$ and $K_{w}$ (this is also called a CUR decomposition):

[TABLE]

[TABLE]

and build $\widehat{X}$ by selecting the elements of $P_{x}$ associated to the largest rows of $R_{x}$ and similarly for $\widehat{Y}$ (if they differ in size, extend the smallest). We denote the rank of this factorization $r_{1}=|\widehat{X}|=|\widehat{Y}|$, and in practice, we observe that $r_{1}\approx r_{SVD}$, where $r_{SVD}$ is the rank the truncated SVD of $\mathscr{K}(X,Y)$ would provide. The resulting factorization is

[TABLE]

Note that, in this process, at no point did we built any Lagrange basis function associated with $\overline{X}$ and $\overline{Y}$. We only evaluate the kernel $\mathscr{K}$ at $\overline{X}\times\overline{Y}$.

This method appears to be very efficient in selecting sets $\widehat{X}$ and $\widehat{Y}$ of minimum sizes. Indeed, instead, one could aim for a simple interpolation of $\mathscr{K}(x,y)$ over both $\mathcal{X}$ and $\mathcal{Y}$ separately. For instance, using the regular polynomial interpolation at Chebyshev nodes $\overline{X}$ and $\overline{Y}$, it would lead to a factorization of the form

[TABLE]

In this expression, we collect the Lagrange basis functions (each one associated to a node of $\overline{X}$) evaluated at $X$ in the columns of $S(X,\overline{X})$ and similarly for $T(Y,\overline{Y})$. This provides a robust way of building a low-rank approximation. The rank $r_{0}=\min(|\overline{X}|,|\overline{Y}|)$, however, is usually much larger than the true rank $r_{SVD}$ and than $r_{1}$ (given a tolerance). Note that even if those factorizations can always be further recompressed to a rank $\approx r_{SVD}$, they incur a high upfront cost because of the rank $r_{0}\gg r_{SVD}$. See subsubsection 1.3.4 for a discussion about this.

1.3.2 Distinguishing features of the method

Since there are many methods that resemble our approach, we point out its distinguishing features. The singular value decomposition (SVD) offers the optimal low-rank representation in the 2-norm. However, its complexity scales like $\mathcal{O}\!\left(n^{3}\right)$. In addition, we will show that the new approach is negligibly less accurate than the SVD in most cases.

The rank-revealing QR and LU factorization, and methods of random projections [19], have a reduced computational cost of $\mathcal{O}\!\left(n^{2}r\right)$, but still scale quadratically with $n$.

Methods like ACA [4], the rank-revealing LU factorization with rook pivoting [11] and techniques that randomly sample from columns and rows of the matrix scale like $\mathcal{O}\!\left(nr\right)$, but they provide no accuracy guarantees. In fact, counterexamples can be found where these methods fail. In contrast, our approach relies on Chebyshev nodes, which offers strong stability and accuracy guarantees. The fact that new interpolation points, $\overline{X}$ and $\overline{Y}$, are introduced (the Chebyshev nodes) in addition to the existing points in $X$ and $Y$ is one of the key elements.

Analytical methods are available, like the fast multipole method, etc., but they are limited to specific kernels. Other techniques, which are more general, like Taylor expansion and Chebyshev interpolation [10], have strong accuracy guarantees and are as general as the method presented. However, their cost is much greater; in fact, the difference in efficiency is measured directly by the reduction from $r_{0}$ to $r_{1}$ in our approach.

1.3.3 Low-rank approximation based on SVD and interpolation

Consider the kernel function $\mathscr{K}$ and its singular value decomposition [26, theorem VI.17]:

Theorem 1.1** (Singular Value Decomposition).**

Suppose $\mathscr{K}:[-1,1]^{d}\times[-1,1]^{d}$ is square integrable. Then there exist two sequences of orthogonal functions $\{u_{i}\}_{i=1}^{\infty}$ and $\{v_{i}\}_{i=1}^{d}$ and a non-increasing sequence of non-negative real number $\{s_{i}\}_{i=1}^{\infty}$ such that

[TABLE]

As one can see, under relatively mild assumptions, any kernel function can be expanded into a singular value decomposition. Hence from any kernel function expansion we find a low-rank decomposition for the matrix $\mathscr{K}(X,Y)$ (which is not the same as the matrix SVD):

[TABLE]

where the sequence $\{s_{i}\}_{i=1}^{\infty}$ was truncated at an appropriate index $r$. As a general rule of thumb, the smoother the function $\mathscr{K}(x,y)$, the faster the decay of the $\sigma_{s}$’s and the lower the rank.

If we use a polynomial interpolation method with Chebyshev nodes, we get a similar form:

[TABLE]

The interpolation functions $S(x,\overline{X})$ and $T(y,\overline{Y})$ have strong accuracy guarantees, but the number of terms required in the expansion is $r_{0}\gg r\approx r_{1}$. This is because Chebyshev polynomials are designed for a broad class of functions. In contrast, the SVD uses basis functions $u_{s}$ and $v_{s}$ that are optimal for the chosen $\mathscr{K}$.

1.3.4 Optimal interpolation methods

We will now discuss a more general problem, then derive our algorithm as a special case. Let’s start with understanding the optimality of the Chebyshev interpolation. With Chebyshev interpolation, $S(x,\overline{X})$ and $T(y,\overline{Y})$ are polynomials. This is often considered one of the best (most stable and accurate) ways to interpolate smooth functions. We know that for general polynomial interpolants we have:

[TABLE]

If we assume that the derivative $f^{(m)}(\xi)$ is bounded, we can focus on finding interpolation points such that

[TABLE]

is minimal, where $r_{\overline{X}}(x)$ is a degree $m-1$ polynomial. Since we are free to vary the interpolation points $\overline{X}$, then we have $m$ parameters (the location of the interpolation points) and $m$ coefficients in $r_{\overline{X}}$. By varying the location of the interpolation points, we can recover any polynomial $r_{\overline{X}}$. Chebyshev points are known to solve this problem optimally. That is, they lead to an $r_{\overline{X}}$ such that $\max_{x}|x^{m}-r_{\overline{X}}(x)|$ is minimal.

Chebyshev polynomials are a very powerful tool because of their generality and simplicity of use. Despite this, we will see that this can be improved upon with relatively minimal effort. Let’s consider the construction of interpolation formulas for a family of functions $\mathscr{K}(x,\lambda)$, where $\lambda$ is a parameter. We would like to use the SVD, but, because of its high computational cost, we rely on the cheaper rank-revealing QR factorization (RRQR, a QR algorithm with column pivoting). RRQR solves the following optimization problem:

[TABLE]

where the 2-norm is computed over $x$—in addition RRQR produces an orthogonal basis for $\{\mathscr{K}(x,\lambda_{s})\}_{s}$ but this is not needed in the current discussion. The vector space span$\{\mathscr{K}(x,\lambda_{s})\}_{s=1,\ldots,m}$ is close to span$\{u_{s}\}_{s=1,\ldots,m}$ [see Eq. 4], and the error can be bounded by $\sigma_{m+1}$.

Define $\widehat{\Lambda}=\{\lambda_{1},\dots,\lambda_{m}\}$. From there, we identify a set of $m$ interpolation nodes $\widehat{X}$ such that the square matrix

[TABLE]

is as well conditioned as possible. We now define our interpolation operator as

[TABLE]

By design, this operator is exact on $\mathscr{K}(x,\lambda_{s})$:

[TABLE]

It is also very accurate for $\mathscr{K}(x,\lambda)$ since

[TABLE]

With Chebyshev interpolation, $S(x,\overline{X})$ is instead defined using order $m-1$ polynomial functions.

A special case that illustrates the difference between SI and Chebyshev, is with rank-1 kernels:

[TABLE]

In this case, we can pick any $x_{1}$ and $\lambda_{1}$ such that $\mathscr{K}(x_{1},\lambda_{1})\neq 0$, and define $\widehat{X}=\{x_{1}\}$ and

[TABLE]

[TABLE]

SI is exact using a single interpolation point $x_{1}$. An interpolation using Chebyshev polynomials would lead to errors, for any expansion order (unless $u$ is fortuitously a polynomial).

So, one of the key differences between SI and Chebyshev interpolation is that SI uses, as basis for its interpolation, a set of nearly optimal functions that approximate the left singular functions of $\mathscr{K}$, rather than generic polynomial functions.

1.3.5 Proposed method

In this paper, we use the framework from subsubsection 1.3.4 to build an interpolation operator for the class of functions $\mathscr{K}(x,y)$, which we view as a family of functions of $x$ parameterized by $y$ (and vice versa to obtain a symmetric interpolation method). The approximation (Eq. 3) can be rewritten

[TABLE]

and by comparing with Eq. 6, we recognize the interpolation operators:

[TABLE]

These interpolation operators are nearly optimal; because of the way these operators are constructed we call the method “Skeletonized Interpolation.” The sets $\widehat{X}$ and $\widehat{Y}$ are the minimal sets such that if we sample $\mathscr{K}$ at these points we can interpolate $\mathscr{K}$ at any other point with accuracy $\epsilon$. In particular, $\widehat{X}$ and $\widehat{Y}$ are much smaller than their Chebyshev-interpolant counterparts $\overline{X}$ and $\overline{Y}$ and their size, $r_{1}$, is very close to $r$ in Eq. 5. The approach we are proposing produces nearly-optimal interpolation functions for our kernel, instead of generic polynomial functions.

Note that none of the previous discussions explains why the proposed scheme is stable; the inverse $\mathscr{K}(\widehat{X},\widehat{\Lambda})^{-1}$ as well as $\mathscr{K}(\widehat{X},\widehat{Y})^{-1}$ in Eq. 3 could become troublesome numerically. We will explain in detail in section 3 why this is not an issue numerically, and we explore the connection with interpolation in more detail in section 4.

1.3.6 Organization of the paper

This paper is organized as follows. In section 2, we present the algorithm in detail and present some theoretical results about its convergence. In section 3, we discuss its numerical stability and in section 4 we revisit the interpolation interpretation on a simple example. Finally, section 5 illustrates the algorithm on more complex geometries, compares its accuracy with other classical algorithms and presents computational complexity results.

2 Skeletonized Interpolation

2.1 The algorithm

Algorithm 1 provides the high-level version of the algorithm. It consists of 3 steps:

• •

Build grids $\overline{X}$ and $\overline{Y}$, tensor grids of Chebyshev nodes. Over $[-1,1]$ in 1D, the $\overline{m}$ Chebyshev nodes of the first kind are defined as

[TABLE]

In higher dimensions, they are defined as the tensor product of one-dimensional grids. The number of points in every dimension should be such that

[TABLE]

provides an $\delta$ uniform approximation over $[-1,1]^{d}\times[-1,1]^{d}$ of $\mathscr{K}(x,y)$. Denote

[TABLE]

• •

Recompress the grid by performing a strong rank-revealing QR factorization [13] of

[TABLE]

and its transpose, up to accuracy $\epsilon$. This factorization is also named CUR decomposition [24, 8]. While our error estimates only hold for strong rank-revealing QR factorizations, in practice, a simple column-pivoted QR factorization based on choosing columns with large norms works as well. In the case of Chebyshev nodes of the first kind in 1D over $[-1,1]$ the integration weights are given by

[TABLE]

The weights in $d$ dimensions are the products of the corresponding weights in 1D, and the $\operatorname{diag}(\overline{W}_{X})$ and $\operatorname{diag}(\overline{W}_{Y})$ matrices are simply the diagonal matrices of the integration weights. Denote

[TABLE]

In case the sets $\widehat{X}$ and $\widehat{Y}$ output by those RRQR’s are of slightly different size (which we rarely noticed in our experiments), extend the smallest to have the same size as the largest.

• •

Given $\widehat{X}$ and $\widehat{Y}$, the low-rank approximation is given by

[TABLE]

of rank $r_{1}\approx r_{SVD}$.

2.2 Theoretical Convergence

2.2.1 Overview

In this section, we prove that the error made during the RRQR is not too much amplified when evaluating the interpolant. We first recall that

From interpolation properties,

[TABLE]

where $T$ and $S$ are small matrices (i.e., bounded by logarithmic factors in $r_{0}$) and $E_{\textrm{INT}}=\mathcal{O}\!\left(\delta\right)$. 2. 2.

From the strong RRQR properties,

[TABLE]

where $\hat{K}_{w}$ has a spectrum similar to that of $K_{w}$ (up to a small polynomial), $\widehat{S}$ and $\widehat{T}$ are bounded by a small polynomial, and $E_{\textrm{QR}}=\mathcal{O}\!\left(\epsilon\right)$.

Then, by combining those two facts and assuming $\delta<\epsilon$, one can show

First, that the interpolation operators are bounded,

[TABLE]

where $p$ is a small polynomial. 2. 2.

Second, that the error $\epsilon$ made in the RRQR is not too much amplified, i.e.,

[TABLE]

where $p^{\prime}$ is another small polynomial.

Finally, if one assume that $\sigma_{i}(K_{w})$ decays exponentially fast, so does $\epsilon$ and the resulting approximation in Eq. 10 converges.

In the following, we present the main lemmas (some proofs are relocated in the appendix for brevity) leading to the above result.

2.2.2 Interpolation-related results

We first consider the interpolation itself. Consider $\overline{X}$ and $\overline{Y}$, constructed such as

[TABLE]

Lemma 2.1** (Interpolation at Chebyshev Nodes).**

$\forall x\in\mathcal{X}$* and $\overline{X}$ tensor grids of Chebyshev nodes of the first kind,*

[TABLE]

where $\mathcal{X}\subset\mathbb{R}^{d}$. In addition, the weights, collected in the weight matrix $\operatorname{diag}(\overline{W}_{X})$ are such that

[TABLE]

*where ${\overline{m}}=|\overline{X}|$. *

2.2.3 Skeletonization results

We now consider the skeletonization step of the algorithm performed through the two successive rank-revealing QR factorizations.

Rank-Revealing QR factorizations

Let us first recall what a rank-revealing QR factorization is. Given a matrix $A\in\mathbb{R}^{m\times n}$, one can compute a rank-revealing QR factorization [11] of the form

[TABLE]

where $\Pi$ is a permutation matrix, $Q$ an orthogonal matrix and $R$ a triangular matrix. Both $R$ and $Q$ are partitioned so that $Q_{1}\in\mathbb{R}^{m\times k}$ and $R_{11}\in\mathbb{R}^{k\times k}$. If $\|R_{22}\|\approx\varepsilon$, this factorization typically indicates that $A$ has an $\varepsilon$-rank of $k$. The converse, however, is not necessarily true [11] in general.

From there, one can also write

[TABLE]

where $T$ is the interpolation operator, $A_{1}$ a set of $k$ columns of $A$ and $E$ the approximation error. This approximation can be achieved by a simple column-pivoted QR algorithm [11]. This algorithm, however, is not guaranteed to always work (i.e., even if $A$ has rapidly decaying singular values, this rank-revealing factorization may fail to exhibit it).

A strong rank-revealing QR, however, has more properties. It has been proven [13, 8] that one can compute in $\mathcal{O}\!\left(mn^{2}\right)$ a rank-revealing QR factorization that guarantees

[TABLE]

where $q_{1}$ and $q_{2}$ are two small polynomials (with fixed constants and degrees). The existence of this factorization is a crucial part of our argument. Using the interlacing property of singular values [11], this implies that we now have both lower and upper bounds on the singular values of $A_{1}$

[TABLE]

From Eq. 11 we can directly relate the error $E$ and $\sigma_{k+1}$ from

[TABLE]

Finally, given a matrix $A$, one can apply the above result to both its rows and columns, leading to a factorization

[TABLE]

with the same properties as detailed above.

Skeletonized Interpolation

We can now apply this results to the $K_{w}$ and $\hat{K}_{w}$ matrices.

Lemma 2.2** (CUR Decomposition of $K_{w}$).**

The partition $\overline{X}=\widehat{X}\cup\widecheck{X}$, $\overline{Y}=\widehat{Y}\cup\widecheck{Y}$ of Algorithm 1 is such that there exist $\widecheck{S}$, $\widecheck{T}$, $E_{\textrm{QR}}(\overline{X},\overline{Y})$ matrices and a slowly-growing polynomial $p(r_{0},r_{1})$ such that

[TABLE]

and where

[TABLE]

Finally, we have

[TABLE]

Proof 2.3**.**

The first three results are direct applications of [8, theorem 3 and remark 5] as explained in the previous paragraph. The last result follows from the properties of the strong rank-revealing QR:

[TABLE]

The first inequality follows from $\sigma_{r_{1}}(K_{w})\leq\sigma_{r_{1}}(\hat{K}_{w})p(r_{0},r_{1})$ (Eq. 12), the second from $\sigma_{r_{1}}(K_{w})\geq\sigma_{r_{1}+1}(K_{w})$ (by definition of singular values) and the last from $\sigma_{r_{1}+1}(K_{w})^{-1}\leq p(r_{0},r_{1})\epsilon^{-1}$ (Eq. 13).

Finally, a less obvious result

Lemma 2.4**.**

There exist a polynomial $q(r_{0},r_{1})$ such that for any $x\in\mathcal{X},y\in\mathcal{Y}$,

[TABLE]

We provide the proof in the appendix; the key ingredient is simply that $\|\hat{K}_{w}^{-1}\|_{2}\leq p(r_{0},r_{1})^{2}\epsilon^{-1}$ from the RRQR properties; hence $\hat{K}_{w}$ is ill-conditioned, but not arbitrarily. Its condition number grows like $\epsilon^{-1}$. Then, when multiplied by quantities like $\epsilon$ or $\delta\ll\epsilon$, the factors cancel out and the resulting product can be bounded.

2.2.4 Link between the node matrix and the continuous SVD

In this section, we link the continuous SVD and the spectrum (singular values) of the matrix $\operatorname{diag}(\overline{W}_{X})^{1/2}K_{w}\operatorname{diag}(\overline{W}_{Y})^{1/2}$. This justifies the use of the weights.

For the sake of simplicity, consider the case where interpolation is performed at Gauss-Legendre nodes $\overline{X},\overline{Y}$ with the corresponding integration weights $\overline{W}_{X},\overline{W}_{Y}$. (A more complete explanation can be found in [31].)

Take the classical discrete SVD of $K_{w}$,

[TABLE]

We then have

[TABLE]

Then, denote the sets of new basis functions

[TABLE]

The key is to note that those functions are orthonormal. Namely, for $\overline{u}$,

[TABLE]

The same result holds for $\overline{v}$. This follows from the fact that a Gauss-Legendre quadrature rule with $n$ points can exactly integrate polynomials up to degree $2n-1$. This shows that we are implicitly building a factorization

[TABLE]

where the approximation error is bounded by the interpolation error $E_{\textrm{INT}}$ and where the sets of basis functions are orthogonal.

Assume now that the kernel $\mathscr{K}$ is square-integrable over $[-1,1]^{d}\times[-1,1]^{d}$. This is called a Hilbert-Schmidt kernel [27, Lemma 8.20]. This implies that the associated linear operator is compact [27, Theorem 8.83]. $\overline{\mathscr{K}}$ is compact as well since it is finite rank [27, Theorem 8.80]. Given the fact that $|E_{\textrm{INT}}(x,y)|\leq\delta$ for all $x,y$, $\|\mathscr{K}-\overline{\mathscr{K}}\|_{L_{2}}\leq C\delta$ for some $C$ and hence, by compactness of both operators [14, Corollary 2.2.14],

[TABLE]

for some $C>0$. Then, from the above discussion, we clearly have $\sigma_{i}(K_{w})=\sigma_{i}(\overline{\mathscr{K}})+\mathcal{O}\!\left(\delta\right)$ and hence

[TABLE]

This result only formally holds for Gauss-Legendre nodes and weights. However, this motivates the use of integration weights in the case of Chebyshev as well.

2.2.5 Convergence of the Skeletonized interpolation

We now present the main result of this paper:

Theorem 2.5** (Convergence of Skeletonized Interpolation).**

If $\widehat{X}$ and $\widehat{Y}$ are constructed following Algorithm 1, then there exist a polynomial $r(r_{0},r_{1})$ such that for any $x\in\mathcal{X}$ and $y\in\mathcal{Y}$,

[TABLE]

The key here is that the error incurred during the CUR decomposition, $\epsilon$ is amplified by, at most, a polynomial of $r_{0}$ and $r_{1}$. Hence, Theorem 2.5 indicates that if the spectrum decays fast enough (i.e., if $\epsilon\to 0$ when $r_{0},r_{1}\to\infty$ faster than $r(r_{0},r_{1})$ grows), the proposed approximation should converge to the true value of $\mathscr{K}(x,y)$.

What is simply left is then linking $\epsilon$, $r_{0}$, and $r_{1}$. We have, from the CUR properties,

[TABLE]

which implies

[TABLE]

Then, following the discussion from subsubsection 2.2.4, we expect

[TABLE]

Hence, if $\mathscr{K}$ has rapidly-decaying singular values, so does $K_{w}$. Assuming the singular values of $K_{w}$ decay exponentially fast, i.e.,

[TABLE]

we find

[TABLE]

as $r_{0},r_{1}\to\infty$, or alternatively, as $\epsilon\to 0$.

3 Numerical stability

3.1 The problem

The previous section indicates that, at least theoretically, we can expect convergence as $\epsilon\to 0$. However, the factorization

[TABLE]

seems to be numerically challenging to compute. Indeed, as we showed in the previous section, we can only really expect at best $\|\hat{K}_{w}^{-1}\|_{2}=\mathcal{O}\!\left(\epsilon^{-1}\right)$ which indicates that, roughly,

[TABLE]

i.e., the condition number grows with the desired accuracy, and convergence beyond a certain threshold (like $10^{-8}$ in double-precision) seems impossible. Hence, we can reasonably be worried about the numerical accuracy of computing Eq. 15 even with a stable algorithm.

Note that this is not a pessimistic upper bound; by construction, $\hat{K}_{w}$ really is ill-conditioned, and experiments show that solving linear systems $\hat{K}_{w}x=b$ with random right-hand sides is numerically challenging and leads to errors of the order $\epsilon^{-1}$.

3.2 Error Analysis

Consider Eq. 15 and let for simplicity

[TABLE]

In this section, our goal is to show why one can expect this formula to be accurately computed if one uses backward stable algorithms. As proved in section 2, we have the following bounds on the interpolation operators

[TABLE]

for some polynomial $p$. The key is that there is no $\epsilon^{-1}$ in this expression. Those bounds essentially follow from the guarantees provided by the strong rank-revealing QR algorithm.

Now, let’s compute the derivative of $K_{x}K^{-1}K_{y}$ with respect to $K_{x}$, $K$ and $K_{y}$ [25]:

[TABLE]

Then, consider perturbing $K_{x}$, $K$, $K_{y}$ by $\varepsilon$ (assume all matrices are of order $\mathcal{O}\!\left(1\right)$ for the sake of simplicity), i.e., let $\delta K_{x}$, $\delta K_{y}$ and $\delta K$ be perturbations of $K_{x}$, $K_{y}$ and $K$, respectively, with

[TABLE]

Then, using the above derivative as a first order approximation, we can write

[TABLE]

We see that the computed result is independent of the condition number of $K=\mathscr{K}(\widehat{X},\widehat{Y})$ and depends on $p(r_{0},r_{1})$ only.

Assume now that we are using backward stable algorithms in our calculations [20]. We then know that the computed result is the result of an exact computation where the inputs have been perturbed by $\varepsilon$. The above result indicates that the numerical result (with roundoff errors) can be expected to be accurate up to $\varepsilon$ times a small polynomial, hence stable.

4 Skeletonized Interpolation as a new interpolation rule

As indicated in the introduction, one can rewrite

[TABLE]

where we recognize two new “cross-interpolation” (because they are build by considering both the $\mathcal{X}$ and $\mathcal{Y}$ space) operators $\widehat{S}(x,\widehat{X})=\mathscr{K}(x,\widehat{Y})\mathscr{K}(\widehat{X},\widehat{Y})^{-1}$ and $\widehat{T}(y,\widehat{Y})=\mathscr{K}(\widehat{X},y)^{\top}\mathscr{K}(\widehat{X},\widehat{Y})^{-\top}$. In this notation, each column of $\widehat{S}(x,\widehat{X})$ and $\widehat{T}(y,\widehat{Y})$ is a Lagrange function associated to the corresponding node in $\widehat{X}$ or $\widehat{Y}$ and evaluated at $x$ or $y$, respectively.

This interpretation is interesting as it allows to “decouple” $x$ and $y$ and analyze them independently. In particular, one can look at the quality of the interpolation of the basis functions $u_{k}(x)$ and $v_{k}(y)$ using $\widehat{S}$ and $\widehat{T}$. Indeed, if this is accurate, it is easy to see that the final factorization is accurate. Indeed,

[TABLE]

To illustrate this, let us consider a simple 1-dimensional example. Let $x,y\in[-1,1]$ and consider

[TABLE]

Then approximate this function up to $\epsilon=10^{-10}$, and obtain a factorization of rank $r$.

Fig. 1 illustrates the 4th Lagrange basis function in $x$, i.e., $\widehat{S}(x,\widehat{X})_{4}$ and the classical Lagrange polynomial basis function associated with the same set $\widehat{X}$. We see that they are both 1 at $\widehat{X}_{4}$ and 0 at the other points. However, $\widehat{S}(x,\widehat{X})_{4}$ is much more stable and small than its polynomial counterpart. In the case of polynomial interpolation at equispaced nodes, the growth of the Lagrange basis function (or, equivalently, of the Lebesgue constant) is the reason for the inaccuracy and instability.

Fig. 2 shows the effect of interpolating $u_{r}(x)$ using $\widehat{S}(x,\widehat{X})$ as well as using the usual polynomial interpolation at the nodes $\widehat{X}$. We see that $\widehat{S}(x,\widehat{X})$ interpolates very well $u_{r}(x)$, showing indeed that we implicitly build an accurate interpolant, even on the last (least smooth) eigenfunctions. The usual polynomial interpolation fails to capture any feature of $u_{r}$ on the other hand. Note that we could have reached a similar accuracy using interpolation at Chebyshev nodes but only by using many more interpolation nodes.

Finally, Fig. 3 shows how well we approximate the various $r$ eigenfunctions. As one can see, interpolation is very accurate on $u_{1}(x)$, but the error grows for less smooth eigenfunctions. The growth is, roughly, similar to the growth of $\frac{\epsilon}{\sigma_{i}}$. Notice how this is just enough so that the resulting factorization is accurate:

[TABLE]

It is also consistent with the analysis of section 3. This illustrates how the algorithm works: it builds an interpolation scheme that allows for interpolating the various eigenfunctions of $\mathscr{K}$ with just enough accuracy so that the resulting interpolation is accurate up to the desired accuracy.

5 Numerical experiments

In this section we present some numerical experiments on various geometries. We study the quality (how far $r_{1}$ is from the optimal SVD-rank $r$ and how accurate the approximation is) of the algorithm in subsection 5.1 and subsection 5.2. We illustrate in subsection 5.3 the improved guarantees of RRQR and justify the use of weights in subsection 5.4. Finally, subsection 5.5 studies the algorithm computational complexity.

The experiments are done using Julia [5] and the code is sequential. For the strong rank-revealing QR algorithm, we use the LowRankApprox.jl Julia package [21]. The code can be downloaded from https://people.stanford.edu/lcambier/publications.

5.1 Simple geometries

We begin this section with an elementary problem, as depicted on Fig. 4b. In this problem, we consider the usual kernel $\mathscr{K}(x,y)=\|x-y\|_{2}^{-1}$ where $x,y\in\mathbb{R}^{2}$. $\mathcal{X}$ and $\mathcal{Y}$ are two squares of side of length 1, centered at $(0.5,0.5)$ and $(2.5,2.5)$ respectively. Finally, $X$ and $Y$ are two uniform meshes of $50\times 50$ mesh points each, i.e., $n=2500$.

We pick the Chebyshev grids $\overline{X}$ and $\overline{Y}$ using a heuristic based on the target accuracy $\varepsilon$. Namely, we pick the number of Chebyshev nodes in each dimension (i.e., $x_{1}$, $x_{2}$, $y_{1}$ and $y_{2}$) independently (by using the midpoint of $\mathcal{X}$ and $\mathcal{Y}$ as reference), such that the interpolation error is approximately less than $\varepsilon^{3/4}$. This values is heuristic, but performs well for those geometries. Other techniques are possible. This choice is based in part on the observation that the algorithm is accurate even when $\delta>\varepsilon$, i.e., when the Chebyshev interpolation is less accurate than the actual final low-rank approximation through Skeletonized Interpolation.

Consider then Fig. 4a. The $r_{0}$ line indicates the rank ($r_{0}=\min(|\overline{X}|,|\overline{Y}|)$) of the low-rank expansion through interpolation. The $r_{1}$ line corresponds to the rank obtained after the RRQR over $\mathscr{K}(\overline{X},\overline{Y})$ and its transpose, i.e., it is the rank of the final approximation

[TABLE]

Finally, “SVD rank” is the rank one would obtain by truncating the SVD of

[TABLE]

at the appropriate singular value, as to ensure

[TABLE]

Similarly, “RRQR” is the rank a rank-revealing QR on $\mathscr{K}(X,Y)$ would obtain. This is usually slightly suboptimal compared to the SVD. Those two values are there as to illustrate that $r_{1}$ is close to the optimal value.

The conclusion regarding Fig. 4a is that Skeletonized Interpolation is nearly optimal in terms of rank. While the rank obtained by the interpolation is clearly far from optimal, the RRQR over $\mathscr{K}(\overline{X},\overline{Y})$ allows us to find subsets $\widehat{X}\subset\overline{X}$ and $\widehat{Y}\subset\overline{Y}$ that are enough to represent well $\mathscr{K}(X,Y)$, and the final rank $r_{1}$ is nearly optimal compared to the SVD-rank $r$. We also see that the rank of a blind RRQR over $\mathscr{K}(X,Y)$ is higher than the SVD-rank and usually closer—if not identical—to $r_{1}$.

We want to re-emphasize that, in practice, the error of the sets $\overline{X}$, $\overline{Y}$ —i.e., the error of the polynomial interpolation based on $\overline{X}\times\overline{Y}$— can be larger than the required tolerance. If they are large enough, the compressed sets $\widehat{X}$, $\widehat{Y}$ will contain enough information so as to properly interpolate $\mathscr{K}$ and the final error will be smaller than the required tolerance. This is important, as the size of the Chebyshev grid for a given tolerance can be fairly large (as indicated in the introduction, and one of the main motivation of this work), even though the final rank is small.

As a sanity check, Fig. 4c gives the relative error measured in the Frobenius norm

[TABLE]

between $\mathscr{K}(X,Y)$ and its interpolation as a function of the tolerance $\epsilon$.111Choosing the Frobenius norm is not critical—very similar results are obtained in the 2-norm for instance. We see that both lines are almost next to each other, meaning our approximation indeed reaches the required tolerance. This is important as it means that one can effectively control the accuracy.

Finally, Fig. 4b also shows the resulting $\widehat{X}$ and $\widehat{Y}$. It is interesting to notice how the selected points cluster near the close corners, as one could expect since this is the area where the kernel is the least smooth.

We then consider results for the same Laplacian kernel $\mathscr{K}(x,y)=\|x-y\|_{2}^{-1}$ between two plates in 3D (Fig. 5b). We observe overall very similar results as for the previous case on Fig. 5a, where the initial rank $r_{0}$ is significantly decreased to $r_{1}$ while keeping the resulting accuracy close to the required tolerance as Fig. 5c shows. Finally, one can see on Fig. 5b the selected Chebyshev nodes. They again cluster in the areas where smoothness is the worst, i.e., at the closes edges of the plates.

5.2 Comparison with ACA and Random Sampling

We then compare our method with other standard algorithms for kernel matrix factorization. In particular, we compare it with ACA [4] and ’Random CUR’ where one selects, at random, pivots $\widetilde{X}$ and $\widetilde{Y}$ and builds a factorization

[TABLE]

based on those. As we are interested in comparing the quality of the resulting sets of pivots for a given rank, we compare those algorithms for sets $X$ and $Y$ with variable distance between each other and for a fixed tolerance ($\epsilon=10^{-8}$) and kernel ($\mathscr{K}(x,y)=\|x-y\|_{2}^{-1}$). The geometry is two unit-length squares side-by-side with a variable distance between their closest edges.

The comparison is done in the following way:

Given $r_{1}$, build the ACA factorization of rank $r_{1}$ and compute its relative error in Frobenius norm with $\mathscr{K}(X,Y)$; 2. 2.

Given $r_{1}$, build the random CUR factorization by sampling uniformly at random $r_{1}$ points from $X$ and $Y$ to build $\widetilde{X}$, $\widetilde{Y}$. Then, compute its relative error with $\mathscr{K}(X,Y)$.

We then do so for sets of varying distance, and for a given distance, we repeat the experiment 25 times by building $X$ and $Y$ at random within the two squares. This allows to study the variance of the error and to collect some statistics.

Fig. 6 gives the resulting errors in relative Frobenius norm for the 3 algorithm using box-plots of the errors to show distributions. The rectangular boxes represent the distributions from the 25% to the 75% quantiles, with the median in the center. The thinner lines represent the complete distribution, except outliers depicted using large dots. We observe that the ($\widehat{X}$, $\widehat{Y}$) sets based on Chebyshev-SI are, for a common size $r_{1}$, more accurate than the Random or ACA sets. In addition, by design, they lead to more stable factorizations (as they have very small variance in terms of accuracy) while ACA for instance has a higher variance. We also see, as one may expect, that while ACA is still fairly stable even when the clusters get close, random CUR starts having higher and higher variance. This is understandable as the kernel gets less and less smooth as the clusters get close.

Finally, we ran the same experiments with several other kernels (${r^{-2}}$, ${r^{-3}}$, $\log(r)$, $\exp(-r)$, $\exp(-r^{2})$) and observed quantitatively very similar results.

5.3 Stability guarantees provided by RRQR

In Algorithm 1, in principle, any rank-revealing factorization providing pivots could be used. In particular, ACA itself could be used. In this case, this is the HCAII (without the weights) algorithm as described in [7]. However, ACA is only a heuristic: unlike strong rank-revealing factorizations, it can’t always reveal the rank. In particular, it may have issues when some parts of $X$ and $Y$ have strong interactions while others are weakly coupled. To highlight this, consider the following example. It can be extended to many other situations.

Let us use the rapidly decaying kernel

[TABLE]

and the situation depicted in Fig. 7 with $X=\begin{bmatrix}X_{1}&X_{2}\end{bmatrix}$ and $Y=\begin{bmatrix}Y_{1}&Y_{2}\end{bmatrix}$. We note that, formally, $X$ and $Y$ are not well-separated.

Since $\mathscr{K}$ is rapidly decaying and $X_{1}$/$Y_{2}$ (resp. $X_{2}$/$Y_{1}$) are far away, the resulting matrix is nearly block diagonal, i.e.,

[TABLE]

for some small $\varepsilon$. This is a challenging situation for ACA since it will need to sweep through the initial block completely before considering the other one. In practice heuristics can help alleviate the issue; see ACA+ [7] for instance. Those heuristic, however, do not come with any guarantees. Strong RRQR, on the other hand, does not suffer from this and picks optimal nodes in each cluster from the start. It guarantees stability and convergence.

5.4 The need for weights

Another characteristic of Algorithm 1 is the presence of weights. We illustrate here why this is necessary in general. Algorithm 1 uses $\overline{X}$ and $\overline{Y}$ both to select interpolation points (the “columns” of the RRQR) and to evaluate the resulting error (the “rows”). Hence, a non-uniform distribution of points leads to over- or under-estimated $L_{2}$ error and to a biased interpolation point selection. The weights, roughly equal to the (square-root of) the inverse points density, alleviate this effect. This is a somewhat small effect in the case of Chebyshev nodes & weights as the weights have limited amplitudes.

To illustrate this phenomenon more dramatically, consider the situation depicted on Fig. 8. We define $\overline{X}$ and $\overline{Y}$ in the following way. Align two segments of $N$ points, separated by a small interval of length $\epsilon=1/N$ with $N=200$. At the close extremities we insert $25N$ additional points inside small 3D spheres of diameter $\epsilon$. As a result, $|\overline{X}|=|\overline{Y}|=26N=5,200$, and $\overline{X}$, $\overline{Y}$ are strongly non-uniform. We see that the small spheres hold a large number of points in an interval of length $N^{-1}$. As a result, their associated weight should be proportional to $N^{-1/2}$, while the weight for the points on the segments should be proportional to $1$. Then we apply Algorithm 1 with and without weights, and evaluate the error on the segments using $|X|=|Y|=10,000$ equispaced points as a proxy for the $L_{2}$ error.

When using a rank $r_{0}=200$, the CUR decomposition picks only 6 more points on the segments (outside the spheres) for the weighted case compared to the unweighted. However, this is enough to dramatically improve the accuracy, as Fig. 8b shows. Overall, the presence of weights has a large effect, and this shows that in general, one should appropriately weigh the node matrix $K_{w}$ to ensure maximum accuracy.

5.5 Computational complexity

We finally study the computational complexity of the algorithm. It’s important to note that two kinds of operations are involved: kernel evaluations and classical flops. As they may potentially differ in cost, we keep those separated in the following analysis.

The cost of the various parts of the algorithm is the following :

• •

$\mathcal{O}\!\left(r_{0}^{2}\right)$ kernel evaluations for the interpolation, i.e., the construction of $\overline{X}$ and $\overline{Y}$ and the construction of $\mathscr{K}(\overline{X},\overline{Y})$

• •

$\mathcal{O}\!\left(r_{0}^{2}r_{1}\right)$ flops for the RRQR over $\mathscr{K}(\overline{X},\overline{Y})$ and $\mathscr{K}(\overline{X},\overline{Y})^{\top}$

• •

$\mathcal{O}\!\left((m+n)r_{1}\right)$ kernel evaluations for computing $\mathscr{K}(X,\widehat{Y})$ and $\mathscr{K}(\widehat{X},Y)$, respectively (with $m=|X|$ and $n=|Y|$)

• •

$\mathcal{O}\!\left(r_{1}^{3}\right)$ flops for $\mathscr{K}(\widehat{X},\widehat{Y})^{-1}$ (through, say, an LU factorization)

So the total complexity of building the three factor is $\mathcal{O}\!\left((m+n)r_{1}\right)$ kernel evaluations. If $m=n$ and $r_{1}\approx r$, the total complexity is

[TABLE]

Also note that the memory requirements are, clearly, of order $\mathcal{O}\!\left((m+n)r_{1}\right)$.

When applying this low-rank matrix on a given input vector $f(Y)\in\mathbb{R}^{n}$, the cost is

• •

$\mathcal{O}\!\left(r_{1}n\right)$ flops for computing $w_{1}=\mathscr{K}(\widehat{X},Y)f(Y)$

• •

$\mathcal{O}\!\left(r_{1}^{2}\right)$ flops for computing $w_{2}=\mathscr{K}(\widehat{X},\widehat{Y})^{-1}w_{1}$ assuming a factorization of

$\mathscr{K}(\widehat{X},\widehat{Y})$ has already been computed

• •

$\mathcal{O}\!\left(mr_{1}\right)$ flops for computing $w_{3}=\mathscr{K}(X,\widehat{Y})w_{2}$

So the total cost is

[TABLE]

flops if $m=n$ and $r_{1}\approx r$.

To illustrate those results, Fig. 9a shows, using the same setup as in the 2D square example of subsection 5.1, the time (in seconds) taken by our algorithm versus the time taken by a naive algorithm that would first build $\mathscr{K}(X,Y)$ and then perform a rank-revealing QR on it. Time is given as a function of $n$ for a fixed accuracy $\epsilon=10^{-8}$. One should not focus on the absolute values of the timing but rather the asymptotic complexities. In this case, the $\mathcal{O}\!\left(n\right)$ and $\mathcal{O}\!\left(n^{2}\right)$ complexities clearly appear, and our algorithm scales much better than the naive one (or, really, that any algorithm that requires building the full matrix first). Note that we observe no loss of accuracy as $n$ grows. Also note that the plateau at the beginning of the Skeletonized Interpolation curve is all the overhead involved in selecting the Chebyshev points $\overline{X}$ and $\overline{Y}$ using some heuristic. This is very implementation-dependent and could be reduced significantly with a better or more problem-tailored algorithm. However, since this is by design independent of $X$ and $Y$ (and, hence, $n$) it does not affect the asymptotic complexity.

Fig. 9b shows the time as a function of the desired accuracy $\epsilon$, for a fixed number of mesh points $n=10^{5}$. Since the singular values of $\mathscr{K}(X,Y)$ decay exponentially, one has $r\approx\mathcal{O}\!\left(\log\left(\frac{1}{\epsilon}\right)\right)$. The complexity of the algorithm being $\mathcal{O}\!\left(nr\right)$, we expect the time to be proportional to $\log\left(\frac{1}{\epsilon}\right)$. This is indeed what we observe.

Fig. 9c depicts the time as a function of the rank $r$ for a fixed accuracy $\epsilon=10^{-8}$ and number of mesh points $n=10^{5}$. In that case, to vary the rank and keep $\epsilon$ fixed, we change the geometry and observe the resulting rank. This is done by moving the top-right square (see Fig. 4b) towards the bottom-left one (keeping approximately one cluster diameter between them) or away from it (up to 6 diameters). The rank displayed is the rank obtained by the factorization. As expected, the algorithm scales linearly as a function of $r$.

6 Conclusion

In this work, we built a kernel matrix low-rank approximation based on Skeletonized interpolation. This can be seen as an optimal way to interpolate families of functions using a custom basis.

This type of interpolation, by design, is always at least as good as polynomial interpolation as it always requires the minimal number of basis functions for a given approximation error. We proved in this paper the asymptotic convergences of the scheme for kernels exhibiting fast (i.e., faster than polynomial) decay of singular values. We also proved the numerical stability of general Schur-complement types of formulas when using a backward stable algorithm.

In practice, the algorithm exhibits a low computational complexity of $\mathcal{O}\!\left(nr\right)$ with small constants and is very simple to use. Furthermore, the accuracy can be set a priori and in practice, we observe nearly optimal convergence of the algorithm. Finally, the algorithm is completely insensible to the mesh point distribution, leading to more stable sets of “pivots” than Random Sampling or ACA.

Acknowledgements

We would like to thank Cleve Ashcraft for his ideas and comments on the paper, as well as the anonymous reviewer for his careful reading and pertinent suggestions that greatly improved the paper.

Appendix A Proofs of the theorems

Proof A.1** (Lemma 2.1).**

This bound on the Lagrange basis is a classical result related to the growth of the Lebesgue constant in polynomial interpolation. For ${\overline{m}}$ Chebyshev nodes of the first kind on $[-1,1]$ and the associated Lagrange basis functions $\ell_{1},\dots,\ell_{\overline{m}}$ we have the following result [23, equation 13]

[TABLE]

This implies that in one dimension,

[TABLE]

Going from one to $d$ dimensions can be done using Kronecker products. Indeed, for $x\in\mathbb{R}^{d}$,

[TABLE]

where $x=(x_{1},\dots,x_{d})$ and $\overline{X}_{1},\dots,\overline{X}_{d}$ are the one-dimensional Chebyshev nodes. Since for all $a\in\mathbb{R}^{m},b\in\mathbb{R}^{n}$, $\|a\otimes b\|_{2}=\sqrt{\sum_{i,j}(a_{i}b_{j})^{2}}=\|ab^{\top}\|_{F}=\|a\|_{2}\|b\|_{2}$, it follows that

[TABLE]

This implies, using a fairly loose bound,

[TABLE]

The same argument can be done for $T(y,\overline{Y})$.

In 1D, the weights are

[TABLE]

for $k=1,\dots,{\overline{m}}$. Obviously, $w_{k}>0$. Clearly, $w_{k}<\frac{\pi}{{\overline{m}}}$. Also, the minimum being reached at $k=1$ or $k={\overline{m}}$,

[TABLE]

Since the nodes in $d$ dimensions are products of the nodes in 1D, it follows that

[TABLE]

*The result follows. *

Proof A.2** (Lemma 2.4).**

We show the result for the second equation. This requires using, consecutively, the interpolation result and the CUR decomposition one. First, one can write from Lemma 2.1 and the interpolation,

[TABLE]

Then, introducing the weight matrices and applying Lemma 2.2 on the interpolation matrix,

[TABLE]

Finally, combining and distributing all the factors gives us

[TABLE]

Here, we can bound all terms:

• •

For the first term, Lemma 2.1 and Lemma 2.2 show that the expression is bounded by a polynomial;

• •

For the second term use the fact that

[TABLE]

hence, since $\|E_{QR}(\overline{X},\overline{Y})\|_{2}=\epsilon$, the product is again bounded by a polynomial since the $\epsilon$ cancel out ;

• •

The last term can be bounded in a similar way using

[TABLE]

We conclude that there exists a polynomial $q$ such that

[TABLE]

*The proof is similar in $x$. *

Proof A.3** (Theorem 2.5).**

Combining interpolation and CUR decomposition results one can write

[TABLE]

Distributing everything, factoring the weights matrices and simplifying, we obtain the following, where we indicate the bounds on each term on the right,

[TABLE]

*This concludes the proof. *