# Spectral properties of kernel matrices in the flat limit

Simon Barthelme, Konstantin Usevich

► To cite this version:

Simon Barthelme, Konstantin Usevich. Spectral properties of kernel matrices in the flat limit. SIAM Journal on Matrix Analysis and Applications, 2021, 42 (1), pp.17-57. 10.1137/19M129677X. hal-03017947v2

HAL Id: hal-03017947

https://hal.science/hal-03017947v2

Submitted on 27 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# SPECTRAL PROPERTIES OF KERNEL MATRICES
## IN THE FLAT LIMIT*

SIMON BARTHELMÉ† AND KONSTANTIN USEVICH‡

**Abstract.** Kernel matrices are of central importance to many applied fields. In this manuscript, we focus on spectral properties of kernel matrices in the so-called “flat limit”, which occurs when points are close together relative to the scale of the kernel. We establish asymptotic expressions for the determinants of the kernel matrices, which we then leverage to obtain asymptotic expressions for the main terms of the eigenvalues. Analyticity of the eigenprojectors yields expressions for limiting eigenvectors, which are strongly tied to discrete orthogonal polynomials. Both smooth and finitely smooth kernels are covered, with stronger results available in the finite smoothness case.

**Key words.** kernel matrices, eigenvalues, eigenvectors, radial basis functions, perturbation theory, flat limit, discrete orthogonal polynomials

**AMS subject classifications.** 15A18, 47A55, 47A75, 47B34, 60G15, 65D05

**1. Introduction.** For an ordered set of points $\mathcal{X} = (x_1, \dots, x_n)$, $x_k \in \Omega \subset \mathbb{R}^d$, not lying in general on a regular grid, and a kernel $K: \Omega \times \Omega \to \mathbb{R}$, the kernel matrix $K$ is defined as

$$K = K(\mathcal{X}) = [K(x_i, x_j)]_{i,j=1}^{n,n}.$$

These matrices occur in approximation theory (kernel-based approximation and interpolation, [33, 43]), statistics and machine learning (Gaussian process models [44], Support Vector Machines and kernel PCA [34]). Often, a positive scaling parameter is introduced<fnref n="1" />, and the scaled kernel matrix becomes

$$(1) \qquad K_{\varepsilon} = K_{\varepsilon}(\mathcal{X}) = [K_{\varepsilon}(x_i, x_j)]_{i,j=1}^{n,n},$$

where typically $K_{\varepsilon}(x, y) = K(\varepsilon x, \varepsilon y)$. If the kernel is radial (the most common case), then its value depends only on the Euclidean distance between $x$ and $y$, and $\varepsilon$ determines how quickly the kernel decays with distance.

Understanding spectral properties of kernel matrices is essential in statistical applications (e.g., for selecting hyperparameters), as well as in scientific computing (e.g., for preconditioning [13, 41]). Because the spectral properties of kernel matrices are not directly tractable in the general case, one usually needs to resort to asymptotic results. The most common form of asymptotic analysis takes $n \to \infty$. Three cases are typically considered: (a) when the distribution of points in $\mathcal{X}$ converges to some

<footnote>*Submitted to the editors on DATE.</footnote>

<footnote>**Funding:** This work was supported by ANR project GenGP (ANR-16-CE23-0008) and ANR project LeaFleT (ANR-19-CE23-0021-01).</footnote>

<footnote>†CNRS, Univ. Grenoble Alpes, Grenoble INP, GIPSA-lab, 38000 Grenoble, France (simon.barthelme@gipsa-lab.fr).</footnote>

<footnote>‡Université de Lorraine and CNRS, CRAN (Centre de Recherche en Automatique en Nancy), UMR 7039, Campus Sciences, BP 70239, 54506 Vandœuvre-lès-Nancy cedex, France (konstantin.usevich@univ-lorraine.fr).</footnote>

<footnote n="1">¹In this paper, for simplicity, we consider only the case of isotropic scaling (i.e., all the variables are scaled with the same parameter $\varepsilon$). The results should hold for the constant anisotropic case, by rescaling the set of points $\mathcal{X}$ in advance.</footnote>
continuous measure on $\Omega$, the kernel matrix tends in some sense to a linear operator in a Hilbert space, whose spectrum is then studied [39]; (b) recently, some authors have obtained asymptotic results in a regime where both $n$ and the dimension $d$ tend to infinity [11, 7, 40], using the tools of random matrix theory; (c) in a special case of $\mathcal{X}$ lying on a regular grid, stationary kernel matrices become Toeplitz or more generally multilevel Toeplitz, whose asymptotic spectral distribution is determined by their symbol (or the Fourier transform of the sequence) [18, 2, 38, 29, 3].

Driscoll & Fornberg [10] pioneered a new form of asymptotic analysis for kernel methods, in the context of Radial Basis Function (RBF) interpolation. The point set $\mathcal{X}$ is considered fixed, with arbitrary geometry (i.e., not lying in general on a regular grid), and the scaling parameter $\varepsilon$ approaches 0. Driscoll & Fornberg called this the “flat limit”, as kernel functions become flat over the range of $\mathcal{X}$ as $\varepsilon \to 0$. Very surprisingly, they showed that for certain kernels the RBF interpolant stays well-defined in the flat limit, and tends to the Lagrange polynomial interpolant. Later, a series of papers extended their results to the multivariate case [31, 24, 32] and established similar convergence results for various types of radial functions (see [35, 27] and references therein). In particular, [35] showed that for kernels of finite smoothness the limiting interpolant is a spline rather than a polynomial.

The flat limit is interesting for several reasons. In contrast to other asymptotic analyses, it is deterministic ($\mathcal{X}$ is fixed), and makes very few assumptions on the geometry of the point set. In addition, kernel methods are plagued by the problem of picking a scale parameter [34]. One either uses burdensome procedures like cross-validation or maximum likelihood [44] or sub-optimal but cheap heuristics like the median distance heuristic [16]. The flat limit analysis may shed some light on the problem. Finally, the results derived here can be thought of as perturbation results, in the sense that they are formally exact in the limit, but useful approximations when the scale is not too small.

Despite its importance, little was known until recently about the eigenstructure of kernel matrices in the flat limit. The difficulty comes from the fact that $K_\varepsilon = K(0,0)\mathbf{1}\mathbf{1}^\mathrm{T} + \mathcal{O}(\varepsilon)$, where $\mathbf{1} = [1 \cdots 1]^\mathrm{T}$, i.e., we are dealing with a singular perturbation problem$^2$.

Schaback [31, Theorem 6], and, more explicitly, Wathen & Zhu [42] obtained results on the orders of eigenvalues of kernel matrices for smooth analytic radial basis kernels, based on the Courant-Fischer minimax principle. A heuristic analysis of the behaviour of the eigenvalues in the flat limit was also performed in [15]. However, the main terms in the expansion of the eigenvalues have never been obtained, and the results in [31, 42] apply only to smooth kernels. In addition, they hold no direct information on the limiting eigenvectors.

In this paper, we try filling this gap by characterising both the eigenvalues and eigenvectors of kernel matrices in the flat limit. We consider both completely smooth kernels and finitely-smooth kernels. The latter (Matérn-type kernels) are very popular in spatial statistics. For establishing asymptotic properties of eigenvalues, we use the expression for the limiting determinants of $K_\varepsilon$ (obtained only for the smooth case), and Binet-Cauchy formulae. As a special case, we recover the results of Schaback [31] and Wathen & Zhu [42], but for a wider class of kernels.

**1.1. Overview of the results.** Some of the results are quite technical, so the goal of this section is to serve as a reader-friendly summary of the contents.

<footnote>$^2$Seen from the point of view of the characteristic polynomial, the equation $\det(K_\varepsilon - \lambda I) = 0$ has a solution of multiplicity $n-1$ at $\varepsilon = 0$, but these roots immediately separate when $\varepsilon > 0$.</footnote>
**1.1.1. Types of kernels.** We begin with some definitions.

A kernel is called translation invariant if

$$ (2) \qquad K(\mathbf{x}, \mathbf{y}) = \varphi(\mathbf{x} - \mathbf{y}) $$

for some function $\varphi$. A kernel is radial, or stationary if, in addition, we have:

$$ (3) \qquad K(\mathbf{x}, \mathbf{y}) = f(\|\mathbf{x} - \mathbf{y}\|_2) $$

i.e., $\varphi(t) = f(|t|)$ and the value of the kernel depends only on the Euclidean distance between $\mathbf{x}$ and $\mathbf{y}$. The function $f$ is an RBF. Finally, $K(\mathbf{x}, \mathbf{y})$ may be positive (semi) definite, in which case the kernel matrix is positive (semi) definite [9] for all sets of distinct points $\mathcal{X} = (\mathbf{x}_1, \dots, \mathbf{x}_n)$, and all $n > 0$.

All of our results are valid for stationary, positive semi-definite kernels. In addition, some are also valid for translation-invariant kernels, or even general, non-radial kernels. For simplicity, we focus on radial kernels in this introductory section.

An important property of a radial kernel is its order of smoothness, which we call $r$ throughout this paper. The definition is at first glance not very enlightening: formally, if the first $p$ odd-order derivatives of the RBF $g$ are zero, and the $(p+1)$-th is nonzero, then $r = p + 1$. To understand the definition, some Fourier analysis is required [37], but for the purposes of this article we will just note two consequences. When interpolating using a radial kernel of smoothness $r$, the resulting interpolant is $r - 1$ times differentiable. When sampling from a Gaussian process with covariance function of smoothness order $r$, the sampled process is also $r - 1$ times differentiable (almost surely). $r$ may equal $\infty$, which is the case we call infinitely smooth. If $r$ is finite we talk about a finitely smooth kernel. We treat the two cases separately due to importance of infinitely smooth kernels, and because proofs are simpler in that case.

Finally, the points are assumed to lie in some subset of $\mathbb{R}^d$, and if $d = 1$ we call this the univariate case, as opposed to the multivariate case ($d > 1$).

**1.1.2. Univariate results.** In the univariate case, we can give simple closed-form expressions for the eigenvalues and eigenvectors of kernel matrices as $\varepsilon \to 0$. What form these expressions take depends essentially on the order of smoothness.

We shall contrast two kernels that are at opposite ends of the smoothness spectrum. One, the Gaussian kernel, is infinitely smooth, and is defined as:

$$ K_{\varepsilon}(x, y) = \exp(-\varepsilon(x - y)^2). $$

The other has smoothness order 1, and is known as the “exponential” kernel (and is also a Matérn kernel):

$$ K_{\varepsilon}(x, y) = \exp(-\varepsilon|x - y|). $$

Both kernels are radial and positive definite. However, the small-$\varepsilon$ asymptotics of these two kernels are strikingly different.

In the case of the Gaussian kernel, the eigenvalues go to 0 extremely fast, except for the first one, which goes to $n$. Specifically, the first eigenvalue is $O(1)$, the second is $O(\varepsilon^2)$, the third is $O(\varepsilon^4)$, etc$^3$. Fig. 1.1 shows the eigenvalues of the Gaussian kernel for a fixed set $\mathcal{X}$ of randomly-chosen nodes in the unit interval ($n = 10$ here). The eigenvalues are shown as a function of $\varepsilon$, under log-log scaling. As expected from Theorem 4.2 (see also [31, 42]), for each $i$, $\log \lambda_i$ is approximately linear as a function

<footnote>$^3$In fact, Theorem 4.2 guarantees the same behaviour for general kernels of sufficient smoothness.</footnote>
<figure><img src="image_2.png" /><figcaption>FIG. 1.1. Eigenvalues of the Gaussian kernel ($d = 1$). The set of $n = 10$ nodes was drawn uniformly from the unit interval. In black, eigenvalues of the Gaussian kernel, for different values of $\epsilon$. The dashed red curves are our small-$\epsilon$ expansions. Note that both axes are scaled logarithmically. The noise apparent for small $\epsilon$ values in the low range is due to loss of precision in the numerical computations.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>FIG. 1.2. Eigenvalues of the exponential kernel ($d = 1$), using the same set of points as in Fig. 1.1. The largest eigenvalue has a slope of 0 for small $\epsilon$, the others have unit slope, as in Theorem 4.5.</figcaption></figure>

of $\log \epsilon$. In addition, the main term in the scaling of $\log \lambda_i$ in $\epsilon$ (i.e., the offsets of the various lines) is also given by Theorem 4.2, and the corresponding asymptotic approximations are plotted in red. They show very good agreement with the exact eigenvalues, up to $\epsilon \approx 1$.

Contrast that behaviour with the one exhibited by the eigenvalues of the exponential kernel. Theorem 4.6 describes the expected behaviour: the top eigenvalue is again $O(1)$ and goes to $n$, while all remaining eigenvalues are $O(\epsilon)$. Fig. 1.2 is the counterpart of the previous figure, and shows clearly that all eigenvalues except for the top one go to 0 at unit rate. The main term in the expansions of eigenvalues determines again the offsets shown in Fig. 1.2, which can be computed from the eigenvalues of the centered distance matrix as shown in Theorem 4.6.
<figure><img src="image_4.png" /><figcaption>FIG. 1.3. First four eigenvectors of the Gaussian kernel. We used the same set of points as in the previous two figures. In blue, eigenvectors of the Gaussian kernel, for different values of ε. The dashed red curve shows the theoretical limit as ε → 0 (i.e., the first four orthogonal polynomials of the discrete measure)</figcaption></figure>

To sum up: except for the top eigenvalue, which behaves in the same way for both kernels, the rest scale quite differently. More generally, Theorem 4.6 states that for kernels of smoothness order $r < n$ ($r = 1$ for the exponential, $r = \infty$ for the Gaussian), the eigenvalues are divided into two groups. The first group of eigenvalues is of size $r$, and have orders 1, $\epsilon^2$, $\epsilon^4$, etc. The second group is of size $n-r$, and all have the same order, $\epsilon^{2r-1}$.

The difference between the two kernels is also reflected in how the eigenvectors behave. For the Gaussian kernel, the limiting eigenvectors (shown in Fig. 1.3) are columns of the *Q* matrix of the QR factorization of the Vandermonde matrix (i.e., the orthogonal polynomials with respect to the discrete uniform measure on *X*). For instance, the top eigenvector equals the constant vector $\frac{1}{\sqrt{n}}\mathbf{1}$, and the second eigenvector equals $[x_1 \cdots x_n]^T - \frac{1}{n}\sum_{i=1}^n x_i$ (up to normalisation). Each successive eigenvector depends on the geometry of *X* via the moments $m_p(X) = \sum_{i=1}^n x_i^p$. In fact, this result is valid for any positive definite smooth analytic in $\epsilon$ kernel as shown by Corollary 4.3.

In the case of finite smoothness, the two groups are associated with different groups of eigenvectors. The first group of $r$ eigenvectors are again orthogonal polynomials. The second group are splines of order $2r-1$. Convergence of eigenvectors is shown in Fig. 1.4. This general result for eigenvectors is shown in Theorem 8.1.

#### 1.1.3. The multivariate case.
The general, multivariate case requires more care. Polynomials continue to play a central role in the flat limit, but when $d > 1$ they appear naturally in groups of equal degree. For instance, in $d = 2$, we may write $\mathbf{x} = [yz]^T$ and the first few monomials are as follows:

* Degree 0: 1

* Degree 1: $y, z$

* Degree 2: $y^2, yz, z^2$

* Degree 3: $y^3, y^2z, yz^2, z^3$

etc. Note that there is one monomial of degree 0, two degree-1 monomials, three degree-2 monomials, and so on. If $d = 1$ there is a single monomial in each group,
<figure><img src="image_5.png" /><figcaption>FIG. 1.4. First four eigenvectors of the exponential kernel. The same set of nodes is used as in Fig. 1.3. In blue, the eigenvectors of the kernel, for different values of ε. The dashed red curve shows the theoretical limit as ε → 0. From Theorem 8.1, these are (1) the vector $\frac{1}{\sqrt{n}}$ and (2-4) the first three eigenvectors of $(I - \frac{11}{n}T) D_{(1)} (I - \frac{11}{n}T)$, where $D_{(1)}$ is the Euclidean distance matrix.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIG. 1.5. Eigenvalues of the Gaussian kernel in the multivariate case ($d = 2$). The set of $n = 10$ nodes was drawn uniformly from the unit square. In black, eigenvalues of the Gaussian kernel, for different values of $\epsilon$. The dashed red curves are our small-$\epsilon$ expansions. Eigenvalues appear in groups of different orders, recognizable here as having the same slopes as $\epsilon \to 0$. A single eigenvalue is of order $O(1)$. Two eigenvalues are of order $O(\epsilon^2)$, as many as there are monomials of degree 1 in two dimensions. Three eigenvalues are of order $O(\epsilon^4)$, as many as there are monomials of degree 2 in two dimensions, etc.</figcaption></figure>

and here lies the essence of the difference between the univariate and multivariate cases.

In infinitely smooth kernels like the Gaussian kernel, as shown in [31, 42], there are as many eigenvalues of order $O(\epsilon^{2k})$ as there are monomials of degree $k$ in dimension $d$: for instance, there are 4 monomials of degree 3 in dimension 2, and so 4 eigenvalues of order $O(\epsilon^8)$. An example is shown in Fig. 1.5. In finitely smooth kernels like the exponential kernel, there are $r$ successive groups of eigenvalues with order $O(1)$, $O(\epsilon^2)$,
<figure><img src="image_7.png" /><figcaption>FIG. 1.6. Eigenvalues of a finitely smooth kernel in $d = 2$. Here we picked a kernel with smoothness index $r = 3$ (the exponential kernel has order $r = 1$). The set of $n = 20$ nodes was drawn uniformly from the unit square. The first $r = 3$ groups of eigenvalues have the same behaviour as in the Gaussian kernel: one is $O(1)$, two are $O(\varepsilon^2)$, three are $O(\varepsilon^4)$. All the rest are $O(\varepsilon^5)$.</figcaption></figure>

..., up to $O(\varepsilon^{2r-2})$. Following that, all remaining eigenvalues have order $O(\varepsilon^{2r-1})$, just like in the one-dimensional case.

We show in Theorem 6.2, that the main terms of these eigenvalues in each group can be computed from the QR factorization of the Vandermonde matrix and a Schur complement associated with the Wronskian matrix of the kernel. In the finite smoothness case (see Fig. 1.6), the same expansions are valid until the smoothness order, and the last group of eigenvalues is given by the eigenvalues of the projected distance matrix, as in the one-dimensional case.

Finally, in the multivariate case, Theorem 6.2 characterises the eigenprojectors. In a nutshell, the invariant subspaces associated with each group of eigenvalues are spanned by orthogonal polynomials of a certain order. The eigenvectors are the subject of a conjecture given in section 8, which we believe to be quite solid.

**1.2. Overview of tools used in the paper.** For finding the orders of the eigenvalues, as in [31, 42], we use the Courant-Fischer minimax principle (see Theorem 3.7). However, unlike [31, 42], we do not use directly the results of Micchelli [28], but rather rotate the kernel matrices using the $Q$ factor in the QR factorization of the Vandermonde matrix, and find the expansion of rotated matrices from the Taylor expansion of the kernel.

The key results are the expansions for the determinants of $K_\varepsilon$, which use the expansions of rotated kernel matrices. Our results on determinants (Theorems 4.1, 4.5, 6.1, and 6.3) generalize those of Lee & Micchelli [26]. The next key observation is that principal submatrices of $K_\varepsilon$ are also kernel matrices, hence the results on determinants imply the results on expansions of elementary symmetric polynomials of eigenvalues (via the correspondence between elementary symmetric polynomials, see Theorem 3.1, and the Binet-Cauchy formula). Finally, the main terms of the eigenvalues can be retrieved from the main terms of the elementary symmetric polynomials, as shown in Lemma 3.3. An important tool for the multivariate and finite smoothness case is Lemma 3.12 on low-rank perturbation of elementary symmetric polynomials that we could not find elsewhere in the literature
To study the properties of the eigenvectors, we use analyticity of the eigenvalues and eigenprojectors for Hermitian analytic matrix-valued functions [22]. By using an extension of the Courant-Fischer principle (see Lemma 3.10), we can fully charac- terise the limiting eigenvectors in the univariate case, and obtain the limiting invariant subspaces for the groups of the limiting eigenvectors in the multivariate case. More- over, by using the perturbation expansions from [22], we can find the last individual eigenvectors in the finitely smooth case.

We note that the multivariate case requires a number of technical assumptions on the arrangement of points, which are typical for multivariate polynomial interpolation. For the results on determinants, no assumptions are required. However, for getting the results on the eigenvalues or eigenvectors at a certain order of ε, we need a relaxed version to the well-known unisolvency condition [31, 14], namely the multivariate Vandermonde matrix up to the corresponding degree to be full rank.

**1.3. Organisation of the paper.** In an attempt to make the paper reader-friendly, it is organized as follows. In Section 2 we recall the main terminology for 1D kernels. In Section 3 we gather well-known (or not so well-known) results on eigenvalues, determinants, elementary symmetric polynomials and their perturbations, which are key tools used in the paper. Section 4 contains the main results on determinants, eigenvalues, and eigenvectors in the univariate ($d = 1$) case. While these results are special cases of the multivariate case ($d > 1$), the latter is burdened with heavier notation due to the complexity of dealing with multivariate polynomials. To get a gist of the results and techniques, the reader is advised to first consult the case $d = 1$. In Section 5, we introduce all the needed notation to handle the multivariate case. Section 6 contains the main results of the paper on determinants, eigenvalues, and groups of eigenvectors in the multivariate case. Thus Sections 2 to 6 are quite self-contained and contain most of the results in the paper, except the result on precise locations of the last group of eigenvectors in the finite smoothness case. In Section 7, we provide a brief summary on analytic perturbation theory, needed only for proving the stronger result on eigenvectors for finitely smooth kernels in Section 8.

**2. Background and main notation.** This section contains main definitions and examples of kernels, mostly in the 1D case (with multivariate extensions given in Section 5). We assume that the kernel $K: \Omega \times \Omega \to \mathbb{R}$ is in the class $C^{(\ell,\ell)}(\Omega)$, $\Omega = (-a; a)$, i.e. all the partial derivatives $\frac{\partial^2}{\partial x^i \partial y^j} K$ exist and are continuous for $0 \le i, j \le \ell$ on $\Omega \times \Omega$.

We will often need the following short-hand notation for partial derivatives

$$K^{(i,j)} \stackrel{\text{def}}{=} \frac{\partial^{i+j}}{\partial x^i \partial y^j} K,$$

and we introduce the so-called Wronskian matrix

$$ (4) \qquad W_{\le k} \stackrel{\text{def}}{=} \begin{bmatrix} \frac{K^{(0,0)}(0,0)}{0!0!} & \cdots & \frac{K^{(0,k)}(0,0)}{0!k!} \\ \vdots & \ddots & \vdots \\ \frac{K^{(k,0)}(0,0)}{k!0!} & \cdots & \frac{K^{(k,k)}(0,0)}{k!k!} \end{bmatrix}. $$

**2.1. Translation invariant and radial kernels.** Let us consider an important example of a translation invariant (2) kernel, which in the univariate case becomes $K(x, y) = \varphi(x - y)$. We assume that $\varphi \in C^{2r}(-2a, 2a)$; hence, $\varphi$ has a Taylor
expansion around 0

$$ \varphi(t) = \alpha_0 + \alpha_1 t + \alpha_2 t^2 + \dots + \alpha_{2r} t^{2r} + o(t^{2r}), \quad \text{where } \alpha_k \stackrel{\text{def}}{=} \frac{\varphi^{(k)}(0)}{k!}. $$

Therefore $K \in C^{(r,r)}(\Omega \times \Omega)$ and its derivatives at 0 are

$$ (5) \qquad \frac{K^{(i,j)}(0,0)}{i!j!} = \frac{(-1)^j \varphi^{(i+j)}(0)}{i!j!} = (-1)^j \binom{i+j}{j} \alpha_{i+j}, $$

i.e., the Wronskian matrix has the form:

$$ W_{\le k} = \begin{bmatrix} \alpha_0 & -\alpha_1 & \alpha_2 & -\alpha_3 & \cdots \\ \alpha_1 & -2\alpha_2 & 3\alpha_3 & -4\alpha_4 & \cdots \\ \alpha_2 & -3\alpha_3 & 6\alpha_4 & -10\alpha_5 & \cdots \\ \alpha_3 & -4\alpha_4 & 10\alpha_5 & -20\alpha_6 & \cdots \\ \vdots & \vdots & \vdots & \vdots & \ddots \end{bmatrix}. $$

A special case of the translational kernels are smooth and radial, which will be considered in the next subsection. The simplest example is the Gaussian kernel with

$$ \varphi(t) = \exp(-t^2) = 1 - t^2 + \frac{t^4}{2!} - \frac{t^6}{3!} + \dots, $$

i.e. for all integer $\nu \ge 0$

$$ \alpha_{2\nu} = \frac{(-1)^\nu}{\nu!}, \quad \alpha_{2\nu+1} = 0. $$

In this case, the Wronskian matrix becomes

$$ W_{\le k} = \begin{bmatrix} 1 & 0 & -1 & 0 & \cdots \\ 0 & 2 & 0 & -2 & \cdots \\ -1 & 0 & 3 & 0 & \cdots \\ 0 & -2 & 0 & -\frac{10}{3} & \cdots \\ \vdots & \vdots & \vdots & \vdots & \ddots \end{bmatrix}. $$

An important subclass consists of radial kernels (3), i.e., $K(x, y) = f(|x - y|)$ in the univariate case. We put the following assumptions on $f$:

* $f \in C^{2r}(-2a, 2a)$;

* the highest derivative of odd order is not zero, i.e., $f^{(2r-1)}(0) \neq 0$;

* the lower derivatives with odd order vanish, i.e., $f^{(2\ell-1)}(0) = 0$ for $\ell < r$.

If the function $f(t)$ satisfies these assumptions, then $r$ is called the order of smoothness of $K$. Note that $f$ admits a Taylor expansion

$$ (6) \qquad f(t) = f_0 + f_2 t^2 + \dots + f_{2r-2} t^{2r-2} + t^{2r-1} (f_{2r-1} + \mathcal{O}(t)), $$

where $f_k = f^{(k)}(0)/k!$ is a shorthand notation for the scaled derivative at 0. For example, for the exponential kernel $\exp(-|x - y|)$, we have $f_0 = 1$ and $f_1 = -1$, so the smoothness order is $r = 1$. For the $C^2$ Matérn kernel $(1 + |x - y|) \exp(-|x - y|)$, we have $f_0 = 1, f_1 = 0, f_2 = -1/2, f_3 = 1/3$, so the smoothness order is $r = 2$.
## 2.2. Distance matrices, Vandermonde matrices, and their properties.

In the general multivariate case, we define $D_{(k)}$ to be the k-th Hadamard power of the Euclidean distance matrix $D_{(1)}$, i.e.,

$$D_{(k)} = \left[ \|\mathbf{x}_i - \mathbf{x}_j\|_2^k \right]_{i,j=1}^{n,n}.$$

Next, we focus on the univariate case (the multivariate case will be discussed later in Section 5) and denote by $V_{\le k}$ the univariate Vandermonde matrix up to degree $k$

$$ (7) \qquad V_{\le k} = \begin{bmatrix} 1 & x_1 & \cdots & x_1^k \\ \vdots & \vdots & & \vdots \\ 1 & x_n & \cdots & x_n^k \end{bmatrix}, $$

which has rank $\min(n, k+1)$ if the nodes are distinct. In particular, the matrix $V_{\le n-1}$ is square and invertible for distinct nodes.

For even $k$, the Hadamard powers $D_{(k)}$ of the distance matrix can be expressed via the columns of the Vandermonde matrix using the binomial expansion

$$ (8) \qquad D_{(2\ell)} = \sum_{j=0}^{2\ell} (-1)^j \binom{2\ell}{j} v_j v_{2\ell-j}^{\top}, $$

where $v_k \stackrel{\text{def}}{=} [x_1^k \cdots x_n^k]^{\top}$ are columns of Vandermonde matrices. Therefore, rank $D_{(2\ell)} = \min(2\ell + 1, n)$ if all points $x_k$ are distinct.

On the other hand, for $k$ odd, the matrices $D_{(k)}$ exhibit an entirely different set of properties. Of most interest here is conditional positive-definitess, which guarantees that the distance matrices are positive definite when projected onto a certain subspace. The following result appears e.g., in [12, Chapter 8] but follows directly from an earlier paper by Micchelli [28].

**LEMMA 2.1.** For a distinct node set $\mathcal{X} \subset \mathbb{R}$ of size $n > r \ge 1$, we let $\mathbf{B}$ be a full column rank matrix such that $\mathbf{B}^{\top}V_{\le r-1} = 0$. Then the matrix $(-1)^r \mathbf{B}^{\top}D_{(2r-1)}\mathbf{B}$ is positive definite.

For instance, if $r = 1$, we may pick any basis $\mathbf{B}$ orthogonal to the vector $\mathbf{1}$, and the lemma implies that $(-\mathbf{B}^{\top}D_{(1)}\mathbf{B})$ has $n-1$ positive eigenvalues. We note for future reference that the result generalises almost unchanged to the multivariate case.

## 2.3. Scaling and expansions of kernel matrices.

In this subsection, we consider the general multivariate case. Given a general kernel $K$, we define its scaled version as

$$ (9) \qquad K_{\varepsilon}(\mathbf{x}, \mathbf{y}) = K(\varepsilon\mathbf{x}, \varepsilon\mathbf{y}), $$

while for the specific case of radial kernels we use the following form

$$ (10) \qquad K_{\varepsilon}(\mathbf{x}, \mathbf{y}) = f(\varepsilon \|\mathbf{x} - \mathbf{y}\|_2). $$

Note that the definitions (9) and (10) coincide for nonnegative $\varepsilon$, but differ if we formally take other values (e.g., complex) of $\varepsilon$. Depending on context, we will use one or the other of the definitions later on, especially when we talk about analyticity of kernel matrices.
Using (6) and (10), we may write the scaled kernel matrix (1), as

$$ (11) \quad \mathbf{K}_{\varepsilon} = f_0 \mathbf{D}_{(0)} + \varepsilon^2 f_2 \mathbf{D}_{(2)} + \dots + \varepsilon^{2r-2} f_{2r-2} \mathbf{D}_{(2r-2)} + \varepsilon^{2r-1} f_{2r-1} \mathbf{D}_{(2r-1)} + \mathcal{O}(\varepsilon^{2r}). $$

In the univariate case, (8) gives a way to rewrite the expansion (11) as

$$ (12) \quad \mathbf{K}_{\varepsilon} = \sum_{\ell=0}^{r-1} \varepsilon^{2j} \underbrace{\mathbf{V}_{\le 2\ell} \mathbf{W}_{/2\ell} \mathbf{V}_{\le 2\ell}^{\top}}_{f_{2\ell} \mathbf{D}_{(2\ell)}} + \varepsilon^{2r-1} f_{2r-1} \mathbf{D}_{(2r-1)} + \mathcal{O}(\varepsilon^{2r}), $$

where $\mathbf{W}_s \in \mathbb{R}^{(s+1) \times (s+1)}$ has nonzero elements only on its antidiagonal:

$$ (13) \quad \left( \mathbf{W}_s \right)_{i+1, s-i+1} \stackrel{\text{def}}{=} f_s(-1)^i \binom{s}{i}. $$

For example,

$$ \mathbf{W}_{/2} = f_2 \begin{bmatrix} 0 & 0 & 1 \\ 0 & -2 & 0 \\ 1 & 0 & 0 \end{bmatrix}. $$

In fact, from (5), non-zero elements of $\mathbf{W}_s$ are scaled derivatives of the kernel:

$$ \left( \mathbf{W}_s \right)_{i+1, s-i+1} = \frac{K^{(i,s-i)}(0,0)}{i!(s-i)!}; $$

this justifies the notation $\mathbf{W}_s$ (i.e., an antidiagonal of the Wronskian matrix).

**3. Determinants and elementary symmetric polynomials.** In this paper, we will heavily use the elementary symmetric polynomials of eigenvalues, and we collect in this section some useful (more or less known) facts about them.

**3.1. Eigenvalues, principal minors and elementary symmetric polynomials.** The k-th elementary symmetric polynomial of $\lambda_1, \dots, \lambda_n$ is defined as:

$$ (14) \quad e_k(\lambda_1, \dots, \lambda_n) = \sum_{\substack{\mathcal{Y}=k \\ \mathcal{Y} \subset \{1, \dots, n\}}} \prod_{i \in \mathcal{Y}} \lambda_i, $$

i.e., the sum is running over all possible subsets of $\{1, \dots, n\}$ of size $k$. In particular, $e_1(\lambda_1, \dots, \lambda_n) = \sum_{i=1}^n \lambda_i$, $e_2(\lambda_1, \dots, \lambda_n) = \sum_{i<j} \lambda_i \lambda_j$, and $e_n(\lambda_1, \dots, \lambda_n) = \prod_{i=1}^n \lambda_i$.

Next, we consider the elementary symmetric polynomials applied to eigenvalues of matrices, and define (with some abuse of notation):

$$ e_k(\mathbf{A}) \stackrel{\text{def}}{=} e_k(\lambda_1(\mathbf{A}), \dots, \lambda_n(\mathbf{A})). $$

Then the first and the last polynomials are the trace and determinant of $\mathbf{A}$:

$$ e_1(\mathbf{A}) = \operatorname{Tr} \mathbf{A}, \quad e_n(\mathbf{A}) = \det \mathbf{A}. $$

This fact is a special case of a more general result on sums of principal minors.

**THEOREM 3.1 ([19, Theorem 1.2.12]).**

$$ (15) \quad e_k(\mathbf{A}) = \sum_{\substack{\mathcal{Y} \subset \{1, \dots, n\} \\ \#\mathcal{Y}=k}} \det(\mathbf{A}_{\mathcal{Y}, \mathcal{Y}}) $$

where $\mathbf{A}_{\mathcal{Y}, \mathcal{Y}}$ is a submatrix of $\mathbf{A}$ with rows and columns indexed by $\mathcal{Y}$, i.e. the sum runs over all principal minors of size $k \times k$.

**REMARK 3.2.** The scaled symmetric polynomials $(-1)^k e_k(\mathbf{A})$ are the coefficients of the characteristic polynomial $\det(\lambda \mathbf{I} - \mathbf{A})$ at the coefficient $\lambda^{n-k}$.
**3.2. Orders of elementary symmetric polynomials.** Next, we assume that $\lambda_1(\varepsilon), \dots, \lambda_n(\varepsilon)$ are functions of some small parameter $\varepsilon$ and we are interested in the orders of the corresponding elementary symmetric polynomials

$$e_s(\varepsilon) \stackrel{\text{def}}{=} e_s(\lambda_1(\varepsilon), \dots, \lambda_n(\varepsilon)), \quad 1 \le s \le n,$$

as $\varepsilon \to 0$. The following obvious observation will be important.

**LEMMA 3.3.** *Assume that*

$$\lambda_1(\varepsilon) = \mathcal{O}(\varepsilon^{L_1}), \quad \lambda_2(\varepsilon) = \mathcal{O}(\varepsilon^{L_2}), \dots, \quad \lambda_n(\varepsilon) = \mathcal{O}(\varepsilon^{L_n}),$$

as $\varepsilon \to 0$ and $0 \le L_1 \le \dots \le L_n$ are some integers. Then it holds that

$$e_1(\varepsilon) = \mathcal{O}(\varepsilon^{L_1}), \quad e_2(\varepsilon) = \mathcal{O}(\varepsilon^{L_1+L_2}), \dots, \quad e_n(\varepsilon) = \mathcal{O}(\varepsilon^{L_1+\dots+L_n}).$$

*Proof.* The proof follows from the definition of $e_s$, the fact that the product $f(\varepsilon) = \mathcal{O}(\varepsilon^a)$ and $g(\varepsilon) = \mathcal{O}(\varepsilon^b)$ is of the order $f(\varepsilon)g(\varepsilon) = \mathcal{O}(\varepsilon^{a+b})$. $\square$

We will need a refinement of Lemma 3.3 concerning the main terms of such functions. We distinguish two situations: when the orders of $\lambda_k$ are separated, and when they form groups. For example, if $\lambda_1(\varepsilon) = a + \mathcal{O}(\varepsilon)$, $\lambda_2(\varepsilon) = \varepsilon(b + \mathcal{O}(\varepsilon))$, $\lambda_3(\varepsilon) = \varepsilon^2(c + \mathcal{O}(\varepsilon))$, then the main terms of elementary symmetric polynomials are products of the main terms of $\lambda_k$:

$$e_1(\varepsilon) = a + \mathcal{O}(\varepsilon), \quad e_2(\varepsilon) = \varepsilon(ab + \mathcal{O}(\varepsilon)), \quad e_3(\varepsilon) = \varepsilon^3(abc + \mathcal{O}(\varepsilon)).$$

On the other hand, in the case $\lambda_1(\varepsilon) = a + \mathcal{O}(\varepsilon)$, $\lambda_{k+1}(\varepsilon) = \varepsilon(b_k + \mathcal{O}(\varepsilon))$, $k \in \{1, 2, 3\}$, the behaviour of the main terms of elementary symmetric polynomials is different:

$$e_1(\varepsilon) = a + \mathcal{O}(\varepsilon), \quad e_2(\varepsilon) = \varepsilon(a(b_1 + b_2 + b_3) + \mathcal{O}(\varepsilon)), \\ e_3(\varepsilon) = \varepsilon^2(a(b_1b_2 + b_1b_3 + b_2b_3) + \mathcal{O}(\varepsilon)), \quad e_4(\varepsilon) = \varepsilon^3(ab_1b_2b_3 + \mathcal{O}(\varepsilon)).$$

The following two lemmas generalize these observations to arbitrary orders.

**LEMMA 3.4.** *Suppose that $\lambda_1(\varepsilon), \dots, \lambda_n(\varepsilon)$ have the form*

$$(16) \qquad \lambda_1(\varepsilon) = \varepsilon^{L_1}(\tilde{\lambda}_1 + \mathcal{O}(\varepsilon)), \dots, \lambda_n(\varepsilon) = \varepsilon^{L_n}(\tilde{\lambda}_n + \mathcal{O}(\varepsilon)),$$

for some integers $0 \le L_1 \le \dots \le L_n$. Then

1. The elementary symmetric polynomials have the form

$$\begin{align} e_1(\varepsilon) &= \varepsilon^{L_1}(\tilde{e}_1 + \mathcal{O}(\varepsilon)) \tag{17} \\ e_2(\varepsilon) &= \varepsilon^{L_1+L_2}(\tilde{e}_2 + \mathcal{O}(\varepsilon)) \\ &\vdots \\ e_n(\varepsilon) &= \varepsilon^{L_1+\dots+L_n}(\tilde{e}_n + \mathcal{O}(\varepsilon)). \end{align}$$

2. If either $s = n$ or $L_s < L_{s+1}$, the main term $\tilde{e}_s$ can be expressed as

$$\tilde{e}_s = \tilde{\lambda}_1 \cdots \tilde{\lambda}_s.$$

In particular, if $s > 1$ and $\tilde{e}_{s-1} \neq 0$, then the main term $\tilde{\lambda}_s$ can be found as

$$\tilde{\lambda}_s = \frac{\tilde{e}_s}{\tilde{e}_{s-1}}.$$
We also need a generalization of Lemma 3.4 to the case of a group of equal $L_k$.

LEMMA 3.5. Let $\lambda_k(\varepsilon)$ and $L_k$ be as in Lemma 3.4 (i.e., $\lambda_k(\varepsilon)$ have the form (16) and the corresponding $e_k(\varepsilon)$ have the form (17)), and define $L_0 = -1$, $L_{n+1} = +\infty$, for an easier treatment of border cases.

If for $1 \le m \le n - s$, there is a separated group of m functions

$$ \lambda_{s+1}(\varepsilon), \dots, \lambda_{s+m}(\varepsilon) $$

of repeating degree, i.e.,

$$ (18) \qquad L_s < L_{s+1} = \dots = L_{s+m} < L_{s+m+1}, $$

then the main terms $\tilde{e}_{s+k}$, $1 \le k \le m$ in (17), are connected with ESPs of the main terms $\tilde{\lambda}_{s+k}$, $1 \le k \le m$ as follows:

$$ \tilde{e}_{s+k} = \begin{cases} e_k(\tilde{\lambda}_{s+1}, \dots, \tilde{\lambda}_{s+m}), & s = 0, \\ \tilde{\lambda}_1 \dots \tilde{\lambda}_s e_k(\tilde{\lambda}_{s+1}, \dots, \tilde{\lambda}_{s+m}) = \tilde{e}_s e_k(\tilde{\lambda}_{s+1}, \dots, \tilde{\lambda}_{s+m}), & s > 0. \end{cases} $$

In particular, if $s > 1$ and $\tilde{e}_s \ne 0$, the ESPs for the main terms are equal to

$$ e_k(\tilde{\lambda}_{s+1}, \dots, \tilde{\lambda}_{s+m}) = \frac{\tilde{e}_{s+k}}{\tilde{e}_s}, $$

hence $\tilde{\lambda}_{s+1}, \dots, \tilde{\lambda}_{s+m}$ are the roots of the polynomial (see Remark 3.2)

$$ q(\lambda) = \tilde{e}_s \lambda^m - \tilde{e}_{s+1} \lambda^{m-1} + \tilde{e}_{s+2} \lambda^{m-2} + \dots + (-1)^m \tilde{e}_{s+m}. $$

Proofs of Lemmas 3.4 to 3.5 are contained in Appendix A. Note that for analytic $\lambda_k(\varepsilon)$, Lemma 3.4 and Lemma 3.5 follow from the Newton-Puiseux theorem (see, for example, [36, section 2.1]), but we prefer to keep a more general formulation in this paper.

REMARK 3.6. Assumptions in Lemmas 3.4 to 3.5 can be relaxed (when expansions (16) are valid up to a certain order), but we keep the current statement for simplicity.

**3.3. Orders of eigenvalues and eigenvectors.** First, we recall a corollary of the Courant-Fischer “min-max” principle giving a bound on smallest eigenvalues.

THEOREM 3.7 ([19, Theorem 4.3.21]). Let $\mathbf{A} \in \mathbb{R}^{n \times n}$ be symmetric, and its eigenvalues arranged in non-increasing order $\lambda_1 \ge \lambda_2 \ge \dots \ge \lambda_n$. If there exist an m-dimensional subspace $\mathcal{L} \subset \mathbb{R}^n$ and a constant $c_1$ such that

$$ \frac{\mathbf{u}^\top \mathbf{A} \mathbf{u}}{\mathbf{u}^\top \mathbf{u}} \le c_1 $$

for all $\mathbf{u} \in \mathcal{L} \setminus \{\mathbf{0}\}$, then the smallest m eigenvalues are bounded as

$$ \lambda_n \le \lambda_{n-1} \le \dots \le \lambda_{n-m+1} \le c_1. $$

For determining the orders of eigenvalues, we will need the following corollary of Theorem 3.7 (see the proof of [42, Theorem 8]). We provide a short proof for completeness.
LEMMA 3.8. Suppose that $\mathbf{K}(\varepsilon) \in \mathbb{R}^{n \times n}$ is symmetric positive semidefinite for $\varepsilon \in [0, \varepsilon_0]$ and its eigenvalues are ordered as

$$ \lambda_1(\varepsilon) \geq \lambda_2(\varepsilon) \geq \dots \geq \lambda_n(\varepsilon) \geq 0. $$

Suppose that there exists a matrix $\mathbf{U} \in \mathbb{R}^{n \times m}$, $\mathbf{U}^{\top}\mathbf{U} = \mathbf{I}_m$, such that

$$ (19) \qquad \mathbf{U}^{\top}\mathbf{K}(\varepsilon)\mathbf{U} = \mathcal{O}(\varepsilon^L) $$

in $[0, \varepsilon_0]$, $\varepsilon_0 > 0$. Then the last $m$ eigenvalues of $\mathbf{K}(\varepsilon)$ are of order at least $L$, i.e.,

$$ (20) \qquad \lambda_j(\varepsilon) = \mathcal{O}(\varepsilon^L), \quad \text{for } n-m < j \leq n. $$

*Proof.* Assumption (19) and equivalence of matrix norms implies that

$$ \|\mathbf{U}^{\top}\mathbf{K}(\varepsilon)\mathbf{U}\|_2 \leq C\varepsilon^L. $$

for some constant $C$. Hence, we have that for any $\mathbf{z} \in \mathbb{R}^m \setminus \{0\}$,

$$ \frac{\mathbf{z}^{\top}\mathbf{U}^{\top}\mathbf{K}(\varepsilon)\mathbf{U}\mathbf{z}}{\mathbf{z}^{\top}\mathbf{U}^{\top}\mathbf{U}\mathbf{z}} \leq \|\mathbf{U}^{\top}\mathbf{K}(\varepsilon)\mathbf{U}\|_2 \leq C\varepsilon^L. $$

By choosing the range of $\mathbf{U}$ as $\mathcal{L}$ and applying Theorem 3.7, we complete the proof.□

Next, we recall a classic result on eigenvalues/eigenvectors for analytic perturbations.

THEOREM 3.9 ([22, Chapter II, Theorem 1.10]). Let $\mathbf{K}(\varepsilon)$ be an $n \times n$ matrix depending analytically on $\varepsilon$ in the neighborhood of 0 and symmetric for real $\varepsilon$. Then all the eigenvalues $\lambda_k(\varepsilon)$, $1 \leq k \leq n$ can be chosen analytic; moreover, the orthogonal projectors $\mathbf{P}_k(\varepsilon)$ on the corresponding rank-one eigenspaces can be also chosen analytic, so that

$$ (21) \qquad \mathbf{K}(\varepsilon) = \sum_{k=1}^{n} \lambda_k(\varepsilon) \mathbf{P}_k(\varepsilon) $$

is the eigenvalue decomposition of $\mathbf{K}(\varepsilon)$ in a neighborhood of 0.

We will be interested in finding the limiting rank-one projectors in Theorem 3.9, i.e.,

$$ \mathbf{P}_k = \lim_{\varepsilon \to 0} \mathbf{P}_k(\varepsilon) = \mathbf{P}_k(0), $$

where the last equality follows from analyticity of $\mathbf{P}_k(\varepsilon)$ at 0. Note that, for kernel matrices, it is impossible to retrieve the information about the limiting projectors just from $\mathbf{K}(0)$ (which is rank-one). In what follows, instead of $\mathbf{P}_k$ we will talk about limiting eigenvectors $\mathbf{p}_k$ (i.e., $\mathbf{P}_k = \mathbf{p}_k \mathbf{p}_k^{\top}$), although the latter are defined only up to a change of sign.

Armed with Theorem 3.9, we can obtain an extension of Lemma 3.8, which also gives us information about limiting eigenvectors.

LEMMA 3.10. Let $\mathbf{K}(\varepsilon)$ and $\mathbf{U} \in \mathbb{R}^{n \times m}$ be as in Lemma 3.8, and, moreover, $\mathbf{K}(\varepsilon)$ be analytic in the neighborhood of 0. Then $\mathrm{span}(\mathbf{U})$ contains all the limiting eigenvectors corresponding to the eigenvalues of the order at least $L$ (i.e., $\mathcal{O}(\varepsilon^L)$).
Proof. Let $L_1 \le \dots \le L_k \le +\infty$ be the exact orders of the eigenvalues, i.e., either

$$ \lambda_k(\varepsilon) = \varepsilon^{L_k} (\tilde{\lambda}_k + \mathcal{O}(\varepsilon)), \quad \tilde{\lambda}_k > 0, $$

or $L_k = +\infty$ in the case $\lambda_k(\varepsilon) \equiv 0$. Then, from orthogonality of the eigenvectors, we only need to prove that, $p_k^\top U = 0$ for $L_k < L$. Note that $\tilde{\lambda}_k > 0$ due to positive semidefiniteness of $K(\varepsilon)$.

Suppose that, on the contrary, there exists $k$ with $L_k < L$ such that $p_k^\top U \ne 0$, and let it be the minimal such $k$. Then from the factorization (21) we have that

$$ U^\top K(\varepsilon) U = \varepsilon^{L_k} \left( \sum_{\ell: L_\ell = L_k} \tilde{\lambda}_\ell U^\top p_k p_k^\top U + \mathcal{O}(\varepsilon) \right), $$

with a nonzero main term (from $\tilde{\lambda}_k > 0$), which contradicts the assumption (19). □

**3.4. Saddle point matrices and elementary symmetric polynomials.** In this subsection, we will be interested in determinants and elementary symmetric polynomials for so-called saddle point matrices [4]. Let $V \in \mathbb{R}^{n \times r}$ be a full column rank matrix. For a matrix $A \in \mathbb{R}^{n \times n}$, we define the saddle-point matrix as

$$ \begin{bmatrix} A & V \\ V^\top & 0 \end{bmatrix}. $$

Consider a full QR decomposition of $V \in \mathbb{R}^{n \times r}$, i.e.,

$$ (22) \qquad V = Q_{\text{full}} R_{\text{full}} = QR, \quad Q_{\text{full}} = [Q \ Q_{\perp}], \quad R_{\text{full}} = \begin{bmatrix} R \\ 0 \end{bmatrix}, $$

where $Q_{\text{full}} \in \mathbb{R}^{n \times n}$, $Q_{\perp} \in \mathbb{R}^{n \times (n-r)}$, $Q_{\text{full}}^\top Q_{\text{full}} = I$, $R$ is upper-triangular, and $Q_{\perp}^\top V = 0$.

LEMMA 3.11. For any $A \in \mathbb{R}^{n \times n}$ and $V \in \mathbb{R}^{n \times r}$, it holds that

$$ \det \begin{bmatrix} A & V \\ V^\top & 0 \end{bmatrix} = (-1)^r \det(V^\top V) \det(Q_\perp^\top A Q_\perp) = (-1)^r (\det R)^2 \det(Q_\perp^\top A Q_\perp), $$

where the matrices $R$ and $Q_\perp$ are from the full QR factorization given in (22).

The proof of Lemma 3.11 is contained in Appendix B.

We will also need to evaluate the elementary symmetric polynomials of matrices of the form $Q_\perp^\top A Q_\perp$. For a power series (or a polynomial)

$$ a(t) = a_0 + a_1 t + \dots + a_N t^N + \dots, $$

we use the following notation, standard in combinatorics, for its coefficients:

$$ [t^r] \{a(t)\} \stackrel{\text{def}}{=} a_r = \frac{a^{(r)}(0)}{r!}. $$

With this notation, the following lemma on low-rank perturbations of $A$ holds true.

LEMMA 3.12. Let $A \in \mathbb{R}^{n \times n}$ and $V \in \mathbb{R}^{n \times r}$, with the QR decomposition of $V$ given as before by eq. (22). Then for $k \ge r$ the polynomial $a(t) = e_k(A + tVV^\top)$ has degree at most $r$, and its leading coefficient is given by

$$ [t^r] \{e_k(A + tVV^\top)\} = \det(V^\top V)e_{k-r}(Q_\perp^\top A Q_\perp); $$

In particular, if $k = n$, we get

$$ [t^r] \{\det(A + tVV^\top)\} = \det(V^\top V) \det(Q_\perp^\top A Q_\perp) = (-1)^r \det \begin{bmatrix} A & V \\ V^\top & 0 \end{bmatrix}. $$
The proof of Lemma 3.12 is also contained in Appendix B.

REMARK 3.13. Alternative expressions for perturbations of elementary symmetric polynomials are available in [20, Theorem 2.16], and [21, Corollary 3.3], but they do not lead directly to the compact expression in Lemma 3.12 that we need.

## 4. Results in the 1D case.

**4.1. Smooth kernels.** We begin this section by generalizing the result of [26, Corollary 2.9] on determinants of scaled kernel matrix $K_\epsilon$ in the smooth case.

**THEOREM 4.1.** Let $K \in \mathcal{C}^{(n,n)}(\Omega \times \Omega)$ be the kernel, and (9) be its scaled version. Then for small $\epsilon$,

1. the determinant of $K_\epsilon$ from (1) has the expansion

$$ \det(K_{\epsilon}) = \epsilon^{n(n-1)} (\det(V_{\le n-1})^2 \det W_{\le n-1} + \mathcal{O}(\epsilon)). $$

2. if $K_\epsilon$ is positive semidefinite on $[0, \epsilon_0]$, $\epsilon_0 > 0$, then eigenvalues have orders

$$ \lambda_k(\epsilon) = \mathcal{O}(\epsilon^{2(k-1)}). $$

While the proof of 1) is given in [26, Corollary 2.9], and the proof of 2) for the radial analytic kernels is contained in [31, 42], we provide a short proof of Theorem 4.1 in Subsection 4.3, which also illustrates the main ideas behind other proofs in the paper.

Theorem 4.1, together with Theorem 3.1 allows us to find the main terms of the eigenvalues for analytic-in-parameter $K_\epsilon$.

**THEOREM 4.2.** Let $K \in \mathcal{C}^{(n,n)}(\Omega)$ be the kernel such that $K_\epsilon$ is symmetric positive semidefinite on $[0, \epsilon_0]$, $\epsilon_0 > 0$, and analytic in $\epsilon$ in the neighborhood of 0. Then for $s \le n$ it holds that

$$ \lambda_s = \epsilon^{2(s-1)} (\tilde{\lambda}_s + \mathcal{O}(\epsilon)), $$

where the main terms satisfy

$$ (23) \qquad \tilde{\lambda}_1 \dots \tilde{\lambda}_s = \det(V_{\le s-1}^\top V_{\le s-1}) \det(W_{\le s-1}). $$

*Proof.* First, due to analyticity and Theorem 4.1, we have that expansions (16) are valid for $L_s = 2(s-1)$. Second, the submatrices of $K_\epsilon$ are also kernel matrices (of smaller size), which, in turn can be found from Theorem 4.1. More precisely,

$$ \begin{align*} e_s(\mathbf{K}_\epsilon) &= \sum_{1 \le k_1 < \dots < k_s \le n} \det(\mathbf{K}_\epsilon(x_{k_1}, \dots, x_{k_s})) \\ &= \epsilon^{s(s-1)} \left( \det \mathbf{W}_{\le s-1} \sum_{1 \le k_1 < \dots < k_s \le n} (\det \mathbf{V}_{\le s-1}(x_{k_1}, \dots, x_{k_s}))^2 + \mathcal{O}(\epsilon) \right) \end{align*} $$

$$ (24) \qquad = \epsilon^{s(s-1)} \left( \det \mathbf{W}_{\le s-1} \det(\mathbf{V}_{\le s-1}^\top \mathbf{V}_{\le s-1}) + \mathcal{O}(\epsilon) \right), $$

where the penultimate equality follows from Theorem 4.1, and the last equality follows from the Binet-Cauchy formula. □

Finally, we employ Lemma 3.4 on relations between the main terms in (16) and (17).

**COROLLARY 4.3.** Let $K_\epsilon$ be as in Theorem 4.2, and the points in $\mathcal{X}$ be distinct.
1. If $1 < s \le n$, $\det W_{\le s-2} \ne 0$, the main term of the s-th eigenvalue can be obtained as

$$ (25) \qquad \tilde{\lambda}_s = \frac{\det(V_{\le s-1}^\top V_{\le s-1})}{\det(V_{\le s-2}^\top V_{\le s-2})} \cdot \frac{\det(W_{\le s-1})}{\det(W_{\le s-2})}. $$

2. If $1 \le s < n$, $\det W_{\le s-1} \ne 0$, then the limiting eigenvectors $p_1, \dots, p_s$ are the first $s$ columns of the $Q_{\text{full}}$ factor of the QR factorization (22) of $V_{\le s-1}$.

3. In particular, if $\det W_{\le n-2} \ne 0$, then all the main terms of the eigenvalues are given by (25), and all the limiting eigenvectors are given by the columns of the $Q$ matrix in the QR factorization of $V_{\le n-1}$.

The proof of Corollary 4.3 is also contained in Subsection 4.3.

REMARK 4.4. In Corollary 4.3, by Cramer's rule, the individual ratios in (25) can be computed in the following way:

$$ \frac{\det(V_{\le s-1}^\top V_{\le s-1})}{\det(V_{\le s-2}^\top V_{\le s-2})} = ((V_{\le s-1}^\top V_{\le s-1})^{-1})_{s,s}^{-1} = R_{s,s}^2, $$

where $R_{s,s}$ is the last ((s,s)-th) diagonal element of the $R \in \mathbb{R}^{s \times s}$ matrix in the thin QR decomposition of $V_{\le s-1}$. Similarly,

$$ \frac{\det(W_{\le s-1})}{\det(W_{\le s-2})} = ((W_{\le s-1})^{-1})_{s,s}^{-1} = C_{s,s}^2, $$

where $C_{s,s}$ is the last diagonal element of the Cholesky factor of $W_{\le s-1} = CC^\top$.

**4.2. Finite smoothness.** Next, we provide an analogue of Theorem 4.1 for a radial kernel with the order of smoothness $r$, which is smaller or equal to the number of points (i.e., Theorem 4.1 cannot be applied.).

THEOREM 4.5. For a radial kernel (10) with order of smoothness $r \le n$,

1. the determinant can be expressed as

$$ \det(K_{\varepsilon}) = \varepsilon^{n(2r-1)-r^2} (\tilde{k} + \mathcal{O}(\varepsilon)), $$

where the main term is given by

$$ (26) \qquad \tilde{k} = (-1)^r \det W_{\le r-1} \det \begin{bmatrix} f_{2r-1} D_{(2r-1)} & V_{\le r-1} \\ V_{\le r-1}^\top & 0 \end{bmatrix} $$

$$ (27) \qquad = \det W_{\le r-1} \det(V_{\le r-1}^\top V_{\le r-1}) \det(f_{2r-1} Q_{\perp}^\top D_{(2r-1)} Q_{\perp}), $$

where $Q_{\perp} \in \mathbb{R}^{n \times (n-r)}$ is the semi-orthogonal matrix such that $Q_{\perp}^\top V_{\le r-1} = 0$ (e.g., the matrix $Q_{\perp}$ in the full QR decomposition (22) of $V = V_{\le r-1}$).

2. If $K_{\varepsilon}$ is positive semidefinite on $[0, \varepsilon_0]$, $\varepsilon_0 > 0$, the eigenvalues have orders

$$ \lambda_s(\varepsilon) = \begin{cases} \mathcal{O}(\varepsilon^{2(s-1)}), & s \le r, \\ \mathcal{O}(\varepsilon^{2r-1}), & s > r. \end{cases} $$

The proof of Theorem 4.5 again is postponed to Subsection 4.4 in order to present a more straightforward corollary on eigenvalues.
As an example, we have for $r = 1$ (exponential kernel):

$$ (28) \qquad \det(\mathbf{K}_{\varepsilon}) = -\varepsilon^{n-1}nK(0,0) \det \begin{bmatrix} \mathbf{D}_{(1)} & \mathbf{1} \\ \mathbf{1}^{\top} & 0 \end{bmatrix}, $$

where $\mathbf{1}$ denotes a vector with all entries equal to 1.

Combining Theorem 4.5 with Lemmas 3.5 and 3.10, we get the following result.

**THEOREM 4.6.** Let $K$ be a kernel satisfying the assumptions of Theorem 4.5, where $\mathbf{K}_{\varepsilon}$ is positive semidefinite on $[0, \varepsilon_0]$, $\varepsilon_0 > 0$, and analytic in $\varepsilon$. Then it holds that

1. The main terms of first $r$ eigenvalues $\tilde{\lambda}_1, \dots, \tilde{\lambda}_r$ satisfy (23), for $1 \le s \le r$. In particular, if for $1 < s \le r$, $\det \mathbf{W}_{\le s-2} \ne 0$, then $\tilde{\lambda}_s$ is given by (25).

2. If the points in $\mathcal{X}$ are distinct and for $1 \le s \le r$, $\det \mathbf{W}_{\le s-1} \ne 0$, then the first $s$ limiting eigenvectors are as in Corollary 4.3. In particular, for the case $\det \mathbf{W}_{\le r-1} \ne 0$, the last $n-r$ limiting eigenvectors span the column space of $\mathbf{Q}_{\perp}$.

3. If $\det \mathbf{W}_{\le r-1} \ne 0$ and the points in $\mathcal{X}$ are distinct, then $\lambda_{r+1}, \dots, \tilde{\lambda}_n$ are the eigenvalues of

$$ f_{2r-1}(\mathbf{Q}_{\perp}^{\top} \mathbf{D}_{(2r-1)} \mathbf{Q}_{\perp}). $$

The proof of Theorem 4.6 is also contained in Subsection 4.4. Note that we obtain results on the precise locations of the last $n-r$ limiting eigenvectors in Section 8.

**4.3. Proofs for the 1D smooth case.** We need the following technical lemma.

**LEMMA 4.7.** For any upper triangular matrix $\mathbf{R} \in \mathbb{R}^{(k+1) \times (k+1)}$ it holds that

$$ \mathbf{\Delta}^{-1} \mathbf{R} \mathbf{\Delta} = \mathrm{diag}(\mathbf{R}) + \mathcal{O}(\varepsilon), $$

where $\mathrm{diag}(\mathbf{R})$ is the diagonal part of $\mathbf{R}$ and $\mathbf{\Delta} = \mathbf{\Delta}_k(\varepsilon) \in \mathbb{R}^{(k+1) \times (k+1)}$ is defined as

$$ (29) \qquad \mathbf{\Delta}_k(\varepsilon) \stackrel{\mathrm{def}}{=} \begin{bmatrix} 1 & \varepsilon & \cdots & \varepsilon^k \\ \vdots & \vdots & \ddots & \vdots \\ \mathbf{R}_k(\varepsilon) & \mathbf{R}_k(\varepsilon) & \cdots & \mathbf{R}_k(\varepsilon) \end{bmatrix} $$

*Proof.* A direct calculation gives

$$ \mathbf{\Delta}^{-1} \mathbf{R} \mathbf{\Delta} = \begin{bmatrix} R_{1,1} & \varepsilon R_{1,2} & \cdots & \varepsilon^k R_{1,k+1} \\ 0 & R_{2,2} & \cdots & \vdots \\ \vdots & \vdots & \ddots & \varepsilon R_{k,k+1} \\ 0 & \cdots & 0 & R_{k+1,k+1} \end{bmatrix}. \quad \square $$

*Proof of Theorem 4.1.* We will use a special form of the Maclaurin expansion (i.e., the Taylor expansion at 0) for bivariate functions, differentiable with respect to a “rectangular” set of multi-indices. Let us take $x, y \in \Omega$ and apply first the Maclaurin expansion with respect to $x$:

$$ (30) \qquad K(x, y) = K(0, y) + xK^{(1,0)}(0, y) + \dots + \frac{x^{n-1}K^{(n-1,0)}(0, y)}{(n-1)!} + r_n(x, y). $$
Then the Maclaurin expansion of (30) with respect to $y$ yields

$$
\begin{aligned}
K(x, y) ={}& K^{(0,0)} + \dots + \frac{y^{n-1} K^{(0,n-1)}(0, \theta_{y,1} y)}{(n-1)!} + \frac{y^n K^{(0,n)}(0, \theta_{y,1} y)}{n!} \\
& + xK^{(1,0)} + \dots + \frac{xy^{n-1} K^{(1,n-1)}(0, \theta_{y,2} y)}{(n-1)!} + \frac{xy^n K^{(1,n)}(0, \theta_{y,2} y)}{n!} \\
& \vdots \\
& + \frac{x^{n-1} K^{(n-1,0)}(0, \theta_{y,n} y)}{(n-1)!} + \dots + \frac{x^{n-1} y^{n-1} K^{(n-1,n-1)}(0, \theta_{y,n} y)}{(n-1)!(n-1)!} + \frac{x^{n-1} y^n K^{(n-1,n)}(0, \theta_{y,n} y)}{(n-1)!n!} \\
& + r_n(x, y),
\end{aligned}
\qquad
\begin{aligned}
&= K(0, y) \\
&= xK^{(1,0)}(0, y) \\
&\vdots \\
&= \frac{x^{n-1} K^{(n-1,0)}(0, y)}{(n-1)!}
\end{aligned}
$$

where $0 \le \theta_{y,1}, \dots, \theta_{y,n} \le 1$ depend only on $y$, the corresponding terms of (30) are given on the right, and a shorthand notation $K^{(i,j)} = K^{(i,j)}(0,0)$ is used. From the integral form of $r_n(x, y)$ and Taylor expansion for $g(y) \stackrel{\text{def}}{=} K^{(n,0)}(t, y)$ (as in (30)), we get

$$
\begin{aligned}
r_n(x, y) &= \int_0^x \frac{(x-t)^{n-1} K^{(n,0)}(t, y)}{(n-1)!} dt \\
&= \sum_{\ell=0}^{n-1} \left( \int_0^x \frac{(x-t)^{n-1}}{(n-1)!} \frac{y^\ell K^{(n,\ell)}(t, 0)}{\ell!} dt \right) + \int_0^x \int_0^y \frac{(x-t)^{n-1}(y-s)^{n-1} K^{(n,n)}(t, s)}{(n-1)!(n-1)!} dt ds.
\end{aligned}
$$

By the mean value theorem (as in [45, §6.3.3]), there exist $0 \le \eta_{x,1}, \dots, \eta_{x,n} \le 1$ (depending only on $x$) and $0 \le \zeta_{x,y}, \xi_{x,y} \le 1$, such that $r_n(x, y) =$

$$ \frac{x^n K^{(n,0)}(\eta_{x,1}x, 0)}{n!} + \dots + \frac{x^n y^{n-1} K^{(n,n-1)}(\eta_{x,n}x, 0)}{n!(n-1)!} + \frac{x^n y^n K^{(n,n)}(\zeta_{x,y}x, \xi_{x,y}y)}{n!n!}, $$

Next, with some abuse of notation, let $\varepsilon_0 > 0$ be such that $\varepsilon_0 x, \varepsilon_0 y \in \Omega$. Then, after replacing $(x, y)$ with $(\varepsilon x, \varepsilon y)$, $\varepsilon \in [0, \varepsilon_0]$ in the expansion of $K(x, y)$, we obtain

$$
\begin{equation} \tag{31}
\begin{split}
K(\varepsilon x, \varepsilon y) = {}& [1, \varepsilon x, \dots, (\varepsilon x)^{n-1}] \mathbf{W}_{\le n-1} [1, \varepsilon y, \dots, (\varepsilon y)^{n-1}]^\mathrm{T} \\
& + \varepsilon^n ([1, \varepsilon x, \dots, (\varepsilon x)^{n-1}] \mathbf{w}_{1,y}(\varepsilon) + \mathbf{w}_{2,x}^\mathrm{T}(\varepsilon) [1, \varepsilon y, \dots, (\varepsilon y)^{n-1}]^\mathrm{T}) + \varepsilon^{2n} w_{3,x,y}(\varepsilon),
\end{split}
\end{equation}
$$

where $\mathbf{W}_{\le n-1}$ is the Wronskian matrix defined in (4),

$$ w_{3,x,y}(\varepsilon) \stackrel{\text{def}}{=} \frac{x^n y^n K^{(n,n)}(\zeta_{\varepsilon x, \varepsilon y}(\varepsilon x), \xi_{\varepsilon x, \varepsilon y}(\varepsilon y))}{n!n!}, $$

such that $w_{3,x,y} : [0, \varepsilon_0] \to \mathbb{R}^n$ is bounded and $w_{1,y}, w_{2,x} : [0, \varepsilon_0] \to \mathbb{R}^n$ are bounded vector functions (depending only on $y$ and $x$ respectively), defined as

$$ \mathbf{w}_{1,y}(\varepsilon) \stackrel{\text{def}}{=} \frac{y^n}{n!} \left[ K^{(0,n)}(0, \theta_{\varepsilon y,1}(\varepsilon y)) \quad K^{(1,n)}(0, \theta_{\varepsilon y,2}(\varepsilon y)) \quad \dots \quad \frac{K^{(n-1,n)}(0, \theta_{\varepsilon y,n}(\varepsilon y))}{(n-1)!} \right]^\mathrm{T}, $$

$$ \mathbf{w}_{2,x}(\varepsilon) \stackrel{\text{def}}{=} \frac{x^n}{n!} \left[ K^{(n,0)}(\eta_{\varepsilon x,1}(\varepsilon x), 0) \quad K^{(n,1)}(\eta_{\varepsilon x,2}(\varepsilon x), 0) \quad \dots \quad \frac{K^{(n,n-1)}(\eta_{\varepsilon x,n}(\varepsilon x), 0)}{(n-1)!} \right]^\mathrm{T}. $$

Let $\varepsilon_0 > 0$, such that $\{\varepsilon_0 x_1, \dots, \varepsilon_0 x_n\} \subset \Omega$. From (31), the scaled kernel matrix admits for $\varepsilon \in [0, \varepsilon_0]$ the expansion

$$ (32) \quad \mathbf{K}_{\varepsilon} = \mathbf{V} \Delta \mathbf{W} \Delta \mathbf{V}^{\mathrm{T}} + \varepsilon^{k+1} (\mathbf{V} \Delta \mathbf{W}_1(\varepsilon) + \mathbf{W}_2(\varepsilon) \Delta \mathbf{V}^{\mathrm{T}}) + \varepsilon^{2(k+1)} \mathbf{W}_3(\varepsilon) $$
with $k = n - 1$, $W = W_{\le k}$, $V = V_{\le k}$, $\Delta = \Delta_k(\varepsilon)$ as in (29), and $W_1(\varepsilon)$, $W_2(\varepsilon)$, $W_3(\varepsilon)$ are $\mathcal{O}(1)$ matrices defined respectively as $W_3(\varepsilon) \stackrel{\text{def}}{=} [w_{3,x_i,x_j}(\varepsilon)]^{n,n}_{i,j=1}$,

$$W_1(\varepsilon) \stackrel{\text{def}}{=} [w_{1,x_1}(\varepsilon) \cdots w_{1,x_n}(\varepsilon)], W_2(\varepsilon) \stackrel{\text{def}}{=} [w_{2,x_1}(\varepsilon) \cdots w_{2,x_n}(\varepsilon)]^\top.$$

Let $V = QR$ be the (full) QR factorization. Then from Lemma 4.7, we have that

$$\Delta^{-1}Q^T V \Delta = \tilde{R} + \mathcal{O}(\varepsilon), \quad \text{where } \tilde{R} = \text{diag}(\mathbf{R}).$$

By pre-/post-multiplying (32) by $\Delta^{-1}Q^T$ and its transpose, we get

$$\begin{align}
\Delta^{-1}Q^T K_{\varepsilon} Q \Delta^{-1} &= (\tilde{R} + \mathcal{O}(\varepsilon))W(\tilde{R}^\top + \mathcal{O}(\varepsilon)) \nonumber \\
&\quad + \varepsilon^{k+1}(\tilde{R} + \mathcal{O}(\varepsilon))W_1(\varepsilon)Q\Delta^{-1} + \varepsilon^{k+1}\Delta^{-1}Q^T W_2(\varepsilon)(\tilde{R}^\top + \mathcal{O}(\varepsilon)) \nonumber \\
&\quad + (\varepsilon^{k+1}\Delta^{-1})Q^T W_3(\varepsilon)Q(\varepsilon^{k+1}\Delta^{-1}) \tag{33} \\
&= \tilde{R}W\tilde{R}^\top + \mathcal{O}(\varepsilon). \nonumber
\end{align}$$

where the last equality follows from $\varepsilon^{k+1}\Delta^{-1} = \mathcal{O}(\varepsilon)$. Now we are ready to prove the statements of the theorem.

1. From (33) we immediately get

$$\varepsilon^{-n(n-1)} \det K_{\varepsilon} = (\det \tilde{R})^2 \det W + \mathcal{O}(\varepsilon) = (\det \mathbf{R})^2 \det W + \mathcal{O}(\varepsilon).$$

2. Since $k = n - 1$, we can also rewrite (33) as

$$(34) \qquad Q^\top K_\varepsilon Q = \begin{bmatrix} \mathcal{O}(1) & \mathcal{O}(\varepsilon) & \cdots & \mathcal{O}(\varepsilon^k) \\ \mathcal{O}(\varepsilon) & \mathcal{O}(\varepsilon^2) & \cdots & \mathcal{O}(\varepsilon^{k+1}) \\ \vdots & \vdots & \ddots & \vdots \\ \mathcal{O}(\varepsilon^k) & \mathcal{O}(\varepsilon^{k+1}) & \cdots & \mathcal{O}(\varepsilon^{2k}) \end{bmatrix},$$

whose lower right submatrices have orders

$$(35) \qquad (Q_{s:k})^T K_{\varepsilon} Q_{s:k} = \mathcal{O}(\varepsilon^{2s}),$$

where $Q_{s:k}$ denotes the matrix of the last $n-s$ columns of $Q$. By Lemma 3.8, (35) implies the required orders of the eigenvalues. $\square$

*Proof of Corollary 4.3.* 1. Follows from (23) and Lemma 3.4.

2. From Theorem 4.2, we have $\tilde{\lambda}_{\ell} \neq 0$ for $1 \le \ell \le s$. Let $Q \in \mathbb{R}^{n \times n}$ be the factor of the QR factorization of $V_{\le n-1}$ (which can also be taken as $Q_{\text{full}}$ factor in the full QR factorization for any $V_{\le \ell-1}, \ell < n$). From (35) and Lemma 3.10, $p_1, \dots, p_\ell$ are orthogonal to $\text{span}(Q_{:,\ell:k})$, for each $1 \le \ell \le s$. Due to orthonormality of the columns of $Q$, the vectors $p_1, \dots, p_s$ must coincide with its first $s$ columns.

3. The result on eigenvalues follows from 1. For the result on eigenvectors, we note that if the first $n-1$ columns of $Q$ are limiting eigenvectors, then the last column should be the remaining limiting eigenvector. $\square$

**4.4. Proofs for the 1D finite smoothness case.**

*Proof of Theorem 4.5.* First, we will rewrite the expansion (11) in a convenient form. We will group the elements in (12) to get

$$K_{\varepsilon} = V_{\le 2r-2} \Delta_{2r-2}(\varepsilon) W_{\perp} \Delta_{2r-2}(\varepsilon) V_{\le 2r-2}^{\top} + \varepsilon^{2r-1} (f_{2r-1} D_{(2r-1)} + \mathcal{O}(\varepsilon)),$$
where $\mathbf{W}_\bullet \in \mathbb{R}^{(2r-1) \times (2r-1)}$ is the antitriangular matrix defined as

$$ (36) \qquad \mathbf{W}_\bullet = \mathbf{W}_{/0} + \mathbf{W}_{/2} + \dots + \mathbf{W}_{/2(r-2)}, $$

and $\mathbf{W}_\bullet$ are defined<fnref n="4" /> in (13). For example, in case when $r=2$

$$ \mathbf{W}_\bullet = \begin{bmatrix} f_0 & 0 & f_1 \\ 0 & -2f_1 & 0 \\ f_1 & 0 & 0 \end{bmatrix}. $$

Next, we note that $\mathbf{W}_\bullet$ can be split as

$$ \mathbf{W}_\bullet = \begin{bmatrix} \mathbf{W}_{\le r-1} & \tilde{\mathbf{W}}_1 \\ \tilde{\mathbf{W}}_2 & \end{bmatrix}, $$

where $\mathbf{W}_{\le r-1}$ is exactly the Wronskian matrix defined in (4). Therefore, since the matrices $\tilde{\mathbf{V}}_{\le 2r-2}$ and $\mathbf{\Delta}_{2r-2}(\varepsilon)$ can be partitioned as

$$ \mathbf{V}_{\le 2r-2} = [\mathbf{V}_{\le r-1} \quad \mathbf{V}_{\mathrm{rest}}], \quad \mathbf{\Delta}_{2r-2}(\varepsilon) = \begin{bmatrix} \mathbf{\Delta}_{r-1}(\varepsilon) & \varepsilon^r \mathbf{\Delta}_{r-2}(\varepsilon) \end{bmatrix}, $$

we get

$$ \mathbf{V}_{\le 2r-2} \mathbf{\Delta}_{2r-2}(\varepsilon) \tilde{\mathbf{W}} \mathbf{\Delta}_{2r-2}(\varepsilon) \mathbf{V}_{\le 2r-2}^{\top} = \mathbf{V}_{\le r-1} \mathbf{\Delta}_{r-1}(\varepsilon) \mathbf{W}_{\le r-1} \mathbf{\Delta}_{r-1}(\varepsilon) \mathbf{V}_{\le r-1}^{\top} \\ + \varepsilon^r \mathbf{V}_{\le r-1} \mathbf{\Delta}_{r-1}(\varepsilon) \underbrace{\tilde{\mathbf{W}}_1 \mathbf{\Delta}_{r-2}(\varepsilon) \mathbf{V}_{\mathrm{rest}}^{\top}}_{\mathbf{W}_1(\varepsilon)} + \varepsilon^r \underbrace{\mathbf{V}_{\mathrm{rest}} \mathbf{\Delta}_{r-2}(\varepsilon) \tilde{\mathbf{W}}_2}_{\mathbf{W}_2(\varepsilon)} \mathbf{\Delta}_{r-1}(\varepsilon) \mathbf{V}_{\le r-1}^{\top}, $$

which, after denoting $\mathbf{W} = \mathbf{W}_{\le r-1}$, $\mathbf{V} = \mathbf{V}_{\le r-1}$, $\mathbf{\Delta} = \mathbf{\Delta}_{r-1}(\varepsilon)$, gives

$$ \mathbf{K}_\varepsilon = \mathbf{V} \mathbf{\Delta} \mathbf{W} \mathbf{\Delta} \mathbf{V}^{\top} + \varepsilon^r (\mathbf{V} \mathbf{\Delta} \mathbf{W}_1(\varepsilon) + \mathbf{W}_2(\varepsilon) \mathbf{\Delta} \mathbf{V}^{\top}) + \varepsilon^{2r-1} (f_{2r-1} \mathbf{D}_{(2r-1)} + \mathcal{O}(\varepsilon)). $$

Next, we take the QR decomposition $\mathbf{V}$ (22) and consider a diagonal scaling matrix

$$ \tilde{\mathbf{\Delta}} = \begin{bmatrix} \mathbf{\Delta} & \\ & \varepsilon^{r-1} \mathbf{I}_{n-r} \end{bmatrix} \in \mathbb{R}^{n \times n}. $$

After pre-/post-multiplying $\mathbf{K}_\varepsilon$ by $\tilde{\mathbf{\Delta}}^{-1} \mathbf{Q}_{\mathrm{full}}^{\top}$ and its transpose, we get

$$ (37) \qquad \tilde{\mathbf{\Delta}}^{-1} \mathbf{Q}_{\mathrm{full}}^{\top} \mathbf{K}_{\varepsilon} \mathbf{Q}_{\mathrm{full}} \tilde{\mathbf{\Delta}}^{-1} = \begin{bmatrix} \mathbf{\Delta}^{-1} \mathbf{Q}^{\top} \mathbf{K}_{\varepsilon} \mathbf{Q} \mathbf{\Delta}^{-1} & \varepsilon^{1-r} \mathbf{\Delta}^{-1} \mathbf{Q}^{\top} \mathbf{K}_{\varepsilon} \mathbf{Q}_{\perp} \\ \varepsilon^{1-r} \mathbf{Q}_{\perp}^{\top} \mathbf{K}_{\varepsilon} \mathbf{Q} \mathbf{\Delta}^{-1} & \varepsilon^{2-2r} \mathbf{Q}_{\perp}^{\top} \mathbf{K}_{\varepsilon} \mathbf{Q}_{\perp} \end{bmatrix}, $$

where $\mathbf{Q}_{\mathrm{full}} = [\mathbf{Q} \quad \mathbf{Q}_{\perp}]$ as in (22). For the upper-left block we get, by Lemma 4.7

$$ \mathbf{\Delta}^{-1} \mathbf{Q}^{\top} \mathbf{K}_{\varepsilon} \mathbf{Q} \mathbf{\Delta}^{-1} = \mathrm{diag}(\mathbf{R}) \mathbf{W} \mathrm{diag}(\mathbf{R}) + \mathcal{O}(\varepsilon). $$

The lower-left block (which is a transpose of the upper-right one) becomes

$$ \varepsilon^{1-r} \mathbf{Q}_{\perp}^{\top} \mathbf{K}_{\varepsilon} \mathbf{Q} \mathbf{\Delta}^{-1} = \varepsilon (\mathbf{Q}_{\perp}^{\top} (\mathbf{W}_2(\varepsilon) \mathrm{diag}(\mathbf{R}) + \varepsilon^{r-1} f_{2r-1} \mathbf{D}_{(2r-1)} \mathbf{Q} \mathbf{\Delta}^{-1}) + \mathcal{O}(\varepsilon)). $$

Finally, the lower right block is

$$ \varepsilon^{2-2r} \mathbf{Q}_{\perp}^{\top} \mathbf{K}_{\varepsilon} \mathbf{Q}_{\perp} = \varepsilon (f_{2r-1} \mathbf{Q}_{\perp}^{\top} \mathbf{D}_{(2r-1)} \mathbf{Q}_{\perp} + \mathcal{O}(\varepsilon)). $$

<footnote n="4">⁴In the sum (36), $\mathbf{W}_{/2\ell}$ are padded by zeros to $(2r-1) \times (2r-1)$ matrices.</footnote>
1. Combining the blocks in (37) gives

$$
\begin{align*}
\varepsilon^{-r(r-1)-2(n-r)(r-1)} \det \mathbf{K}_{\varepsilon} &= \det(\tilde{\Delta}^{-1} \mathbf{Q}_{\text{full}}^{\top} \mathbf{K}_{\varepsilon} \mathbf{Q}_{\text{full}} \tilde{\Delta}^{-1}) = \\
&= \varepsilon^{n-r} ((\det \mathbf{R})^2 \det \mathbf{W} \det(f_{2r-1} \mathbf{Q}_{\perp}^{\top} \mathbf{D}_{(2r-1)} \mathbf{Q}_{\perp}) + \mathcal{O}(\varepsilon)), \\
&= (-1)^r \varepsilon^{n-r} \left( \det \mathbf{W} \det \left[ \begin{array}{cc} f_{2r-1} \mathbf{D}_{(2r-1)} & \mathbf{V}_{\le r-1} \\ \mathbf{V}_{\le r-1}^{\top} & 0 \end{array} \right] + \mathcal{O}(\varepsilon) \right),
\end{align*}
$$

where the last equality follows by Lemma 3.11.

2. From (37) it follows that

$$
(38) \quad \mathbf{Q}_{\text{full}}^{\top} \mathbf{K}_{\varepsilon} \mathbf{Q}_{\text{full}} =
\begin{bmatrix}
\mathcal{O}(1) & \cdots & \mathcal{O}(\varepsilon^{r-1}) & \mathcal{O}(\varepsilon^r) & \cdots & \mathcal{O}(\varepsilon^r) \\
\vdots & & \vdots & \vdots & & \vdots \\
\mathcal{O}(\varepsilon^{r-1}) & \cdots & \mathcal{O}(\varepsilon^{2(r-1)}) & \mathcal{O}(\varepsilon^{2r-1}) & \cdots & \mathcal{O}(\varepsilon^{2r-1}) \\
\mathcal{O}(\varepsilon^r) & \cdots & \mathcal{O}(\varepsilon^{2r-1}) & \mathcal{O}(\varepsilon^{2r-1}) & \cdots & \mathcal{O}(\varepsilon^{2r-1}) \\
\vdots & & \vdots & \vdots & & \vdots \\
\mathcal{O}(\varepsilon^r) & \cdots & \mathcal{O}(\varepsilon^{2r-1}) & \mathcal{O}(\varepsilon^{2r-1}) & \cdots & \mathcal{O}(\varepsilon^{2r-1})
\end{bmatrix},
$$

thus Lemma 3.8 implies the orders of the eigenvalues, as in the proof of Theorem 4.1.□

*Proof of Theorem 4.6.* 1. Note that $K \in C^{(r,r)}$, hence for $s \le r$ we can proceed as in Theorem 4.2, taking into account the fact that the orders of the eigenvalues are given by Theorem 4.5. Therefore, (23) holds true for $1 \le s \le n$.

2. For $k = n - 1$, from (38), we have that

$$
(39) \quad (\mathbf{Q}_{s:k})^{\top} \mathbf{K}_{\varepsilon} \mathbf{Q}_{s:k} = \begin{cases} \mathcal{O}(\varepsilon^{2s}), & s < r, \\ \mathcal{O}(\varepsilon^{2r-1}), & s \ge r. \end{cases}
$$

Then the statement follows from Lemma 3.10, as in the proof of Corollary 4.3.

3. Now let us consider the case $s > r$. In this case, we have

$$
\begin{align*}
\varepsilon^{-n(2r-1)+r^2} e_s(\mathbf{K}_\varepsilon) &= \varepsilon^{-n(2r-1)+r^2} \sum_{\substack{\mathcal{Y} \subset \{1, \dots, n\} \\ |\mathcal{Y}|=s}} \det(\mathbf{K}_{\varepsilon, \mathcal{Y}}) \\
&= \det \mathbf{W}_{\le r-1} \sum_{|\mathcal{Y}|=s} \det(f_{2r-1} \mathbf{Q}_{\perp, \mathcal{Y}}^{\top} \mathbf{D}_{(2r-1), \mathcal{Y}} \mathbf{Q}_{\perp, \mathcal{Y}}) \det(\mathbf{V}_{\le r-1, \mathcal{Y}}^{\top} \mathbf{V}_{\le r-1, \mathcal{Y}}) + \mathcal{O}(\varepsilon) \\
&= \det \mathbf{W}_{\le r-1} \sum_{|\mathcal{Y}|=s} [t^r] \{\det(f_{2r-1} \mathbf{D}_{(2r-1), \mathcal{Y}} + t \mathbf{V}_{\le r-1, \mathcal{Y}} \mathbf{V}_{\le r-1, \mathcal{Y}}^{\top})\} + \mathcal{O}(\varepsilon) \\
&= \det \mathbf{W}_{\le r-1} [t^r] \{e_s(f_{2r-1} \mathbf{D}_{(2r-1)} + t \mathbf{V}_{\le r-1} \mathbf{V}_{\le r-1}^{\top})\} + \mathcal{O}(\varepsilon) \\
&= \det \mathbf{W}_{\le r-1} \det(\mathbf{V}_{\le r-1}^{\top} \mathbf{V}_{\le r-1}) e_{s-r} (f_{2r-1} \mathbf{Q}_{\perp}^{\top} \mathbf{D}_{(2r-1)} \mathbf{Q}_{\perp}) + \mathcal{O}(\varepsilon),
\end{align*}
$$

where the individual steps follow from Theorem 4.5 and Lemma 3.12. Therefore, by
Lemma 3.5 we get that for all $1 \le j \le n-r$

$$
\tilde{e}_j(\tilde{\lambda}_{r+1}, \dots, \tilde{\lambda}_n) = e_j(f_{2r-1} \mathbf{Q}_{\perp}^{\top} \mathbf{D}_{(2r-1)} \mathbf{Q}_{\perp}),
$$

which together with Remark 3.2 completes the proof.

5. Multidimensional case: Preliminary facts and notations. The multidimensional case requires introducing heavier notation, which we review in this section.
<figure><img src="image_8.png" /><figcaption>FIG. 5.1. Sets of multiindices, $d=2$. Black dots: $\mathbb{P}_2$, grey dots: $\mathbb{H}_3$.</figcaption></figure>

**5.1. Multi-indices and sets.** For a multi-index $\alpha = (\alpha_1, \dots, \alpha_d) \in \mathbb{Z}_+^d$, denote

$$ \alpha! \stackrel{\text{def}}{=} \alpha_1! \cdots \alpha_d!, \qquad |\alpha| \stackrel{\text{def}}{=} \sum_{k=1}^{d} \alpha_k, $$

where $0! = 1$ by convention. For example, $|(2, 1, 3)| = 6$ and $(2, 1, 3)! = 12$.

We will frequently use the following notations

$$ \mathbb{P}_k \stackrel{\text{def}}{=} \{\alpha \in \mathbb{Z}_+^d : |\alpha| \le k\}, \qquad \mathbb{H}_k \stackrel{\text{def}}{=} \{\alpha \in \mathbb{Z}_+^d : |\alpha| = k\} = \mathbb{P}_k \setminus \mathbb{P}_{k-1}. $$

The cardinalities of these sets are given by the following well-known formulas:

$$ \#\mathbb{P}_k = \binom{k+d}{d}, \qquad \#\mathbb{H}_k = \binom{k+d-1}{d-1} = \#\mathbb{P}_k - \#\mathbb{P}_{k-1}, $$

and will be used throughout this paper.

**EXAMPLE 1.** For $d=1$, we have $\mathbb{H}_k = \{k\}$ and $\mathbb{P}_k = \{0, 1, \dots, k\}$. For $d=2$, an example is shown in Fig. 5.1.

An important class of multi-index sets is the lower sets. An $\mathcal{A} \subset \mathbb{Z}_+^d$ is called a lower set [17] if for any $\alpha \in \mathcal{A}$ all “lower” multi-indices are also in the set, i.e.,

$$ \alpha \in \mathcal{A}, \beta \le \alpha \Rightarrow \beta \in \mathcal{A}, $$

$$ \text{where } (\beta_1, \dots, \beta_d) \le (\alpha_1, \dots, \alpha_d) \iff \beta_1 \le \alpha_1, \dots, \beta_d \le \alpha_d. $$

Note that all $\mathbb{P}_k$ are lower sets.

**5.2. Monomials and orderings.** For a vector of variables $\mathbf{x} = [x_1 \cdots x_d]^\mathrm{T}$, the monomial $\mathbf{x}^\alpha$ is defined as

$$ \mathbf{x}^{\alpha} \stackrel{\text{def}}{=} x_1^{\alpha_1} \cdots x_d^{\alpha_d}. $$

**REMARK 5.1.** Note that $|\alpha|$ is the total degree of the monomial $\mathbf{x}^\alpha$. The sets of multi-indices $\mathbb{P}_k$ and $\mathbb{H}_k$ therefore correspond to the sets of monomials of degree $\le k$ and $= k$ respectively.

$$ \{\mathbf{x}^{\alpha} : |\alpha| \le k\}, \qquad \{\mathbf{x}^{\alpha} : |\alpha| = k\}. $$

In what follows, we assume that an ordering of multi-indices, i.e., all the elements in $\mathbb{Z}_+$ are linearly ordered, i.e. the relation $\prec$ is defined for all pairs of multi-indices. For example, an ordering for $d=2$ is given by

$$ (40) \quad (0,0) \prec (1,0) \prec (0,1) \prec (2,0) \prec (1,1) \prec (0,2) \prec (3,0) \prec (2,1) \prec (1,2) \prec \cdots. $$
In this paper, the ordering will not be important, as the results will not depend on the ordering. The only requirement is that the order is graded [8, Ch.2 §2], i.e.,

$$|\alpha| < |\beta| \Rightarrow \alpha < \beta.$$

### REMARK 5.2. For convenience, in case $d \ge 2$, we can use an ordering satisfying

$$(1, 0, \dots, 0) < (0, 1, \dots, 0) < \dots < (0, 0, \dots, 1),$$

such that the matrix $V_1$ defined in the next subsection is equal to

$$V_1 = [x_1 \ \dots \ x_n]^T.$$

This is not the case for graded lexicographic or reverse lexicographic [8, Ch.2 §2] orderings. Instead, a graded reflected lexicographic order [1] can be used (see (40)).

### 5.3. Multivariate Vandermonde matrices.
Next, for an ordered set of points $\mathcal{X} = \{x_1, \dots, x_n\} \subset \mathbb{R}^d$ and set of multi-indices $\mathcal{A} = \{\alpha_1, \dots, \alpha_m\} \subset \mathbb{Z}_+^d$ ordered according to the chosen ordering, we define the multivariate Vandermonde matrix

$$V_{\mathcal{A}} = V_{\mathcal{A}}(\mathcal{X}) = \left[ (x_i)^{\alpha_j} \right]_{\substack{1 \le i \le n \\ 1 \le j \le m}}.$$

We will introduce a special notation $V_{\le k} \stackrel{\text{def}}{=} V_{\mathbb{P}_k}$ and $V_k \stackrel{\text{def}}{=} V_{\mathbb{F}_k}$. Since the ordering is graded, the matrix $V_{\le k}$ can be split into blocks $V_k$ arranged by increasing degree:

$$(41) \qquad V_{\le k} = [V_0 \ V_1 \ \dots \ V_k].$$

It is easy to see that in the case $d = 1$ the definition coincides with the previous definition of the Vandermonde matrix (7).

An example of the Vandermonde matrix for $d = 2$ and the set of points

$$\mathcal{X} = \{[\begin{smallmatrix} y_1 \\ z_1 \end{smallmatrix}], [\begin{smallmatrix} y_2 \\ z_2 \end{smallmatrix}], [\begin{smallmatrix} y_3 \\ z_3 \end{smallmatrix}]\},$$

with the ordering as in (40) is given below

$$V_{\le 2} = \begin{bmatrix} 1 & y_1 & z_1 & y_1^2 & y_1 z_1 & z_1^2 \\ 1 & y_2 & z_2 & y_2^2 & y_2 z_2 & z_2^2 \\ 1 & y_3 & z_3 & y_3^2 & y_3 z_3 & z_3^2 \end{bmatrix}.$$

A special case is when the Vandermonde matrix is square, i.e., the number of monomials of degree $\le k$ is equal to the number of points:

$$(42) \qquad n = \binom{k+d}{d} = \#\mathbb{P}_k.$$

For example, $n = k + 1$ if $d = 1$ and $n = \frac{(k+2)(k+1)}{2}$ if $d = 2$.

### REMARK 5.3. Unlike the 1D case, even if all the points are different, the Vandermonde matrix $V_{\le k}$ is not necessarily invertible. For example, take the set of points on one of the axes

$$\mathcal{X} = \{[\begin{smallmatrix} -1 \\ 0 \end{smallmatrix}], [\begin{smallmatrix} 0 \\ 0 \end{smallmatrix}], [\begin{smallmatrix} 1 \\ 0 \end{smallmatrix}]\},$$

for which the Vandermonde matrix is rank-deficient:

$$V_{\le 1} = \begin{bmatrix} 1 & -1 & 0 \\ 1 & 0 & 0 \\ 1 & 1 & 0 \end{bmatrix}.$$
This effect is well-known in approximation theory [17]. If the square Vandermonde matrix is nonsingular, then the set of points $\mathcal{X}$ is called unisolvent. It is known [30, Prop. 4] that a general configuration of points (e.g., $\mathcal{X}$ are drawn from an absolutely continuous probability distribution with respect to the Lebesgue measure), is unisolvent almost surely.

**5.4. Kernels and smoothness classes.** For a function $f: \mathcal{U} \to \mathbb{R}$, $\mathcal{U} \in \mathbb{R}^d$ and a multi-index $\alpha = (\alpha_1, \dots, \alpha_d) \in \mathbb{Z}_+^d$ we use a shorthand notation for its partial derivatives (if they exist):

$$f^{(\alpha)}(\boldsymbol{x}) = \frac{\partial f^{|\alpha|}}{\partial x_1^{\alpha_1} \dots \partial x_d^{\alpha_d}}(\boldsymbol{x}).$$

It makes sense to define the smoothness classes with respect to lower sets. For a lower set $\mathcal{A} \subset \mathbb{Z}_+^d$ we define the class of functions $\mathcal{U} \to \mathbb{R}$ which have on $\mathcal{U}$ all continuous derivatives $f^{(\alpha)}$, $\alpha \in \mathcal{A}$. This class is denoted by $C^A(\mathcal{U})$.

We will consider kernels $K: \Omega \times \Omega \to \mathbb{R}$ in the class $C^{(k,k)}(\Omega) \stackrel{\text{def}}{=} C^{\mathbb{P}_k \times \mathbb{P}_k}(\Omega)$, i.e., which has all partial derivatives up to order $k$ for $\boldsymbol{x}$ and $\boldsymbol{y}$ separately.

Next, assume that we are given a kernel $K \in C^{A \times B}(\Omega)$ for lower sets $\mathcal{A}$ and $\mathcal{B}$. We will define the Wronskian matrix for this function as

$$(43) \qquad \mathbf{W}_{\mathcal{A},\mathcal{B}} = \left[ \frac{K^{(\alpha,\beta)}(\mathbf{0}, \mathbf{0})}{\alpha!\beta!} \right]_{\alpha \in \mathcal{A}, \beta \in \mathcal{B}},$$

where the rows and columns are indexed by multi-indices in $\mathcal{A}$ and $\mathcal{B}$, according to the chosen ordering.

As a special case, we will denote $\mathbf{W}_{\le k} = \mathbf{W}_{\mathbb{P}_k} \stackrel{\text{def}}{=} \mathbf{W}_{\mathbb{P}_k, \mathbb{P}_k}$. For example, for $d=2$ and $k=2$ (Example in Fig. 5.1), and ordering (40) we have $\mathbf{W}_{\le 2} =$

$$\left[ \begin{array}{ccccccc} K^{((0,0),(0,0))} & K^{((0,0),(1,0))} & K^{((0,0),(0,1))} & \frac{K^{((0,0),(2,0))}}{K^{((1,0),(2,0))}} & K^{((0,0),(1,1))} & \frac{K^{((0,0),(0,2))}}{K^{((1,0),(0,2))}} \\ K^{((1,0),(0,0))} & K^{((1,0),(1,0))} & K^{((1,0),(0,1))} & \frac{K^{((1,0),(2,0))}}{K^{((0,1),(2,0))}} & K^{((1,0),(1,1))} & \frac{K^{((1,0),(1,2))}}{K^{((0,1),(0,2))}} \\ K^{((0,1),(0,0))} & K^{((0,1),(1,0))} & K^{((0,1),(0,1))} & \frac{K^{((0,1),(2,0))}}{K^{((2,0),(2,0))}} & K^{((0,1),(1,1))} & \frac{K^{((0,1),(0,1))}}{K^{((2,0),(0,1))}} \\ \frac{K^{((2,0),(0,0))}}{2} & \frac{K^{((2,0),(1,0))}}{2} & \frac{K^{((2,0),(0,1))}}{2} & \frac{K^{((2,0),(2,0))}}{K^{((0,1),(2,0))}} & \frac{K^{((2,0),(1,1))}}{2} & \frac{K^{((2,0),(0,2))}}{K^{((1,1),(0,2))}} \\ K^{((1,1),(0,0))} & K^{((1,1),(1,0))} & K^{((1,1),(0,1))} & \frac{K^{((1,1),(2,0))}}{K^{((0,1),(2,0))}} & K^{((1,1),(1,1))} & \frac{K^{((1,1),(0,1))}}{K^{((0,1),(0,2))}} \\ \frac{K^{((0,2),(0,0))}}{2} & \frac{K^{((0,2),(1,0))}}{2} & \frac{K^{((0,2),(0,1))}}{2} & \frac{K^{((0,2),(2,0))}}{4} & \frac{K^{((0,2),(1,1))}}{2} & \frac{K^{((0,2),(0,2))}}{K^{((0,2),(0,2))}} \end{array} \right],$$

where we omit the arguments of $K^{(\alpha,\beta)}$. We will also need block-antidiagonal matrices $\mathbf{W}_s \in \mathbb{R}^{\#\mathbb{P}_s \times \#\mathbb{P}_s}$ defined as follows

$$(44) \qquad \mathbf{W}_s = \left[ \mathbf{W}_{\mathbb{H}_0, \mathbb{H}_s} \quad \dots \quad \mathbf{W}_{\mathbb{H}_s, \mathbb{H}_0} \right],$$

where $\mathbf{W}_{\mathcal{A},\mathcal{B}}$ are blocks of the Wronskian matrix defined in (43). For example

$$\mathbf{W}_{/0} = [W_{0,0}], \quad \mathbf{W}_{/1} = \left[ \mathbf{W}_{\mathbb{H}_1, \mathbb{H}_0} \quad \mathbf{W}_{\mathbb{H}_0, \mathbb{H}_1} \right],$$

and in general $\mathbf{W}_s$ contains the main block antidiagonal of $\mathbf{W}_{\le s}$.
**5.5. Taylor expansions.** The standard Taylor expansion (at 0, i.e., Maclaurin expansion) in the multivariate case is as follows [45, §8.4.4]. Let $f \in C^{k+1}(\Omega)$, where $\Omega$ is an open neighborhood of 0 containing a line segment from 0 to $x$, denoted as $[0, x]$. Then the following the following Taylor expansion holds:

$$ (45) \qquad f(x) = \sum_{\alpha \in \mathbb{P}_k} \frac{x^{\alpha}}{\alpha!} f^{(\alpha)}(0) + r_k(x), $$

where the remainder can be expressed in the Lagrange or integral forms:

$$ r_k(x) = \int_0^1 (k+1)(1-t)^k \left( \sum_{\beta \in \mathbb{H}_{k+1}} \frac{x^\beta}{\beta!} f^{(\beta)}(tx) \right) dt = \sum_{\beta \in \mathbb{H}_{k+1}} \frac{x^\beta}{\beta!} f^{(\beta)}(\theta x) $$

with $\theta \in [0, 1]$. A more general Taylor (Maclaurin) expansion has remainder in the Peano form, and requires smoothness of order one less, i.e., if $f \in C^k(\Omega)$, we have:

$$ (46) \qquad f(x) = \sum_{\alpha \in \mathbb{P}_k} \frac{x^{\alpha}}{\alpha!} f^{(\alpha)}(0) + o(\|x\|_2^k). $$

We also need a “bivariate” version for a function $f: \Omega \times \Omega \to \mathbb{R}$ (the arguments are split into two groups) such that $f \in C^{\mathbb{P}_{k_1+1} \times \mathbb{P}_{k_2+1}}(\Omega \times \Omega)$. Then we can take $x, y$ such that $[0, x], [0, y] \subset \Omega$ and apply the same steps as in the proof of Theorem 4.1 to get

$$ f(x, y) = \sum_{\alpha \in \mathbb{P}_{k_1}, \beta \in \mathbb{P}_{k_2}} \frac{x^{\alpha} y^{\beta}}{\alpha! \beta!} f^{(\alpha, \beta)}(0, 0) + \sum_{\alpha \in \mathbb{P}_{k_1}, \beta \in \mathbb{H}_{k_2+1}} \frac{x^{\alpha} y^{\beta}}{\alpha! \beta!} f^{(\alpha, \beta)}(0, \theta_{y, \alpha} y) \\ + \sum_{\alpha \in \mathbb{H}_{k_1+1}, \beta \in \mathbb{P}_{k_2}} \frac{x^{\alpha} y^{\beta}}{\alpha! \beta!} f^{(\alpha, \beta)}(\eta_{x, \beta} x, 0) + \sum_{\alpha \in \mathbb{H}_{k_1+1}, \beta \in \mathbb{H}_{k_2+1}} \frac{x^{\alpha} y^{\beta}}{\alpha! \beta!} f^{(\alpha, \beta)}(\zeta_{x, y} x, \xi_{x, y} y), $$

where $\{\eta_{x,\beta}\}_{\beta \in \mathbb{P}_{k_2}} \subset [0,1]$ depend on $x$, $\{\theta_{y,\alpha}\}_{\alpha \in \mathbb{P}_{k_1}} \subset [0,1]$ depend on $y$, and $\zeta_{x,y}, \xi_{x,y} \in [0,1]$ depend on both $x$ and $y$.

**5.6. Distance matrices and expansions of radial kernels.** Next, we consider the radial kernel (10) with order of smoothness $r$. For even $k$, as in the univariate case, we will need an expansion in the form similar to (12), which was obtained in the univariate case via the binomial expansion (8). Although we can also use the same approach in the multivariate case, we prefer to derive this expansion directly from Taylor’s formula. Let $K \in C^{2r-2}(\Omega \times \Omega)$ (not necessarily radial). Then the Taylor expansion in Peano’s form (46) yields an expansion in $\varepsilon$

$$ K(\varepsilon x, \varepsilon y) = \sum_{k=0}^{2r-2} \varepsilon^k \sum_{\substack{\alpha, \beta \in \mathbb{P}_k \\ |\alpha| + |\beta| = k}} \frac{x^{\alpha} y^{\beta}}{\alpha! \beta!} K^{(\alpha, \beta)}(0, 0) + o(\varepsilon^{2r-2}), $$

which in matrix form can be written as

$$ (47) \qquad \mathbf{K}_{\varepsilon} = \sum_{k=0}^{2r-2} \varepsilon^k \mathbf{V}_{\leq k} \mathbf{W}_{/k} \mathbf{V}_{\leq k}^{\top} + o(\varepsilon^{2r-2}). $$

For a radial kernel, the two expansions (47) and (11) coincide on $[0, \varepsilon_0]$; therefore, the
distance matrices of even order $D_{(2\ell)}$ have a compact expression as:

$$f_{2\ell}D_{(2\ell)} = V_{\le 2\ell}W_{/2\ell}V_{\le 2\ell}^{\top},$$

and moreover the expansion of $K_{\varepsilon}$ given in (12) is also valid in the multivariate case$^5$.

REMARK 5.4. For k odd, the matrices $D_{(k)}$ in the multivariate case also have the conditional positive-definiteness property (as in Lemma 2.1), except that the number of points should be $n > \#\mathbb{P}_{r-1}$.

# 6. Results in the multivariate case.

## 6.1. Determinants in the smooth case.
For a degree $k$, we will introduce a notation for the sum of all total degrees of monomials with degrees in $\mathbb{P}_k$

$$M = M(k, d) \stackrel{\text{def}}{=} \sum_{\alpha \in \mathbb{P}_k} |\alpha|, \quad \text{such that} \quad \prod_{\alpha \in \mathbb{P}_k} \varepsilon^{|\alpha|} = \varepsilon^{M(k,d)},$$

which is given by$^6$

$$M(k, d) = d \binom{k+d}{d+1}.$$

For example, if $d=1$, then $M(k, 1) = \binom{k+2}{2}$. With this notation, we can formulate the result on determinants in the multivariate case.

THEOREM 6.1. Assume that the kernel is in $C^{(k+1,k+1)}$, the scaled kernel matrix is defined by (9) and (1), and also

$$\#\mathbb{P}_{k-1} < n \le \#\mathbb{P}_k.$$

1. If $n = \#\mathbb{P}_k$, then

$$\det \mathbf{K}_{\varepsilon} = \varepsilon^{2M(k,d)} (\det \mathbf{W}_{\le k} (\det(\mathbf{V}_{\le k}))^2 + \mathcal{O}(\varepsilon)).$$

2. If $n < \#\mathbb{P}_k$, for $\ell = \#\mathbb{P}_k - n$, we have

$$\det \mathbf{K}_{\varepsilon} = \varepsilon^{2(M(k,d)-k\ell)} (\det(\mathbf{Y}\mathbf{W}_{\le k}\mathbf{Y}^{\top}) \det(\mathbf{V}_{\le k-1}^{\top}\mathbf{V}_{\le k-1}) + \mathcal{O}(\varepsilon)),$$

where $\mathbf{Y} \in \mathbb{R}^{n \times \#\mathbb{P}_k}$ is defined as

$$\mathbf{Y} = \begin{bmatrix} I_{\#\mathbb{P}_{k-1}} & \mathbf{Q}_{\perp}^{\top}\mathbf{V}_k \end{bmatrix},$$

$\mathbf{V}_k$ is the Vandermonde matrix (41) for monomials of degree $k$, and $\mathbf{Q}_{\perp} \in \mathbb{R}^{n \times \ell}$ comes from the full QR decomposition of $\mathbf{V} = \mathbf{V}_{\le k-1}$ (see (22)).

3. If $\mathbf{K}_{\varepsilon}$ is positive semidefinite on $[0, \varepsilon_0]$, the eigenvalues split into $k+1$ groups

$$(48) \qquad \underbrace{\tilde{\lambda}_{0,0}}_{\mathcal{O}(1)} \underbrace{\{\tilde{\lambda}_{1,j}\}_{j=1}^d}_{\mathcal{O}(\varepsilon^2)} \dots \underbrace{\{\tilde{\lambda}_{s,j}\}_{j=1}^{\#\mathbb{H}_s}}_{\mathcal{O}(\varepsilon^{2s})} \dots \underbrace{\{\tilde{\lambda}_{k,j}\}_{j=1}^{\#\mathbb{H}_{k-\ell}}}_{\mathcal{O}(\varepsilon^{2k})}.$$

<footnote>$^5$Note that there is an equivalent way of obtaining the expansion of the distance matrices in terms of monomials, simply by writing $\|\mathbf{x} - \mathbf{y}\|^{2p} = (\sum_{i=1}^d (x_i - 2y_i x_i + y_i))^p$ and expanding.</footnote>

<footnote>$^6$See [27, eq. (3.19)-(3.20)], where $M(k, d)$ is given in a slightly different form.</footnote>
The proof of Theorem 6.1 is postponed to Subsection 6.3.

From Theorem 6.1, we can also get a result on eigenvalues and eigenvectors. For this, we partition the $Q_{\text{full}}$ matrix in the full QR factorization of $V_{\le k-1}$ as

$$ (49) \qquad Q_{\text{full}} = [Q_0 \ Q_1 \ \dots \ Q_k] $$

with $Q_s \in \mathbb{R}^{n \times \#\mathbb{H}_s}$, $0 \le s < k$, and $Q_k \in \mathbb{R}^{n \times (\#\mathbb{H}_k - \ell)}$.

**THEOREM 6.2.** Let $K$ be as in Theorem 6.1, such that $K_\epsilon$ is symmetric positive semidefinite on $[0, \epsilon_0]$ and analytic in $\epsilon$ in a neighborhood of 0. Then the eigenvalues in the groups have the form

$$ \lambda_{s,j} = \epsilon^{2s} (\tilde{\lambda}_{s,j} + \mathcal{O}(\epsilon)), $$

where $\tilde{\lambda}_{0,0} = nK^{(0,0)}$ and the other main terms are given as follows.

1. For $1 \le s < k$, if $\det W_{\le s-1} \ne 0$ and $V_{\le s-1}$ is full rank, then

$$ (50) \qquad \tilde{\lambda}_{s,1} \cdots \tilde{\lambda}_{s,\#\mathbb{H}_s} = \frac{\det(V_{\le s}^\top V_{\le s}) \det(W_{\le s})}{\det(V_{\le s-1}^\top V_{\le s-1}) \det(W_{\le s-1})} $$

2. For any $1 \le s \le k$, if $\det W_{\le s-1} \ne 0$ and $V_{\le s-1}$ is full rank, the main terms $\tilde{\lambda}_{s,1}, \dots, \tilde{\lambda}_{s,\#\mathbb{H}_s}$ (or $\tilde{\lambda}_{k,1}, \dots, \tilde{\lambda}_{k,\#\mathbb{H}_k-\ell}$ if $s=k$) are the eigenvalues of

$$ (51) \qquad Q_s^\top V_s \tilde{W} V_s^\top Q_s, $$

where $\tilde{W} = W_\perp - W_\perp(W_{\le s-1})^{-1}W_\perp$ is the Schur complement coming from the following partition of the Wronskian:

$$ W_{\le s} = \begin{bmatrix} W_{\le s-1} & W_\perp \\ W_\perp & W_\perp \end{bmatrix} = \begin{bmatrix} W_{\le s-1} & W_\perp, s \\ W_\perp, s & W_\perp \end{bmatrix}, $$

3. For $0 \le s < k$, if $\det W_{\le s} \ne 0$ and $V_{\le s}$ is full rank, then the limiting eigenvectors from the s-th group $p_{s,1}, \dots, p_{s,\#\mathbb{H}_s}$ span the column space of $Q_s$. Moreover, if $\det W_{\le k-1} \ne 0$, the remaining eigenvectors span the column space of $Q_k$.

The proof of Theorem 6.2 is postponed to Subsection 6.4.

**Theorem 6.2** does not give information on the precise location of limiting eigenvectors in each group. We formulate the following conjecture, which we validated numerically.

**CONJECTURE 1.** For $1 \le s \le k$, if $\det W_{\le s} \ne 0$ and $V_{\le s}$ is full rank, the limiting eigenvectors in the $\mathcal{O}(\epsilon^{2s})$ group are the columns of $Q_s U_s$, where $U_s$ contains the eigenvectors$^7$ of the matrix $Q_s^\top V_s \tilde{W} V_s^\top Q_s$ from (51).

**6.2. Finite smoothness case.** We prove a generalization of Theorem 4.5 to the multivariate case.

**THEOREM 6.3.** For small $\epsilon$ and a radial kernel (10) with order of smoothness $r$: 1. the determinant of $K_\epsilon$ in the case $n = \#\mathbb{P}_{r-1} + N$ given in (1) has the expansion

$$ \det(K_\epsilon) = \epsilon^{2M(r-1,d)+(2r-1)N} (\tilde{k} + \mathcal{O}(\epsilon)), $$

<footnote>$^7$There is a usual issue of ambiguous definition of $U_s$ if the matrix has repeating eigenvalues.</footnote>
where $\tilde{k}$ has exactly the same expression as in (26) or (27) (with $V_{\le r}$ and $D_{(2r-1)}$ replaced with their multivariate counterparts).

2. If $K_{\epsilon}$ is positive semidefinite on $[0, \epsilon_0]$, the eigenvalues are split into $r+1$ groups

$$ \underbrace{\tilde{\lambda}_{0,0}}_{\mathcal{O}(1)}, \underbrace{\{\tilde{\lambda}_{1,j}\}_{j=1}^d}_{\mathcal{O}(\epsilon^2)}, \dots, \underbrace{\{\tilde{\lambda}_{r-1,j}\}_{j=1}^{\#\mathbb{H}_{r-1}}}_{\mathcal{O}(\epsilon^{2(r-1)})}, \underbrace{\{\tilde{\lambda}_{k,j}\}_{j=1}^N}_{\mathcal{O}(\epsilon^{2r-1})}. $$

3. In the analytic in $\epsilon$ case, the main terms for the first $r$ groups are the same as in Theorem 6.2. For the last group, if $\det W_{\le r-1} \ne 0$ and $V_{\le r-1}$ is full rank, the main terms are the eigenvalues of

$$ f_{2r-1}(Q_{\perp}^{\top} D_{(2r-1)} Q_{\perp}), $$

where $Q_{\perp}$ comes from the full QR factorization (22) of $V_{\le r-1}$.

4. For $0 \le s < r$, the subspace spanned by the limiting eigenvectors for the $\mathcal{O}(\epsilon^{2s})$ group of eigenvalues are as in Theorem 6.2. If $\det W_{\le r-1} \ne 0$ and $V_{\le r-1}$ is full rank, the eigenvectors for the last group of $\mathcal{O}(\epsilon^{2r-1})$ eigenvalues span the column space of $Q_{\perp}$.

The proof of Theorem 6.3 is postponed to Subsection 6.4. Note that we obtain a stronger result on the precise location of the last group of eigenvectors in Section 8.

**6.3. Determinants in the smooth case.** Before proving Theorem 6.1, we again need a technical lemma, which is an analogue of (4.7).

**LEMMA 6.4.** Let $\mathbf{R}$ be an upper-block-triangular matrix

$$ \mathbf{R} = \begin{bmatrix} \mathbf{R}_{1,1} & \mathbf{R}_{1,k} & \cdots & \mathbf{R}_{1,k+1} \\ 0 & \ddots & & \vdots \\ \vdots & \ddots & \ddots & \mathbf{R}_{k,k+1} \\ 0 & \cdots & 0 & \mathbf{R}_{k+1,k+1} \end{bmatrix}, $$

where the blocks $R_{i,j} \in \mathbb{R}^{M_i \times N_j}$ are not necessarily square. Then it holds that

$$ \begin{bmatrix} I_{M_1} & & & & \\ & \varepsilon^{-1} I_{M_2} & & & \\ & & \ddots & & \\ & & & \varepsilon^{-k} I_{M_{k+1}} & \\ & & & & \varepsilon^{k} I_{N_{k+1}} \end{bmatrix} \mathbf{R} \begin{bmatrix} I_{N_1} & & & & \\ & \varepsilon I_{N_2} & & & \\ & & \ddots & & \\ & & & \varepsilon^{k} I_{N_{k+1}} & \\ & & & & \varepsilon^{-k} I_{M_{k+1}} \end{bmatrix} = \text{blkdiag}(\mathbf{R}) + \mathcal{O}(\varepsilon), $$

where $\text{blkdiag}(\mathbf{R})$ is just the block-diagonal part of $\mathbf{R}$:

$$ \text{blkdiag}(\mathbf{R}) = \begin{bmatrix} \mathbf{R}_{1,1} \\ & \ddots \\ & & \mathbf{R}_{k+1,k+1} \end{bmatrix}. $$

*Proof.* The proof is analogous to that of Lemma 4.7. □

**Proof of Theorem 6.1.** First, we fix a degree-compatible ordering of multi-indices

$$ \mathbb{P}_k = \{\alpha_1, \dots, \alpha_{\#\mathbb{P}_k}\}, $$
and denote the $\#\mathbb{P}_k \times \#\mathbb{P}_k$ matrix (note that $\#\mathbb{P}_k = n + \ell$)

$$ (52) \quad \Delta = \Delta_k(\varepsilon) \stackrel{\text{def}}{=} \begin{bmatrix} \varepsilon^{\alpha_1} \\ \vdots \\ \varepsilon^{\alpha_{\#\mathbb{P}_k}} \end{bmatrix} = \begin{bmatrix} \varepsilon I_{\#\mathbb{H}_1} \\ \vdots \\ \varepsilon^k I_{\#\mathbb{H}_k} \end{bmatrix}, $$

and by $E_n$ the principal $n \times n$ submatrix of $\Delta$. Note that their determinants are

$$ (53) \quad \det(\Delta) = \varepsilon^{M(k,d)} \quad \text{and} \quad \det(E_n) = \varepsilon^{M(k,d)-\ell k}. $$

From the bimultivariate Taylor expansion, as in the proof of Theorem 4.1, we get:

$$ K(\varepsilon\mathbf{x}, \varepsilon\mathbf{y}) = [(\varepsilon\mathbf{x})^{\alpha_1}, \dots, (\varepsilon\mathbf{x})^{\alpha_{n+\ell}}] W [(\varepsilon\mathbf{y})^{\alpha_1}, \dots, (\varepsilon\mathbf{y})^{\alpha_{n+\ell}}]^T $$

$$ (54) \quad + \varepsilon^{k+1}[(\varepsilon\mathbf{x})^{\alpha_1}, \dots, (\varepsilon\mathbf{x})^{\alpha_{n+\ell}}]w_{1,\mathbf{y}}(\varepsilon) + \varepsilon^{k+1}w_{2,\mathbf{x}}(\varepsilon)^T[(\varepsilon\mathbf{y})^{\alpha_1}, \dots, (\varepsilon\mathbf{y})^{\alpha_{n+\ell}}]^T \\ + \varepsilon^{2(k+1)}w_{3,\mathbf{x},\mathbf{y}}(\varepsilon), $$

where $w_1$, $w_2$ are bounded (and continuous) $[0, \varepsilon_0] \to \mathbb{R}^n$ vector functions depending on $\mathbf{y}$ and $\mathbf{x}$ respectively, and $w_3$ is a bounded (and continuous) function $[0, \varepsilon_0] \to \mathbb{R}$.

Let $\varepsilon_0 > 0$, such that $\{\varepsilon_0\mathbf{x}_1, \dots, \varepsilon_0\mathbf{x}_n\} \in \Omega$ for all $i$. From (54), the scaled kernel matrix admits for $0 \le \varepsilon \le \varepsilon_0$ the expansion

$$ (55) \quad \mathbf{K}_{\varepsilon} = \mathbf{V}_{\le k} \Delta \mathbf{W} \Delta \mathbf{V}_{\le k}^{\top} + \varepsilon^n (\mathbf{V}_{\le k} \Delta \mathbf{W}_1(\varepsilon) + \mathbf{W}_2(\varepsilon) \Delta \mathbf{V}_{\le k}^{\top}) + \varepsilon^{2n} \mathbf{W}_3(\varepsilon), $$

where $\mathbf{W}_1(\varepsilon), \mathbf{W}_2(\varepsilon), \mathbf{W}_3(\varepsilon) = \mathcal{O}(1)$, $\mathbf{W} = \mathbf{W}_{\le k}, \mathbf{V} = \mathbf{V}_{\le k}$.

Next, we take the full QR decomposition $\mathbf{V} = \mathbf{Q}_{\text{full}}\mathbf{R}_{\text{full}}$ of $\mathbf{V} = \mathbf{V}_{\le k-1}$, partitioned as in (22), so that $\mathbf{R} \in \mathbb{R}^{\#\mathbb{P}_{k-1} \times \#\mathbb{P}_{k-1}}$ and $\mathbf{Q}_{\perp} \in \mathbb{R}^{n \times (n-\#\mathbb{P}_{k-1})}$. Note that

$$ \mathbf{Q}_{\text{full}}^{\top} \mathbf{V}_{\le k} = \begin{bmatrix} \mathbf{R} & \mathbf{Q}_{\perp}^{\top} \mathbf{V}_k \\ \mathbf{O} & \mathbf{Q}_{\perp}^{\top} \mathbf{V}_k \end{bmatrix}, $$

and by Lemma 6.4 we have

$$ \mathbf{E}_n^{-1} \mathbf{Q}_{\text{full}}^{\top} \mathbf{V}_{\le k} \mathbf{\Delta} = \underbrace{\left[ \text{blkdiag}(\mathbf{R}) \quad \mathbf{Q}_{\perp}^{\top} \mathbf{V}_k \right]}_{\tilde{\mathbf{R}}} + \mathcal{O}(\varepsilon). $$

Next, we pre-/post- multiply (55) by $\mathbf{E}_n^{-1} \mathbf{Q}_{\text{full}}^{\top}$ and its transpose, to get (as in (33)),

$$ (56) \quad \mathbf{E}_n^{-1} \mathbf{Q}_{\text{full}}^{\top} \mathbf{K}_{\varepsilon} \mathbf{Q}_{\text{full}} \mathbf{E}_n^{-1} = \tilde{\mathbf{R}} \mathbf{W} \tilde{\mathbf{R}}^{\top} + \mathcal{O}(\varepsilon). $$

Finally, we prove the statements of the theorem.

1. From (53) we have that

$$ \varepsilon^{-2(M(k,d)-\ell k)} \det \mathbf{K}_{\varepsilon} = \det(\tilde{\mathbf{R}} \mathbf{W} \tilde{\mathbf{R}}^{\top}). $$

Thus, if $n = \#\mathbb{P}_k$, then we have

$$ \det(\tilde{\mathbf{R}} \mathbf{W} \tilde{\mathbf{R}}^{\top}) = (\det(\tilde{\mathbf{R}}))^{2} \det(\mathbf{W}) = \det(\mathbf{V}_{\le k})^{2} \det(\mathbf{W}), $$

where the last equality follows from the fact that

$$ \det(\tilde{\mathbf{R}}) = \det(\text{blkdiag}(\mathbf{R})) \det(\mathbf{Q}_{\perp}^{\top} \mathbf{V}_k) = \det(\mathbf{Q}_{\text{full}}^{\top} \mathbf{V}_{\le k}) = \det(\mathbf{V}_{\le k}), $$

because $\mathbf{Q}_{\text{full}}^{\top} \mathbf{V}_{\le k}$ is block-triangular.
2. For $n < \#\mathbb{P}_k$, we note that

$$ \tilde{R} = \begin{bmatrix} \text{blkdiag } R & I_\ell \end{bmatrix} Y, $$

hence

$$ \det(\tilde{R}W\tilde{R}^\mathrm{T}) = (\det(R))^2 \det(YWY^\mathrm{T}) = \det(V_{\le k-1}^\mathrm{T}V_{\le k-1}) \det(YWY^\mathrm{T}). $$

3. Finally, as in the proof of Theorem 4.1, (56) implies that (35) holds as well in the multivariate case; this, together with Theorem 3.7 completes the proof.□

### 6.4. Individual eigenvalues, eigenvectors, and finite smoothness.

*Proof of Theorem 6.2.* 1–2. Choose a subset $\mathcal{Y}$ of $\mathcal{X}$ of size $m$, $\#\mathbb{P}_{s-1} < m \le \#\mathbb{P}_s$, $s \le k$. Then we have that

$$ K_{\varepsilon,y} = \varepsilon^{2(M-s(\#\mathbb{P}_s-m))} (\det(YW_{\le s}Y^\mathrm{T}) \det(V_{\le s-1,y}^\mathrm{T}V_{\le s-1,y}) + \mathcal{O}(\varepsilon)), $$

and

$$ Y = \begin{bmatrix} I_{\#\mathbb{P}_{s-1}} & R_{s,y} \end{bmatrix}, $$

where $R_{s,y} = Q_{\perp,y}^\mathrm{T}V_s$. In particular

$$ \begin{aligned} \det(YW_{\le s}Y^\mathrm{T}) &= \det \begin{bmatrix} W_{\le s-1} & W^\mathrm{T} R_{s,y}^\mathrm{T} \\ R_{s,y} W_\perp & R_{s,y} W_\perp R_{s,y}^\mathrm{T} \end{bmatrix} \\ &= \det W_{\le s-1} \det(R_{s,y} W_\perp R_{s,y}^\mathrm{T} - R_{s,y} W_\perp (W_{\le s-1})^{-1} W^\mathrm{T} R_{s,y}^\mathrm{T}) \\ &= \det W_{\le s-1} \det(Q_\perp^\mathrm{T} V_s \tilde{W} V_s^\mathrm{T} Q_\perp, y), \end{aligned} $$

hence by Lemma 3.11

$$ \begin{aligned} & \det(YW_{\le s}Y^\mathrm{T}) \det(V_{\le s-1,y}^\mathrm{T}V_{\le s-1,y}) \\ &= \det W_{\le s-1}[\gamma^{\#\mathbb{P}_{s-1}}] \{\det(V_s,y \tilde{W} V_{s,y}^\mathrm{T} + \gamma V_{\le s-1,y} V_{\le s-1,y}^\mathrm{T})\} \end{aligned} $$

and therefore

$$ \begin{aligned} \tilde{e}_m &= \det W_{\le s-1}[\gamma^{\#\mathbb{P}_{s-1}}] \{e_m(V_s \tilde{W} V_s^\mathrm{T} + \gamma V_{\le s-1} V_{\le s-1}^\mathrm{T})\} \\ &= \det W_{\le s-1} \det(V_{\le s-1}^\mathrm{T} V_{\le s-1}) e_{m-\#\mathbb{P}_{s-1}}(Q_{s:k}^\mathrm{T} V_s \tilde{W} V_s^\mathrm{T} Q_{s:k}) \\ (57) \qquad \qquad \qquad \qquad &= \det W_{\le s-1} \det(V_{\le s-1}^\mathrm{T} V_{\le s-1}) e_{m-\#\mathbb{P}_{s-1}}(Q_s^\mathrm{T} V_s \tilde{W} V_s^\mathrm{T} Q_s), \end{aligned} $$

where $Q_{s:k} \stackrel{\mathrm{def}}{=} [Q_s \quad \dots \quad Q_k]$, the penultimate equality follows from Lemma 3.11, and the last equality follows from the fact that only the top block of $Q_{s:k}^\mathrm{T} V_s$ is nonzero. The rest of the proof follows from Lemma 3.5 to (57).

3. We repeat the same steps as in the proof of Corollary 4.3 (for groups of limiting eigenvectors). Since, from the proof of Theorem 6.1, for any $0 \le \ell \le s$, the order of $(Q_{\ell+1:k})^\mathrm{T} K_\varepsilon Q_{\ell+1:k}$ is $\mathcal{O}(\varepsilon^{2(\ell+1)})$ (as in (35)), and the order of eigenvalues in $\mathcal{O}(1), \dots, \mathcal{O}(\varepsilon^{2s})$ groups is exact, all the eigenvectors

$$ \{\mathbf{p}_{\ell,1}, \dots, \mathbf{p}_{\ell,\#\mathbb{H}_\ell}\}_{0\le\ell\le s}, $$

must be orthogonal to $\mathrm{span}(Q_{\ell+1:k})$, which proves the first part of the statement. If, moreover, $s = k-1$, then the last block of eigenvectors (corresponding to eigenvalues of order $\mathcal{O}(\varepsilon^{2k})$) must be contained in $\mathrm{span}(Q_k)$. □
Proof of Theorem 6.3. 1. The proof repeats that of Theorem 4.5 with the following minor modifications (in order to take into account the multivariate case):

• the matrix $W_{\bullet} \in \mathbb{R}^{\#\mathbb{P}_{2r-2} \times \#\mathbb{P}_{2r-2}}$ is defined as

$$W_{\bullet, \alpha, \beta} = \begin{cases} \frac{K^{(\alpha, \beta)}(0,0)}{\alpha!\beta!}, & |\alpha| + |\beta| \le 2r, \\ 0, & |\alpha| + |\beta| > 2r; \end{cases}$$

i.e., in the sum (36) the matrices $W_{/s}$ are defined according to (44);

• the matrix $\Delta_k(\varepsilon)$ is defined in (52), and Lemma 6.4 is used instead of Lemma 4.7;

• the extended diagonal scaling matrix is

$$\tilde{\Delta} = \begin{bmatrix} \Delta_{r-1} & \\ & \varepsilon^{r-1} I_N \end{bmatrix} \in \mathbb{R}^{n \times n},$$

where $N = n - \#\mathbb{P}_{r-1}$;

• the last displayed formula in the proof of Theorem 4.5 becomes

$$\begin{aligned} \varepsilon^{-2M(r-1,d)-2N(r-1)} \det \mathbf{K}_{\varepsilon} &= \det(\tilde{\Delta}^{-1} \mathbf{Q}_{\text{full}}^{\top} \mathbf{K}_{\varepsilon} \mathbf{Q}_{\text{full}} \tilde{\Delta}^{-1}) \\ &= \varepsilon^N ((\det \mathbf{R})^2 \det \mathbf{W} \det(f_{2r-1} \mathbf{Q}_{\perp}^{\top} \mathbf{D}_{(2r-1)} \mathbf{Q}_{\perp}) + \mathcal{O}(\varepsilon)). \end{aligned}$$

2. Again, as in the proof of Theorem 4.5, the matrix $\mathbf{Q}_{s:k}^{\top} \mathbf{K}_{\varepsilon} \mathbf{Q}_{s:k}$ has order given in (39), hence the orders of the eigenvalues follow from Theorem 3.7.

3. The proof of this statement repeats the proof of Theorem 4.6 without changes.

4. The last statement follows from combining (39) with Theorem 3.7, and proceeding as in the proof of the corresponding statement in Theorem 6.2. □

**7. Perturbation theory: a summary for Hermitian matrices.** This section contains a summary of facts from [22, Ch. II] to deal with analytic perturbations of self-adjoint operators in a finite-dimensional vector space (i.e., Hermitian matrices). Formally, and in keeping with the notation used in [22], we assume that we are given a matrix-valued function depending on $\varepsilon$ such that

$$(58) \qquad \mathbf{T}(\varepsilon) = \mathbf{T}^{(0)} + \varepsilon \mathbf{T}^{(1)} + \varepsilon^2 \mathbf{T}^{(2)} + \dots,$$

where we assume that the matrices $\mathbf{T}^{(k)} \in \mathbb{C}^{n \times n}$ are Hermitian. In a neighborhood of 0, $0 \in \mathcal{D}_0 \subset \mathbb{C}$, $\mathbf{T}(\varepsilon)$ has $s \le n$ semi-simple eigenvalues generically (i.e., except a finite number of exceptional points). For simplicity of presentation$^8$, we assume that $s = n$, which is the case if there exists $\varepsilon_1 \in \mathcal{D}_0$ having all distinct eigenvalues. The interesting case (considered in this paper) is when $\varepsilon = 0$ is an exceptional point, i.e. $\mathbf{T}(0) = \mathbf{T}^{(0)}$ has multiple eigenvalues (e.g., a low-rank matrix with a multiple eigenvalue 0).

**7.1. Perturbation of eigenvalues and group eigenprojectors.** Since all matrices are Hermitian, by [22, Theorem 1.10, Ch. II] (see Theorem 3.9), the eigenvalues $\lambda_1(\varepsilon), \dots, \lambda_n(\varepsilon)$ and the rank-one projectors on the corresponding eigenspaces $\mathbf{P}_1(\varepsilon), \dots, \mathbf{P}_n(\varepsilon)$ are holomorphic functions of $\varepsilon$ in a neighborhood of 0, $0 \in \mathcal{D} \subset \mathbb{C}$.

<footnote>$^8$We can also consider the general case if needed.</footnote>
REMARK 7.1. If the matrix $\mathbf{T}^{(0)}$ has $d$ multiple eigenvalues $\{\mu_k\}_{k=1}^d$ with multiplicities $m_1, \dots, m_d$, i.e., after proper reordering, at $\varepsilon = 0$

$$ (\lambda_1(0), \dots, \lambda_n(0)) = \underbrace{(\mu_1, \dots, \mu_1)}_{m_1 \text{ times}}, \underbrace{(\mu_2, \dots, \mu_2)}_{m_2 \text{ times}}, \dots, \underbrace{(\mu_d, \dots, \mu_d)}_{m_d \text{ times}} $$

then the projectors on the invariant subspaces associated with $\mu_1, \dots, \mu_s$ are sums of the corresponding rank-one projectors on the eigenspaces:

$$ \begin{align*} \mathbf{P}_{\mu,1} &= \mathbf{P}_1(0) + \dots + \mathbf{P}_{m_1}(0) \\ \mathbf{P}_{\mu,2} &= \mathbf{P}_{m_1+1}(0) + \dots + \mathbf{P}_{m_1+m_2}(0) \\ \vdots \\ \mathbf{P}_{\mu,d} &= \mathbf{P}_{m_1+\dots+m_{d-1}+1}(0) + \dots + \mathbf{P}_{n}(0) \end{align*} $$

In this paper, our aim is to obtain a limiting eigenstructure at $\varepsilon = 0$. In case of multiple eigenvalues, this information cannot be retrieved from the spectral decomposition of $\mathbf{T}^{(0)}$ alone (we can only retrieve the group projectors $\mathbf{P}_{\mu,j}$ from the spectral decomposition of $\mathbf{T}^{(0)}$). In what follows, we look in details at perturbation expansions in order to find individual $\mathbf{P}_k(0)$.

As shown in [22, Ch. II], we can consider perturbations of a possibly multiple eigenvalue. Let $\lambda$ be an eigenvalue of $\mathbf{T}^{(0)}$ of multiplicity $m$, and $\mathbf{P}$ is the corresponding orthogonal projector on the $m$-dimensional eigenspace. The projector on the perturbed $m$-dimensional invariant subspace is an analytic matrix-valued function

$$ \mathbf{P}(\varepsilon) = \sum_{k=0}^{\infty} \varepsilon^k \mathbf{P}^{(k)}, $$

with the coefficients given by

$$ \begin{gathered} \mathbf{P}^{(0)} = \mathbf{P}, \quad \mathbf{P}^{(1)} = -\mathbf{P}\mathbf{T}^{(1)}\mathbf{S} - \mathbf{S}\mathbf{T}^{(1)}\mathbf{P}, \\ (59) \quad \mathbf{P}^{(k)} = -\sum_{p=1}^{k} (-1)^p \sum_{\substack{\nu_1+\dots+\nu_p=k \\ k_1+\dots+k_{p+1}=p \\ \nu_j \ge 1, k_j \ge 0}} \mathbf{S}^{(k_1)}\mathbf{T}^{(\nu_1)}\mathbf{S}^{(k_2)}\dots\mathbf{S}^{(k_p)}\mathbf{T}^{(\nu_p)}\mathbf{S}^{(k_{p+1})}, \end{gathered} $$

where $\mathbf{S} = (\mathbf{T} - \lambda\mathbf{I})^{\dagger}$, $\mathbf{S}^{(0)} = -\mathbf{P}$ and $\mathbf{S}^{(k)} = \mathbf{S}^k$.

7.2. Reduction and splitting the groups. In order to find the individual projectors of the eigenspaces corresponding to a multiple $\lambda$, and the expansion of the corresponding eigenvalues, the following reduction (or splitting) procedure [22, Ch. II, §2.3] can be applied, which localises the matrix to the $m$-dimensional subspace corresponding to the perturbations of $\lambda$.

We first define the eigennilpotent matrix $\mathbf{D}(\varepsilon)$ as

$$ \mathbf{D}(\varepsilon) = (\mathbf{T}(\varepsilon) - \lambda\mathbf{I})\mathbf{P}(\varepsilon) = \mathbf{P}(\varepsilon)(\mathbf{T}(\varepsilon) - \lambda\mathbf{I}) = \mathbf{P}(\varepsilon)(\mathbf{T}(\varepsilon) - \lambda\mathbf{I})\mathbf{P}(\varepsilon), $$

which from [22, Ch. II, §2.2] has an expansion

$$ \mathbf{D}(\varepsilon) = 0 + \sum_{k=1}^{\infty} \varepsilon^k \tilde{\mathbf{T}}^{(k)}, $$
where the expressions for $\tilde{T}^{(k)}$ are as follows:

$$ 
\begin{align*}
\tilde{T}^{(1)} &= PT^{(1)}P, \\
\tilde{T}^{(2)} &= PT^{(2)}P - PT^{(1)}PT^{(1)}S - PT^{(1)}ST^{(1)}P - ST^{(1)}PT^{(1)}P, \\
(60) \quad \tilde{T}^{(k)} &= -\sum_{p=1}^{k} (-1)^p \sum_{\substack{\nu_1+\dots+\nu_p=k \\ k_1+\dots+k_{p+1}=p-1 \\ \nu_j \ge 1, k_j \ge 0}} S^{(k_1)}T^{(\nu_1)}S^{(k_2)}\dots S^{(k_p)}T^{(\nu_p)}S^{(k_{p+1})}.
\end{align*}
$$

REMARK 7.2. Note that the matrices $\tilde{T}^{(k)}$ are selfadjoint, which follows from the fact that for real $\varepsilon$ the matrices $T(\varepsilon)$ and $P(\varepsilon)$ (and hence $D(\varepsilon)$) are selfadjoint.

Next, we define the matrix $\tilde{T}(\varepsilon)$ as

$$ \tilde{T}(\varepsilon) = \frac{1}{\varepsilon} D(\varepsilon) = \sum_{k=0}^{\infty} \varepsilon^k \tilde{T}^{(k+1)}, $$

such that $\tilde{T}(0) = \tilde{T}^{(1)}$. Note that by Remark 7.2, all the matrices $\tilde{T}^{(k+1)}$ are Hermitian and all the eigenvalues of $\tilde{T}^{(1)}(\varepsilon)$ are holomorphic functions of $\varepsilon$. The idea of the reduction process is to apply the perturbation theory to the matrix $\tilde{T}^{(1)}(\varepsilon)$.

Let $\tilde{T}^{(1)}$ have $s$ eigenvalues $\tilde{\lambda}_1, \dots, \tilde{\lambda}_s$ with multiplicities

$$ \ell_1 + \dots + \ell_s = m, $$

where we take into account only the $m$ eigenvalues<fnref n="9" /> in the subspace spanned by $P(\varepsilon)$.

Then $\tilde{\lambda}_1, \dots, \tilde{\lambda}_s$ determine the splitting of $\lambda$ in the following way.

LEMMA 7.3 (A summary of [22, Ch. II, §2.3]). Let

$$ \underbrace{\tilde{\lambda}_{1,1}(\varepsilon), \dots, \tilde{\lambda}_{1,\ell_1}(\varepsilon), \dots, \tilde{\lambda}_{s,1}(\varepsilon), \dots, \tilde{\lambda}_{s,\ell_s}(\varepsilon)}_{\text{pert. of } \tilde{\lambda}_1} $$

be the holomorphic functions for the perturbations of the eigenvalues $\tilde{\lambda}_1, \dots, \tilde{\lambda}_s$ of $\tilde{T}(\varepsilon)$. Then the holomorphic functions corresponding to perturbations of the eigenvalue $\lambda$ of the original matrix $T(\varepsilon)$ are given by

$$ \{\lambda_1(\varepsilon), \dots, \lambda_m(\varepsilon)\} = \bigcup_{k=1}^{s} \{\lambda + \varepsilon \tilde{\lambda}_{k,1}(\varepsilon), \dots, \lambda + \varepsilon \tilde{\lambda}_{k,\ell_k}(\varepsilon)\}. $$

Moreover, the expansions $\tilde{P}_{k,j}(\varepsilon)$ of the projectors on the eigenspaces of $\tilde{T}(\varepsilon)$ (corresponding to $\tilde{\lambda}_{k,j}(\varepsilon)$) give the expansions of the projectors on the eigenspaces of $T(\varepsilon)$ corresponding to $\lambda_1(\varepsilon), \dots, \lambda_m(\varepsilon)$.

Lemma 7.3 can be applied recursively: for each individual eigenvalue $\tilde{\lambda}_k$ (of multiplicity $\ell_k > 1$) we can consider the corresponding reduced matrix

$$ \tilde{T}^{(2)}(\varepsilon) = \frac{1}{\varepsilon}(\tilde{T}^{(1)}(\varepsilon) - \tilde{\lambda}_k I)\tilde{P}_k(\varepsilon), $$

<footnote n="9">⁹The other $n-m$ eigenvalues are 0.</footnote>
$$S_0 = \left\{ \boxed{\phantom{\substack{\text{light grey} \\ \text{dark grey}}} \phantom{\substack{\text{light grey} \\ \text{light grey}}} \phantom{\substack{\text{light grey} \\ \text{light grey}} \phantom{\substack{\text{light grey} \\ \text{light grey}} } } \right\}, S_1 = \left\{ \boxed{\phantom{\substack{\text{light grey} \\ \text{light grey}}} \phantom{\substack{\text{light grey} \\ \text{light grey}}} \phantom{\substack{\text{light grey} \\ \text{light grey}}} \phantom{\substack{\text{light grey} \\ \text{light grey}} \phantom{\substack{\text{light grey} \\ \text{light grey}} } } \right\}, S_2 = \left\{ \boxed{\phantom{\substack{\text{light grey} \\ \text{light grey}}} \phantom{\substack{\text{light grey} \\ \text{light grey}}} \phantom{\substack{\text{light grey} \\ \text{light grey}}} \phantom{\substack{\text{light grey} \\ \text{light grey}}} \phantom{\substack{\text{light grey} \\ \text{light grey}} \phantom{\substack{\text{light grey} \\ \text{light grey}} } } \right\}, S_3 = \left\{ \boxed{\phantom{\substack{\text{light grey} \\ \text{light grey}}} \phantom{\substack{\text{light grey} \\ \text{light grey}}} \phantom{\substack{\text{light grey} \\ \text{light grey}}} \phantom{\substack{\text{light grey} \\ \text{light grey}}} \phantom{\substack{\text{light grey} \\ \text{light grey}} \phantom{\substack{\text{light grey} \\ \text{light grey}} } } \right\}, S_4 = \left\{ \boxed{\phantom{\substack{\text{light grey} \\ \text{light grey}}} \phantom{\substack{\text{light grey} \\ \text{light grey}}} \phantom{\substack{\text{light grey} \\ \text{light grey}}} \phantom{\substack{\text{light grey} \\ \text{light grey}}} \phantom{\substack{\text{light grey} \\ \text{light grey}} \phantom{\substack{\text{light grey} \\ \text{light grey}} } } \right\}.$$

<caption>FIG. 8.1. Classes of staircase matrices. White color stands for zero blocks.</caption>

where $\tilde{P}_k(\varepsilon)$ is the perturbation of the total projection on the $\ell_k$-dimensional eigenspace corresponding to $\tilde{\lambda}_k$ (which can be computed as in the previous subsection). Depending on the eigenvalues of the main term of the reduced matrix, either the splitting will occur again, or there will be no splitting; in any case, after a finite number of steps, all the individual eigenvalues will be split into simple (multiplicity 1) eigenvalues$^{10}$.

**8. Results on eigenvectors for finitely smooth kernels.** In this section we are going to prove the following result for the multivariate case.

**THEOREM 8.1.** Let the radial kernel be as in Theorem 6.3, with $K_\varepsilon$ positive semi-definite on $[0, \varepsilon_0]$, $\varepsilon_0 > 0$ and analytic in $\varepsilon$. If $\det W_{\le r-1} \ne 0$, $V_{\le r-1}$ is full rank, and $Q_\perp^T D_{(2r-1)} Q_\perp$ is invertible, then the eigenvectors corresponding to the last group of $O(\varepsilon^{2r-1})$ eigenvalues, are given by the columns of

$$Q_\perp U,$$

where the columns of $U$ are eigenvectors of $Q_\perp^T D_{(2r-1)} Q_\perp$, and the matrix $Q_\perp$ comes from the full QR factorization (22) of $V_{\le r-1}$ (as in Theorem 6.3).

We conjecture that for finitely-smooth kernels as well, the individual eigenvectors for the $O(\varepsilon^{2s})$ groups, $0 \le s \le r-1$, can be obtained as in Conjecture 1.

**8.1. Block staircase matrices.** We first need some facts about a class of so-called block staircase matrices. Let $n = (n_0, \dots, n_s) \in \mathbb{N}^{s+1}$ such that

$$n = n_0 + \dots + n_s$$

and consider the following block partition of a matrix $A \in \mathbb{C}^{m \times m}$

$$A = \begin{bmatrix} A^{(0,0)} & \dots & A^{(0,s)} \\ \vdots & \ddots & \vdots \\ A^{(s,0)} & \dots & A^{(s,s)} \end{bmatrix}$$

where the blocks are of size $A^{(k,l)} \in \mathbb{C}^{n \times n}$. Assuming that the partition is fixed, we define the classes $S_p$ of “staircase” matrices

$$\{0\} = S_{-1} \subset S_0 \subset S_1 \subset S_2 \subset \dots \subset S_{2s-1} = S_{2s} = S_{2s+1} = \dots = \mathbb{C}^{n \times n},$$

such that the matrices in $S_p$ have nonzero blocks only up to the $s$-th antidiagonal

$$S_p = \{ A \in \mathbb{C}^{m \times n} \mid A^{(k,l)} = 0 \text{ for all } k+l > p \}.$$

In Fig. 8.1, we illustrate the classes for $s=2$:

The following obvious property will be useful.

**LEMMA 8.2.** 1. $A \in S_p, B \in S_q \Rightarrow AB \in S_{p+q}$.

2. For any $A \in S_p$ and upper triangular $R$, it follows that $RAR^T \in S_p$.

3. For any $A \in S_p$ and block-diagonal matrix $\Lambda$ it holds that $\Lambda A, A\Lambda \in S_p$.

*Proof.* The proof follows from straightforward verification. □

<footnote>$^{10}$This follows from our assumption that the eigenvalues are simple generically.</footnote>
## 8.2. Proof of Theorem 8.1.

*Proof.* The kernel matrix has an expansion with only even powers of $\varepsilon$ until $2r-2$

$$ K_{\varepsilon} = \sum_{j=0}^{r-1} \varepsilon^{2j} \mathbf{V}_{\le 2j} \mathbf{W}_{/2j} \mathbf{V}_{\le 2j}^{\top} + \varepsilon^{2r-1} f_{2r-1} \mathbf{D}_{(2r-1)} + \mathcal{O}(\varepsilon^{2r}), $$

since odd block diagonals $\mathbf{W}_{/2j-1}$, $j < r$ vanish. We look at the transformed matrix

$$ \mathbf{T}(\varepsilon) = \mathbf{Q}_{\text{full}}^{\top} \mathbf{K}_{\varepsilon} \mathbf{Q}_{\text{full}} = \mathbf{T}_0 + \varepsilon \mathbf{T}_1 + \cdots + \varepsilon^{2r-1} \mathbf{T}_{2r-1} + \mathcal{O}(\varepsilon^{2r}), $$

where $\mathbf{Q}_{\text{full}}$ is the matrix of the full QR decomposition of $\mathbf{V}_{\le 2r-2}$

$$ \mathbf{T}_s = \begin{cases} \mathbf{Q}_{\text{full}}^{\top} \mathbf{V}_{\le s} \mathbf{W}_{/s} \mathbf{V}_{\le s}^{\top} \mathbf{Q}_{\text{full}}, & s \text{ even, } s < 2r-1, \\ 0, & s \text{ odd, } s < 2r-1, \\ f_{2r-1} \mathbf{Q}_{\text{full}}^{\top} \mathbf{D}_{(2r-1)} \mathbf{Q}_{\text{full}}, & s = 2r-1. \end{cases} $$

Due to the fact that that $\mathbf{W}_{/s}$ is block antidiagonal, $\mathbf{Q}_{\text{full}}^{\top} \mathbf{V}_{\le s}$ is upper triangular and by Lemma 8.2, we have that $\mathbf{T}_s$ is block staircase for $s < 2r-1$, i.e. $\mathbf{T}_{/s} \in S_s$.

We proceed by a series of Kato's reductions, according to Lemma 7.3. At each order of $\varepsilon$, a multiple eigenvalue 0 is split into a group of nonzero eigenvalues and eigenvalue 0 of smaller multiplicity. Formally, we consider a sequence of reduced matrices

$$ \mathbf{T}^{s+1} \tilde{\mathbf{T}}(\varepsilon) = \frac{1}{\varepsilon} \mathbf{T}^{s} \tilde{\mathbf{T}}(\varepsilon) \mathbf{P}_s(\varepsilon), \quad \mathbf{T}^{0} \tilde{\mathbf{T}}(\varepsilon) = \mathbf{T}(\varepsilon) $$

where $\mathbf{P}_k(\varepsilon)$ is the projector onto the perturbation of the nullspace of $\mathbf{T}^{s} \tilde{\mathbf{T}}(0)$ (i.e., the invariant subspace associated with the eigenvalue 0). Its power series expansion

$$ \mathbf{T}^{s} \tilde{\mathbf{T}}(\varepsilon) = \sum_{j=s}^{\infty} \varepsilon^{j-s} \cdot \mathbf{T}_{j}^{s} $$

can be computed according to (60). For each matrix we will be interested only in the first $2r-s$ terms, as summarized below,

$$ \begin{array}{c|ccccc} 1 & \varepsilon & \varepsilon^2 & \cdots & \varepsilon^{2r-1} \\ \hline \mathbf{T}_0 & & & & \\ \mathbf{T}_1 & \tilde{\mathbf{T}}_1 & & & \\ \mathbf{T}_2 & \tilde{\mathbf{T}}_2 & \mathbf{T}_2^{2\tilde{\mathbf{T}}_2} & & \\ \vdots & & & \ddots & \\ \mathbf{T}_{2r-1} & \tilde{\mathbf{T}}_{2r-1} & \mathbf{T}_{2r-1}^{2\tilde{\mathbf{T}}_{2r-1}} & \mathbf{T}_{2r-1}^{2r-1\tilde{\mathbf{T}}_{2r-1}} & \end{array} $$

since we are interested only in the terms

$$ \mathbf{T}_0, \mathbf{T}_1, \mathbf{T}_2, \dots, \mathbf{T}_{2k}, $$

whose eigenvectors give limiting eigenvectors for the original matrix $\mathbf{T}(\varepsilon)$.

Next, we will look in detail at the form of coefficients of the reduced matrices. The projector on the image space of $\mathbf{T}_0$ is $\mathbf{\Pi}_0 = \mathbf{e}_0 \mathbf{e}_0^{\top}$ (where $\mathbf{e}_0 = (1, 0, \dots, 0)^{\top}$), the projector on the nullspace is

$$ \mathbf{P}_0 = \begin{bmatrix} 0 & \\ & \mathbf{I}_{n-1} \end{bmatrix}, $$
and the matrix $S = S_0$ in (59). By examining the terms in (59), we have that the coefficients of the reduced matrices preserve the staircase class, i.e. $\tilde{T}_s \in S_s$. This can be seen by verifying that if $A_s \in S_s$ and $S^{(j)}$ are diagonal, then the products

$$ (61) \qquad S^{(k_1)} A^{(\nu_1)} S^{(k_2)} \cdots S^{(k_p)} A^{(\nu_p)} S^{(k_{p+1})} \in S_s $$

if $\nu_1 + \cdots + \nu_p = s$. Next, we note that since $T_1 \in S_1$, then we have

$$ \tilde{T}_1 = P_0 T_1 P_0 = 0. $$

Hence we have that $P_1 = I$ and $S_1 = 0$, and the second step of reduction does not change the matrices, i.e. $\tilde{T}_s = \tilde{T}_s$.

Proceeding by induction, at the step of reduction $(2s-2) \to (2s-1)$ the staircase order of the matrices is preserved due to (61), and block diagonality of $P_{2s-1}$ and $\tilde{T}_{2s-2}$. Since the reduction step $(2s-1) \to 2s$ does not change anything, we get

$$ (62) \qquad P_{2s} = \left[ \begin{array}{c} 0(\#P_s) \times (\#P_s) \\ I_{(n-\#P_s) \times (n-\#P_s)} \end{array} \right], $$

which we know from Theorem 6.3, and

$$ (63) \qquad \tilde{T}_{2s} = \left[ \begin{array}{c} 0(\#P_s) \times (\#P_s) \\ *^{(\#H_{s+1})} \times (\#H_{s+1}) \\ 0^{(n-\#P_{s+1})} \times (n-\#P_{s+1}) \end{array} \right], $$

The last reduction step is different, as we get $\tilde{T}_{2r-1}$ which is not equal to zero. In order to obtain $\tilde{T}_{2r-1}$, we note the following: at the first step of the reduction the matrices $\tilde{T}_{2j-1}$ defined by (60), are the sums of the terms running over multi-indices

$$ \nu_1 + \cdots + \nu_p = 2j-1, $$

where at least one of $\nu_\ell$ should be odd and all $\nu_\ell \le 2j-1$. Therefore, we have $\tilde{T}_{2j-1} = 0$ if $j < r$ and $\tilde{T}_{2r-1} = PT_{2r-1}P$. Proceeding by induction, we get that

$$ \tilde{T}_{2j-1} = \begin{cases} 0, & j < r \text{ and } s < 2r-1, \\ f_{2r-1} P_{2r-2} T_{2r-1} P_{2r-2}, & j = r \text{ and } s = 2r-1, \end{cases} $$

where $P_{2r-2}$ is defined in (62).

Thus we have that the limiting eigenvectors of $Q_{\text{full}}T_\epsilon Q_{\text{full}}^\text{T}$ for the order $O(\epsilon^{2r-1})$ are the limiting eigenvectors (corresponding to non-zero eigenvalues) of the matrix

$$ Q_{\perp} Q_{\perp}^{\text{T}} D_{(2r-1)} Q_{\perp} Q_{\perp}^{\text{T}}, $$

where $Q_{\text{full}} = [Q \ Q_\perp]$ is the splitting of $Q_{\text{full}}$ such that $Q_\perp \in \mathbb{R}^{n \times (n-\#P_{r-1})}$. $\square$

**9. Discussion.** We have shown that kernel matrices become tractable in the flat limit, and exhibit deep ties to orthogonal polynomials. We would like to add some remarks and highlight some open problems.

First, we expect our analysis to generalise in a mostly straightforward manner to the “continuous” case, i.e., to kernel integral operators. This should make it possible
to examine a double asymptotic, in which $n \to \infty$ as $\varepsilon \to 0$. One could then leverage recent results on the asymptotics of orthogonal polynomials, for instance [23].

Second, our results may be used empirically to create preconditioners for kernel matrices. There is already a vast literature on approximate kernel methods, including in the flat limit (e.g., [13, 25]), and future work should examine how effective polynomial preconditioners are compared to other available methods.

Third, many interesting problems (e.g., spectral clustering [39]) involve not the kernel matrix itself but some rescaled variant. We expect that multiplicative perturbation theory could be brought to bear here [36].

Finally, while multivariate polynomials are relatively well-understood objects, our analysis also shows that in the finite smoothness case, a central role is played by a different class of objects: namely, multivariate polynomials are replaced by the eigenvectors of distance matrices of an odd power. To our knowledge, very little has been proved about such objects but some literature from statistical physics [6] points to a link to “Anderson localization”. Anderson localization is a well-known phenomenon in physics whereby eigenvectors of certain operators are localised, in the sense of having fast decay over space. This typically does not hold for orthogonal polynomials, which tend rather to be localised in frequency. Thus, we conjecture that eigenvectors of completely smooth kernels are localised in frequency, contrary to eigenvectors of finitely smooth kernels, which (at low energies) are localised in space. The results in [6] are enough to show that this holds for the exponential kernel in $d = 1$, but extending this to $d > 1$ and higher regularity orders is a fascinating and probably nontrivial problem.

## Appendix A. Proofs for elementary symmetric polynomials.

*Proof of Lemma 3.4.*

$$ 
\begin{align} 
e_s(\varepsilon) &= \sum_{1 \le t_1 < \dots < t_s \le n} \prod_{j=1}^{s} \varepsilon^{L_{t_j}} (\tilde{\lambda}_{t_j} + \mathcal{O}(\varepsilon)) \tag{64} \\ 
&= \varepsilon^{L_1 + \dots + L_s} \underbrace{\left( \sum_{1 \le t_1 < \dots < t_s \le n} \prod_{j=1}^{s} \tilde{\lambda}_{t_j} + \mathcal{O}(\varepsilon) \right)}_{\tilde{e}_s} \nonumber 
\end{align} 
$$

which follows from the fact that $L_{t_1} + \dots + L_{t_s}$ is minimized at $t_j = j$.

2. The case $s = n$ is obvious, because there is only one possible tuple $(t_1, \dots, t_s)$. Consider the case $L_s < L_{s+1}$, $1 \le s < n$. If $t_s > s$, then the sum is increased

$$L_1 + \dots + L_s < L_{t_1} + \dots + L_{t_s},$$

hence $(1, \dots, s)$ is the only possible choice for $(t_1, \dots, t_s)$ in (64). □

*Proof of Lemma 3.5.* When (18) is satisfied, we need to find

$$ \tilde{e}_{s+k} = \sum_{\substack{1 \le t_1 < \dots < t_{s+k} \le n \\ L_1 + \dots + L_{s+k} = L_{t_1} + \dots + L_{t_{s+k}}}} \prod_{j=1}^{s+k} \tilde{\lambda}_{t_j} $$

where, as in the previous case, the minimum sum $L_{t_1} + \dots + L_{t_{s+k}}$ is achieved by

$$L_1 + \dots + L_s + L_{s+1} + \dots + L_{s+k},$$
and the sum increases if $t_s > s$ or if $t_{s+k} > s+m$. Therefore, $(t_1, \dots, t_s) = (1, \dots, s)$ and the main term becomes

$$ \tilde{e}_{s+k} = \prod_{i=1}^{s} \tilde{\lambda}_i \left( \sum_{s+1 \le t_{s+1} < \dots < t_{s+k} \le s+m} \prod_{j=1}^{k} \tilde{\lambda}_{t_{s+j}} \right). \quad \square $$

## Appendix B. Proofs for saddle point matrices.

*Proof of Lemma 3.11.* We note that $(\det \mathbf{R})^2 = \det(\mathbf{V}^\mathrm{T}\mathbf{V})$ and

$$ \begin{aligned} \det \begin{bmatrix} \mathbf{A} & \mathbf{V} \\ \mathbf{V}^\mathrm{T} & 0 \end{bmatrix} &= \det \begin{pmatrix} [\mathbf{Q}_{\mathrm{full}}^\mathrm{T}] & \mathbf{I}_r \\ [\mathbf{V}^\mathrm{T}] & 0 \end{pmatrix} \begin{bmatrix} \mathbf{A} & \mathbf{V} \\ \mathbf{V}^\mathrm{T} & 0 \end{bmatrix} [\mathbf{Q}_{\mathrm{full}}] \\ &= \det \begin{bmatrix} \mathbf{Q}_{\mathrm{full}}^\mathrm{T} \mathbf{A} \mathbf{Q}_{\mathrm{full}}^\mathrm{T} & \mathbf{Q}_{\mathrm{full}}^\mathrm{T} \mathbf{A} \mathbf{Q}_{\mathrm{full}}^\mathrm{T} \\ \mathbf{Q}_{\mathrm{full}}^\mathrm{T} \mathbf{A} \mathbf{Q}_{\mathrm{full}}^\mathrm{T} & \mathbf{Q}_{\mathrm{full}}^\mathrm{T} \mathbf{A} \mathbf{Q}_{\mathrm{full}}^\mathrm{T} \end{bmatrix} = (-1)^r (\det \mathbf{R})^2 \det(\mathbf{Q}_{\mathrm{full}}^\mathrm{T} \mathbf{A} \mathbf{Q}_{\mathrm{full}}^\mathrm{T}). \quad \square \end{aligned} $$

Before proving Lemma 3.12, we need a technical lemma first.

**LEMMA B.1.** For $\mathbf{G} \in \mathbb{R}^{r \times r}$ and $\mathbf{M} = \begin{pmatrix} \mathbf{A} & \mathbf{B} \\ \mathbf{C} & \mathbf{D} \end{pmatrix} \in \mathbb{R}^{n \times n}$, $\mathbf{A} \in \mathbb{R}^{r \times r}$ it holds that

$$ [t^r] \{ \det \begin{pmatrix} \mathbf{M} + t \begin{pmatrix} \mathbf{G} & 0 \\ 0 & 0 \end{pmatrix} \} = \det \mathbf{G} \det \mathbf{D}. \quad \square $$

*Proof.* From [5, Theorem 1], we have that<fnref n="11" />

*Proof of Lemma 3.12.* Due to invariance under similarity transformations of the elementary polynomials, we get

$$ \begin{aligned} [t^r] \{e_k(\mathbf{A} + t\mathbf{V}\mathbf{V}^\mathrm{T})\} &= [t^r] \{e_k(\mathbf{Q}_{\mathrm{full}}^\mathrm{T}(\mathbf{A} + t\mathbf{V}\mathbf{V}^\mathrm{T})\mathbf{Q}_{\mathrm{full}})\} \\ &= [t^r] \left\{ e_k \left( \underbrace{\begin{bmatrix} \mathbf{Q}^\mathrm{T}\mathbf{A}\mathbf{Q} + t\mathbf{R}\mathbf{R}^\mathrm{T} & \mathbf{Q}^\mathrm{T}\mathbf{A}\mathbf{Q}_\perp \\ \mathbf{Q}_\perp^\mathrm{T}\mathbf{A}\mathbf{Q} & \mathbf{Q}_\perp^\mathrm{T}\mathbf{A}\mathbf{Q}_\perp \end{bmatrix}}_{\mathcal{B}(t)} \right) \right\} \\ &= \sum_{\substack{|\mathcal{J}|=k \\ \mathcal{J} \subseteq \{1, \dots, n\}}} [t^r] \{\det(\mathbf{B}_{\mathcal{J},\mathcal{J}}(t))\} = \sum_{\substack{|\mathcal{J}|=k \\ \{1, \dots, r\} \subseteq \mathcal{J} \subseteq \{1, \dots, n\}}} [t^r] \{\det(\mathbf{B}_{\mathcal{J},\mathcal{J}}(t))\}, \end{aligned} $$

where the last equality follows from the fact that the polynomial $\det(\mathbf{B}_{\mathcal{J},\mathcal{J}}(t))$ has degree that is equal to the cardinality of the intersection $\mathcal{J} \cap \{1, \dots, r\}$. Any such $\mathcal{J}$ can be written as $\{1, \dots, r\} \cup \mathcal{J}'$, where $\mathcal{J}' \subseteq \{r+1, \dots, n\}$. Applying Lemma B.1 to each term individually, we get

$$ [t^r] \{\det(\mathbf{B}_{\mathcal{J},\mathcal{J}}(t))\} = \det(\mathbf{R}\mathbf{R}^\mathrm{T}) \det((\mathbf{Q}_\perp^\mathrm{T}\mathbf{A}\mathbf{Q}_\perp)^{\mathcal{J}',\mathcal{J}'})$$

<footnote n="11">¹¹Note that in the sum in [5, eq. (7)] only one term is nonzero.</footnote>
thus, summing over all $\mathcal{J}' \subseteq \{r+1, \dots, n\}$ yields

$$ \det(\mathbf{RR}^\top) \sum_{\substack{\mathcal{J}' \subseteq \{r+1, \dots, n\} \\ |\mathcal{J}'|=k-r}} \det((\mathbf{Q}_\perp^\top \mathbf{A}\mathbf{Q}_\perp)_{\mathcal{J}', \mathcal{J}'}) = \det(\mathbf{RR}^\top)e_{k-r}(\mathbf{Q}_\perp^\top \mathbf{A}\mathbf{Q}_\perp). \quad \square $$

**Acknowledgments.** The authors would like to acknowledge Pierre Comon for his help at the beginning of this project, and Cosme Louart and Malik Tiomoko for their help in improving the presentation of the results. We would also like to thank two anonymous referees for their careful reading of our work and many suggestions.

REFERENCES

[1] OEIS orderings wiki page. https://oeis.org/wiki/Orderings. Accessed: 2019-09-25.

[2] B. BAXTER, Norm estimates for inverses of Toeplitz distance matrices, Journal of Approximation Theory, 79 (1994), pp. 222–242.

[3] B. BAXTER, Preconditioned conjugate gradients, radial basis functions, and Toeplitz matrices, Computers & Mathematics with Applications, 43 (2002), pp. 305–318.

[4] M. BENZI AND V. SIMONCINI, On the eigenvalues of a class of saddle point matrices, Numerische Mathematik, 103 (2006), pp. 173–196.

[5] R. BHATIA AND T. JAIN, Higher order derivatives and perturbation bounds for determinants, Linear Algebra and its Applications, 431 (2009), pp. 2102–2108, https://doi.org/https://doi.org/10.1016/j.laa.2009.07.003.

[6] E. BOGOMOLNY, O. BOHIGAS, AND C. SCHMIT, Spectral properties of distance matrices, Journal of Physics A: Mathematical and General, 36 (2003), p. 3595.

[7] R. COUILLET AND F. BENAYCH-GEORGES, Kernel spectral clustering of large dimensional data, Electronic Journal of Statistics, 10 (2016), pp. 1393–1454.

[8] D. COX, J. LITTLE, AND D. O'SHEA, Ideals, Varieties and Algorithms: An Introduction to Computational Algebraic Geometry and Commutative Algebra, Springer, 2nd ed., 1997.

[9] S. DE MARCHI AND R. SCHABACK, Nonstandard kernels and their applications, Dolomites Research Notes on Approximation, 2 (2009).

[10] T. A. DRISCOLL AND B. FORNBERG, Interpolation in the limit of increasingly flat radial basis functions, Computers & Mathematics with Applications, 43 (2002), pp. 413–422.

[11] N. EL KAROUI, The spectrum of kernel random matrices, The Annals of Statistics, 38 (2010), pp. 1–50.

[12] G. E. FASSHAUER, Meshfree approximation methods with MATLAB, vol. 6, World Scientific, 2007.

[13] B. FORNBERG, E. LARSSON, AND N. FLYER, Stable computations with Gaussian radial basis functions, SIAM Journal on Scientific Computing, 33 (2011), pp. 869–892.

[14] B. FORNBERG, G. WRIGHT, AND E. LARSSON, Some observations regarding interpolants in the limit of flat radial basis functions, Computers & Mathematics with Applications, 47 (2004), pp. 37–55.

[15] B. FORNBERG AND J. ZUEV, The runge phenomenon and spatially variable shape parameters in rbf interpolation, Computers & Mathematics with Applications, 54 (2007), pp. 379–398.

[16] D. GARREAU, W. JITKRITTUM, AND M. KANAGAWA, Large sample analysis of the median heuris- tic, arXiv preprint arXiv:1707.07269, (2017).

[17] M. GASCA AND T. SAUER, Polynomial interpolation in several variables, Advances in Compu- tational Mathematics, 12 (2000), p. 377.

[18] U. GRENANDER AND G. SZEGÖ, Toeplitz Forms and their Applications, Chelsea, New York, 1984.

[19] R. A. HORN AND C. R. JOHNSON, Matrix Analysis, Cambridge University Press, 1990.

[20] I. IPSEN AND R. REHMAN, Perturbation bounds for determinants and characteristic polynomials, SIAM Journal on Matrix Analysis and Applications, 30 (2008), pp. 762–776, https://doi.org/10.1137/070704770.

[21] T. JAIN, Derivatives for antisymmetric tensor powers and perturbation bounds, Linear Algebra and its Applications, 435 (2011), pp. 1111–1121, https://doi.org/https://doi.org/10.1016/j.laa.2011.02.026.

[22] T. KATO, Perturbation theory for linear operators, Springer-Verlag, 2nd corrected ed., 1995.

[23] A. KROÓ AND D. LUBINSKY, Christoffel functions and universality in the bulk for multivariate orthogonal polynomials, Canadian Journal of Mathematics, 65 (2013), pp. 600–620.
[24] E. LARSSON AND B. FORNBERG, *Theoretical and computational aspects of multivariate interpolation with increasingly flat radial basis functions*, Computers & Mathematics with Applications, 49 (2005), pp. 103–130.

[25] S. LE BORNE, *Factorization, symmetrization, and truncated transformation of radial basis function-ga stabilized gaussian radial basis functions*, SIAM Journal on Matrix Analysis and Applications, 40 (2019), pp. 517–541.

[26] Y. J. LEE AND C. A. MICCHELLI, *On collocation matrices for interpolation and approximation*, Journal of Approximation Theory, 174 (2013), pp. 148–181, https://doi.org/https://doi.org/10.1016/j.jat.2013.06.007.

[27] Y. J. LEE, C. A. MICCHELLI, AND J. YOON, *A study on multivariate interpolation by increasingly flat kernel functions*, Journal of Mathematical Analysis and Applications, 427 (2015), pp. 74–87.

[28] C. A. MICCHELLI, *Interpolation of scattered data: Distance matrices and conditionally positive definite functions*, Constructive Approximation, 2 (1986), pp. 11–22.

[29] M. MIRANDA AND P. TILLI, *Asymptotic spectra of Hermitian block Toeplitz matrices and preconditioning results*, SIAM Journal on Matrix Analysis and Applications, 21 (2000), pp. 867–881.

[30] T. SAUER, *Polynomial interpolation in several variables: lattices, differences, and ideals*, 12 (2006), pp. 191–230.

[31] R. SCHABACK, *Multivariate interpolation by polynomials and radial basis functions*, Constructive Approximation, 21 (2005), pp. 293–317.

[32] R. SCHABACK, *Limit problems for interpolation by analytic radial basis functions*, Journal of Computational and Applied Mathematics, 212 (2008), pp. 127–149.

[33] R. SCHABACK AND H. WENDLAND, *Kernel techniques: from machine learning to meshless methods*, Acta Numerica, 15 (2006), pp. 543–639.

[34] B. SCHÖLKOPF AND A. J. SMOLA, *Learning with kernels: support vector machines, regularization, optimization, and beyond*, MIT press, 2002.

[35] G. SONG, J. RIDDLE, G. E. FASSHAUER, AND F. J. HICKERNELL, *Multivariate interpolation with increasingly flat radial basis functions of finite smoothness*, Advances in Computational Mathematics, 36 (2012), pp. 485–501.

[36] F. SOSA AND J. MORO, *First order asymptotic expansions for eigenvalues of multiplicatively perturbed matrices*, SIAM Journal on Matrix Analysis and Applications, 37 (2016), pp. 1478–1504.

[37] M. L. STEIN, *Interpolation of Spatial Data: Some Theory for Kriging*, Springer, 1999.

[38] E. E. TYRTYSHNIKOV, *A unifying approach to some old and new theorems on distribution and clustering*, Linear Algebra and its Applications, 232 (1996), pp. 1–43, https://doi.org/10.1016/0024-3795(94)00025-5.

[39] U. VON LUXBURG, M. BELKIN, AND O. BOUSQUET, *Consistency of spectral clustering*, The Annals of Statistics, (2008), pp. 555–586.

[40] R. WANG, Y. LI, AND E. DARVE, *On the numerical rank of radial basis function kernels in high dimensions*, SIAM Journal on Matrix Analysis and Applications, 39 (2018), pp. 1810–1835.

[41] A. J. WATHEN, *Realistic eigenvalue bounds for Galerkin mass matrix*, IMA Journal on Numerical Analysis, 7 (1987), pp. 449–457.

[42] A. J. WATHEN AND S. ZHU, *On spectral distribution of kernel matrices related to radial basis functions*, Numerical Algorithms, 70 (2015), pp. 709–726, https://doi.org/10.1007/s11075-015-9970-0.

[43] H. WENDLAND, *Scattered Data Approximation*, Cambridge Monographs on Applied and Computational Mathematics, Cambridge University Press, 2004, https://doi.org/10.1017/CBO9780511617539.

[44] C. K. WILLIAMS AND C. E. RASMUSSEN, *Gaussian processes for machine learning*, vol. 2, MIT press Cambridge, MA, 2006.

[45] V. A. ZORICH, *Mathematical Analysis I*, Springer, 2004.