# Local optimisation of Nyström samples
through stochastic gradient descent

Matthew<EMAIL_ADDRESS>444Cardiff University, School of
Mathematics
Abacws, Senghennydd Road, Cardiff, CF24 4AG, United Kingdom Bertrand
<EMAIL_ADDRESS>44footnotemark: 4

###### Abstract

We study a relaxed version of the column-sampling problem for the Nyström
approximation of kernel matrices, where approximations are defined from
multisets of landmark points in the ambient space; such multisets are referred
to as Nyström samples. We consider an unweighted variation of the radial
squared-kernel discrepancy (SKD) criterion as a surrogate for the classical
criteria used to assess the Nyström approximation accuracy; in this setting,
we discuss how Nyström samples can be efficiently optimised through stochastic
gradient descent. We perform numerical experiments which demonstrate that the
local minimisation of the radial SKD yields Nyström samples with improved
Nyström approximation accuracy.

Keywords: Low-rank matrix approximation; Nyström method; reproducing kernel
Hilbert spaces;
stochastic gradient descent.

##  1 Introduction

In Data Science, the Nyström method refers to a specific technique for the
low-rank approximation of symmetric positive-semidefinite (SPSD) matrices; see
e.g. [5, 11, 18, 10, 4]. Given an $N\times N$ SPSD matrix $\mathbf{K}$, with
$N\in\mathbb{N}$, the Nyström method consists of selecting a sample of
$n\in\mathbb{N}$ columns of $\mathbf{K}$, generally with $n\ll N$, and next
defining a low-rank approximation $\hat{\mathbf{K}}$ of $\mathbf{K}$ based on
this sample of columns. More precisely, let
$\boldsymbol{c}_{1},\cdots,\boldsymbol{c}_{N}\in\mathbb{R}^{N}$ be the columns
of $\mathbf{K}$, so that
$\mathbf{K}=(\boldsymbol{c}_{1}|\cdots|\boldsymbol{c}_{N})$, and let
$I=\\{i_{1},\cdots,i_{n}\\}\subseteq\\{1,\cdots,N\\}$ denote the indices of a
sample of $n$ columns of $\mathbf{K}$ (note that $I$ is a multiset, i.e. the
indices of some columns might potentially be repeated). Let
$\mathbf{C}=(\boldsymbol{c}_{i_{1}}|\cdots|\boldsymbol{c}_{i_{n}})$ be the
$N\times n$ matrix defined from the considered sample of columns of
$\mathbf{K}$, and let $\mathbf{W}$ be the $n\times n$ principal submatrix of
$\mathbf{K}$ defined by the indices in $I$, i.e. the $k,l$ entry of
$\mathbf{W}$ is $[\mathbf{K}]_{i_{k},i_{l}}$, the $i_{k},i_{l}$ entry of
$\mathbf{K}$. The Nyström approximation of $\mathbf{K}$ defined from the
sample of columns indexed by $I$ is given by

$\hat{\mathbf{K}}=\mathbf{C}\mathbf{W}^{{\dagger}}\mathbf{C}^{T},$ (1)

with $\mathbf{W}^{{\dagger}}$ the Moore-Penrose pseudoinverse of $\mathbf{W}$.
The column-sampling problem for Nyström approximation consists of designing
samples of columns such that the induced approximations are as accurate as
possible (see Section 1.2 for more details).

###  1.1 Kernel Matrix Approximation

If the initial SPSD matrix $\mathbf{K}$ is a kernel matrix, defined from a
SPSD kernel $K$ and a set or multiset of points
$\mathcal{D}=\\{x_{1},\cdots,x_{N}\\}\subseteq\mathscr{X}$ (and with
$\mathscr{X}$ a general ambient space), i.e. the $i,j$ entry of $\mathbf{K}$
is $K(x_{i},x_{j})$, then a sample of columns of $\mathbf{K}$ is naturally
associated with a subset of $\mathcal{D}$; more precisely, a sample of columns
$\\{\boldsymbol{c}_{i_{1}},\cdots,\boldsymbol{c}_{i_{n}}\\}$, indexed by $I$,
naturally defines a multiset
$\\{x_{i_{1}},\cdots,x_{i_{n}}\\}\subseteq\mathcal{D}$, so that the induced
Nyström approximation can in this case be regarded as an approximation induced
by a subset of points in $\mathcal{D}$. Consequently, in the kernel-matrix
framework, instead of relying only on subsets of columns, we may more
generally consider Nyström approximations defined from a multiset
$\mathcal{S}\subseteq\mathscr{X}$. Using matrix notation, the Nyström
approximation of $\mathbf{K}$ defined by a subset
$\mathcal{S}=\\{s_{1},\cdots,s_{n}\\}$ is the $N\times N$ SPSD matrix
$\hat{\mathbf{K}}(\mathcal{S})$, with $i,j$ entry

$\big{[}\hat{\mathbf{K}}(\mathcal{S})\big{]}_{i,j}=\mathbf{k}^{T}(x_{i})\mathbf{K}_{\mathcal{S}}^{\dagger}\mathbf{k}(x_{j}),$
(2)

where $\mathbf{K}_{\mathcal{S}}$ is the $n\times n$ kernel matrix defined by
the kernel $K$ and the subset $\mathcal{S}$, and where

$\mathbf{k}(x)=\big{(}K(x,s_{1}),\cdots,K(x,s_{n})\big{)}^{T}\in\mathbb{R}^{n}.$

We shall refer to such a set or multiset $\mathcal{S}$ as a Nyström sample,
and to the elements of $\mathcal{S}$ as landmark points; the notation
$\hat{\mathbf{K}}(\mathcal{S})$ emphasises that the considered Nyström
approximation of $\mathbf{K}$ is induced by $\mathcal{S}$. As in the column-
sampling case, the landmark-point-based framework naturally raises questions
related to the characterisation and the design of efficient Nyström samples
$\mathcal{S}$ (i.e. leading to accurate approximations of $\mathbf{K}$). As an
interesting feature, Nyström samples of size $n$ may be regarded as elements
of $\mathscr{X}^{n}$, and if the underlying set $\mathscr{X}$ is regular
enough, they might be directly optimised on $\mathscr{X}^{n}$; the situation
we consider in this work corresponds to the case $\mathscr{X}=\mathbb{R}^{d}$,
with $d\in\mathbb{N}$, but $\mathscr{X}$ may more generally be a
differentiable manifold.

######  Remark 1.1.

If we denote by $\mathcal{H}$ the reproducing kernel Hilbert space (RKHS, see
e.g. [1, 14]) of real-valued functions on $\mathscr{X}$ associated with $K$,
we may then note that the matrix $\hat{\mathbf{K}}(\mathcal{S})$ is the kernel
matrix defined by $K_{S}$ and the set $\mathcal{D}$, with $K_{S}$ the
reproducing kernel of the subspace

$\mathcal{H}_{S}=\operatorname{span}\\{k_{s_{1}},\cdots,k_{s_{n}}\\}\subseteq\mathcal{H},$

where, for $t\in\mathscr{X}$, the function $k_{t}\in\mathcal{H}$ is defined as
$k_{t}(x)=K(x,t)$, for all $x\in\mathscr{X}$. $\triangleleft$

###  1.2 Assessing the Accuracy of Nyström Approximations

In the classical literature on the Nyström approximation of SPSD matrices, the
accuracy of the approximation induced by a Nyström sample $\mathcal{S}$ is
often assessed through the following criteria:

  1. (C.1)

$\big{\|}\mathbf{K}-\hat{\mathbf{K}}(\mathcal{S})\big{\|}_{*}$, with
$\|.\|_{*}$ the trace norm;

  2. (C.2)

$\big{\|}\mathbf{K}-\hat{\mathbf{K}}(\mathcal{S})\|_{\mathrm{F}}$, with
$\|.\|_{\mathrm{F}}$ the Frobenius norm;

  3. (C.3)

$\big{\|}\mathbf{K}-\hat{\mathbf{K}}(\mathcal{S})\big{\|}_{2}$, with
$\|.\|_{2}$ the spectral norm.

Although defining relevant and easily interpretable measures of the
approximation error, these criteria are relatively costly to evaluate. Indeed,
each of them involves the inversion or pseudoinversion of the kernel matrix
$\mathbf{K}_{\mathcal{S}}$, with complexity $\mathcal{O}(n^{3})$. The
evaluation of the criterion (C.1) also involves the computation of the $N$
diagonal entries of $\hat{\mathbf{K}}(\mathcal{S})$, leading to an overall
complexity of $\mathcal{O}(n^{3}+Nn^{2})$. The evaluation of (C.2) involves
the full construction of the matrix $\hat{\mathbf{K}}(\mathcal{S})$, with an
overall complexity of $\mathcal{O}(n^{3}+n^{2}N^{2})$, and the evaluation of
(C.3) in addition requires the computation of the largest eigenvalue of an
$N\times N$ SPSD matrix, leading to an overall complexity of
$\mathcal{O}(n^{3}+n^{2}N^{2}+N^{3})$. If $\mathscr{X}=\mathbb{R}^{d}$, then
the evaluation of the partial derivatives of these criteria (regarded as maps
from $\mathscr{X}^{n}$ to $\mathbb{R}$) with respect to a single coordinate of
a landmark point has a complexity similar to the complexity of evaluating the
criteria themselves. As a result, a direct optimisation of these criteria over
$\mathscr{X}^{n}$ is intractable in most practical applications.

###  1.3 Radial Squared-Kernel Discrepancy

As a surrogate for the criteria (C.1)-(C.3), and following the connections
between the Nyström approximation of SPSD matrices, the approximation of
integral operators with SPSD kernels and the kernel embedding of measures, we
consider the following radial squared-kernel discrepancy criterion (radial
SKD, see [9, 7]), denoted by $R$ and given by, for
$\mathcal{S}=\\{s_{1},\cdots,s_{n}\\}$,

$R(\mathcal{S})=\|\mathbf{K}\|_{\mathrm{F}}^{2}-\frac{1}{\|\mathbf{K}_{\mathcal{S}}\|_{\mathrm{F}}^{2}}\bigg{(}\sum_{i=1}^{N}\sum_{j=1}^{n}K^{2}(x_{i},s_{j})\bigg{)}^{2},\text{
if $\|\mathbf{K}_{\mathcal{S}}\|_{\mathrm{F}}>0$},$ (3)

and $R(\mathcal{S})=\|\mathbf{K}\|_{\mathrm{F}}^{2}$ if
$\|\mathbf{K}_{\mathcal{S}}\|_{\mathrm{F}}=0$; the notation
$K^{2}(x_{i},s_{j})$ stands for $\big{(}K(x_{i},s_{j})\big{)}^{2}$. We may
note that $R(\mathcal{S})\geqslant 0$. In (3), the evaluation of the term
$\|\mathbf{K}\|_{\mathrm{F}}^{2}$ has complexity $\mathcal{O}(N^{2})$;
nevertheless, this term does not depend on the Nyström sample $\mathcal{S}$,
and may thus be regarded as a constant. The complexity of the evaluation of
the term $R(\mathcal{S})-\|\mathbf{K}\|_{\mathrm{F}}^{2}$, i.e. of the radial
SKD up to the constant $\|\mathbf{K}\|_{\mathrm{F}}^{2}$, is
$\mathcal{O}(n^{2}+nN)$, and the same holds for the complexity of the
evaluation of the partial derivative of $R(\mathcal{S})$ with respect to a
coordinate of a landmark point, see equation (5) below. We may in particular
note that the evaluation of the radial SKD criterion or its partial
derivatives does not involve the inversion or pseudoinversion of the $n\times
n$ matrix $\mathbf{K}_{\mathcal{S}}$.

######  Remark 1.2.

From a theoretical standpoint, the radial SKD criterion measures the distance,
in the Hilbert space of all Hilbert-Schmidt operators on $\mathcal{H}$,
between the integral operator corresponding to the initial matrix
$\mathbf{K}$, and the projection of this operator onto the subspace spanned by
an integral operator defined from the kernel $K$ and a uniform measure on
$\mathcal{S}$. The radial SKD may also be defined for non-uniform measures,
and the criterion in this case depends not only on $\mathcal{S}$, but also on
a set of relative weights associated with each landmark point in
$\mathcal{S}$; in this work, we only focus on the uniform-weight case. See [9,
7] for more details. $\triangleleft$

The following inequalities hold:

$\big{\|}\mathbf{K}-\hat{\mathbf{K}}(\mathcal{S})\big{\|}_{2}^{2}\leqslant\big{\|}\mathbf{K}-\hat{\mathbf{K}}(\mathcal{S})\big{\|}_{\mathrm{F}}^{2}\leqslant
R(\mathcal{S})\leqslant\|\mathbf{K}\|_{\mathrm{F}}^{2},\quad\text{and}\quad\frac{1}{N}\big{\|}\mathbf{K}-\hat{\mathbf{K}}(\mathcal{S})\big{\|}_{*}^{2}\leqslant\big{\|}\mathbf{K}-\hat{\mathbf{K}}(\mathcal{S})\big{\|}_{\mathrm{F}}^{2},$

which, in complement to the theoretical properties enjoyed by the radial SKD,
further support the use of the radial SKD as a numerically affordable
surrogate for (C.1)-(C.3) (see also the numerical experiments in Section 4).

From now on, we assume that $\mathscr{X}=\mathbb{R}^{d}$. Let $[s]_{l}$, with
$l\in\\{1,\cdots,d\\}$, be the $l$-th coordinate of $s$ in the canonical basis
of $\mathscr{X}=\mathbb{R}^{d}$. For $x\in\mathscr{X}$, we denote by (assuming
they exist)

$\partial^{[\mathrm{l}]}_{[s]_{l}}K^{2}(s,x)\quad\text{ and
}\quad\partial^{[\mathrm{d}]}_{[s]_{l}}K^{2}(s,s)$ (4)

the partial derivatives of the maps $s\mapsto K^{2}(s,x)$ and $s\mapsto
K^{2}(s,s)$ at $s$ and with respect to the $l$-th coordinate of $s$,
respectively; the notation $\partial^{[\mathrm{l}]}$ indicates that the left
entry of the kernel is considered, while $\partial^{[\mathrm{d}]}$ refers to
the diagonal of the kernel; we use similar notations for any kernel function
on $\mathscr{X}\times\mathscr{X}$.

For a fixed number of landmark points $n\in\mathbb{N}$, the radial SKD
criterion can be regarded as a function from $\mathscr{X}^{n}$ to
$\mathbb{R}$. For a Nyström sample
$\mathcal{S}=\\{s_{1},\cdots,s_{n}\\}\in\mathscr{X}^{n}$, and for
$k\in\\{1,\cdots,n\\}$ and $l\in\\{1,\cdots,d\\}$, we denote by
$\partial_{[s_{k}]_{l}}R(\mathcal{S})$ the partial derivative of the map
$R:\mathscr{X}^{n}\to\mathbb{R}$ at $\mathcal{S}$ with respect to the $l$-th
coordinate of the $k$-th landmark point $s_{k}\in\mathscr{X}$. We have

$\displaystyle\begin{split}\partial_{[s_{k}]_{l}}R(\mathcal{S})&=\frac{1}{\|\mathbf{K}_{\mathcal{S}}\|_{\mathrm{F}}^{4}}\bigg{(}\sum_{i=1}^{N}\sum_{j=1}^{n}K^{2}(s_{j},x_{i})\bigg{)}^{2}\bigg{(}\partial_{[s_{k}]_{l}}^{[\mathrm{d}]}K^{2}(s_{k},s_{k})+2\sum_{\begin{subarray}{c}j=1,\\\
j\neq
k\end{subarray}}^{n}\partial_{[s_{k}]_{l}}^{[\mathrm{l}]}K^{2}(s_{k},s_{j})\bigg{)}\\\
&\quad\quad\phantom{x}-\frac{2}{\|\mathbf{K}_{\mathcal{S}}\|_{\mathrm{F}}^{2}}\bigg{(}\sum_{i=1}^{N}\sum_{j=1}^{n}K^{2}(s_{j},x_{i})\bigg{)}\bigg{(}\sum_{i=1}^{N}\partial_{[s_{k}]_{l}}^{[\mathrm{l}]}K^{2}(s_{k},x_{i})\bigg{)}.\end{split}$
(5)

In this work, we investigate the possibility to use the partial derivatives
(5), or stochastic approximations of these derivatives, to directly optimise
the radial SKD criterion $R$ over $\mathscr{X}^{n}$ via gradient or stochastic
gradient descent; the stochastic approximation schemes we consider aim at
reducing the burden of the numerical cost induced by the evaluation of the
partial derivatives of $R$ when $N$ is large.

The document is organised as follows. In Section 2, we discuss the convergence
of a gradient descent with fixed step size for the minimisation of $R$ over
$\mathscr{X}^{n}$. The stochastic approximation of the gradient of the radial
SKD criterion (3) is discussed in Section 3, and some numerical experiments
are carried out in Section 4. Section 5 consists of a concluding discussion,
and the Appendix contains a proof of Theorem 2.1.

##  2 A Convergence Result

We use the same notation as in Section 1.3 (in particular, we still assume
that $\mathscr{X}=\mathbb{R}^{d}$), and by analogy with (4), for $s$ and
$x\in\mathscr{X}$, and for $l\in\\{1,\cdots,d\\}$, we denote by
$\partial^{[\mathrm{r}]}_{[s]_{l}}K^{2}(x,s)$ the partial derivative of the
map $s\mapsto K^{2}(x,s)$ with respect to the $l$-th coordinate of $s$. Also,
for a fixed $n\in\mathbb{N}$, we denote by $\nabla
R(\mathcal{S})\in\mathscr{X}^{n}=\mathbb{R}^{nd}$ the gradient of
$R:\mathscr{X}^{n}\to\mathbb{R}$ at $\mathcal{S}$; in matrix notation, we have

$\nabla
R(\mathcal{S})=\Big{(}\big{(}\nabla_{s_{1}}R(\mathcal{S})\big{)}^{T},\cdots,\big{(}\nabla_{s_{1}}R(\mathcal{S})\big{)}^{T}\Big{)}^{T},$

with
$\nabla_{s_{k}}R(\mathcal{S})=\big{(}\partial_{[s_{k}]_{1}}R(\mathcal{S}),\cdots,\partial_{[s_{k}]_{d}}R(\mathcal{S})\big{)}^{T}\in\mathbb{R}^{d}$
for $k\in\\{1,\cdots,n\\}$.

######  Theorem 2.1.

We make the following assumptions on the squared-kernel $K^{2}$, which we
assume hold for all $x$ and $y\in\mathscr{X}=\mathbb{R}^{d}$, and all $l$ and
$l^{\prime}\in\\{1,\cdots,d\\}$, uniformly:

  1. (C.1)

there exists $\alpha>0$ such that $K^{2}(x,x)\geqslant\alpha$;

  2. (C.2)

there exists $M_{1}>0$ such that
$\bBigg<EMAIL_ADDRESS>M_{1}$ and
$\bBigg<EMAIL_ADDRESS>M_{1}$;

  3. (C.3)

there exists $M_{2}>0$ such that
$\bBigg<EMAIL_ADDRESS>M_{2}$,
$\bBigg<EMAIL_ADDRESS>M_{2}$ and
$\bBigg<EMAIL_ADDRESS>M_{2}$.

Let $\mathcal{S}$ and $\mathcal{S}^{\prime}\in\mathbb{R}^{nd}$ be two Nyström
samples; under the above assumptions, there exists $L>0$ such that

$\big{\|}\nabla R(\mathcal{S})-\nabla R(\mathcal{S}^{\prime})\big{\|}\leqslant
L\big{\|}\mathcal{S}-\mathcal{S}^{\prime}\big{\|}$

with $\|.\|$ the Euclidean norm of $\mathbb{R}^{nd}$; in other words, the
gradient of $R:\mathbb{R}^{nd}\to\mathbb{R}$ is Lipschitz-continuous with
Lipschitz constant $L$.

Since $R$ is bounded from below, for $0<\gamma\leqslant 1/L$ and independently
of the considered initial Nyström sample $\mathcal{S}^{(0)}$, Theorem 2.1
entails that a gradient descent from $\mathcal{S}^{(0)}$, with fixed stepsize
$\gamma$ for the minimisation of $R$ over $\mathscr{X}^{n}$, produces a
sequence of iterates that converges to a critical point of $R$. Barring some
specific and largely pathological cases, the resulting critical point is
likely to be a local minimum of $R$, see for instance [12]. See the Appendix
for a proof of Theorem 2.1.

The conditions considered in Theorem 2.1 ensure the existence of a general
Lipschitz constant $L$ for the gradient of $R$; they, for instance, hold for
all sufficiently regular Matérn kernels (thus including the Gaussian or
squared-exponential kernel). These conditions are only sufficient conditions
for the convergence of a gradient descent for the minimisation of $R$. By
introducing additional problem-dependent conditions, some convergence results
might be obtained for more general squared kernels $K^{2}$ and adequate
initial Nyström samples $\mathcal{S}^{(0)}$. For instance, the condition (C.1)
simply aims at ensuring that
$\|\mathbf{K}_{\mathcal{S}}\|_{\mathrm{F}}^{2}\geqslant n\alpha>0$ for all
$\mathcal{S}\in\mathscr{X}^{n}$; this condition might be relaxed to account
for kernels with vanishing diagonal, but one might then need to introduce ad
hoc conditions to ensure that $\|\mathbf{K}_{\mathcal{S}}\|_{\mathrm{F}}^{2}$
remains large enough during the minimisation process.

##  3 Stochastic Approximation of the Radial SKD Gradient

The complexity of evaluating a partial derivative of
$R:\mathscr{X}^{n}\to\mathbb{R}$ is $\mathcal{O}(n^{2}+nN)$, which might
become prohibitive for large values of $N$. To overcome this limitation,
stochastic approximations of the gradient of $R$ might be considered (see e.g.
[2]).

The evaluation of (5) involves, for instance, terms of the form
$\sum_{i=1}^{N}K^{2}(s,x_{i})$, with $s\in\mathscr{X}$ and
$\mathcal{D}=\\{x_{1},\cdots,x_{N}\\}$. Introducing a random variable $X$ with
uniform distribution on $\mathcal{D}$, we can note that

$\sum_{i=1}^{N}K^{2}(s,x_{i})=N\mathbb{E}\big{[}K^{2}(s,X)\big{]},$

and the mean $\mathbb{E}[K^{2}(s,X)]$ may then, classically, be approximated
by random sampling. More precisely, if $X_{1},\cdots,X_{b}$ are
$b\in\mathbb{N}$ copies of $X$, we have

$\mathbb{E}\big{[}K^{2}(s,X)\big{]}=\frac{1}{b}\sum_{j=1}^{b}\mathbb{E}\big{[}K^{2}(s,X_{j})\big{]}\quad\text{and}\quad\mathbb{E}\big{[}\partial_{[s]_{l}}^{[\mathrm{l}]}K^{2}(s,X)\big{]}=\frac{1}{b}\sum_{j=1}^{b}\mathbb{E}\big{[}\partial_{[s]_{l}}^{[\mathrm{l}]}K^{2}(s,X_{j})\big{]},$

so that we can easily define unbiased estimators of the various terms
appearing in (5). We refer to the sample size $b$ as the batch size.

Let $k\in\\{1,\ldots,n\\}$ and $l\in\\{1,\ldots,d\\}$; the partial derivative
(5) can be rewritten as

$\partial_{[s_{k}]_{l}}R(\mathcal{S})=\frac{T_{1}^{2}}{\|\mathbf{K}_{\mathcal{S}}\|_{\mathrm{F}}^{4}}\Upsilon(\mathcal{S})-\frac{2T_{1}T_{2}^{k,l}}{\|\mathbf{K}_{\mathcal{S}}\|_{\mathrm{F}}^{2}},$

with $T_{1}=\sum_{i=1}^{N}\sum_{j=1}^{n}K^{2}(s_{j,}x_{i})$ and
$T_{2}^{k,l}=\sum_{i=1}^{N}\partial_{[s_{k}]_{l}}^{[\mathrm{l}]}K^{2}(s_{k},x_{i})$,
and

$\Upsilon(\mathcal{S})=\partial_{[s_{k}]_{l}}^{[\mathrm{d}]}K^{2}(s_{k},s_{k})+2\sum_{\begin{subarray}{c}j=1,\\\
j\neq
k\end{subarray}}^{n}\partial_{[s_{k}]_{l}}^{[\mathrm{l}]}K^{2}(s_{k},s_{j}).$

The terms $T_{1}$ and $T_{2}^{k,l}$ are the only terms in (5) that depend on
$\mathcal{D}$. From a uniform random sample
$\mathbf{X}=\\{X_{1},\cdots,X_{b}\\}$, we define the unbiased estimators
$\hat{T}_{1}(\mathbf{X})$ of $T_{1}$, and $\hat{T}_{2}^{k,l}(\mathbf{X})$ of
$T_{2}^{k,l}$, as

$\hat{T}_{1}(\mathbf{X})=\frac{N}{b}\sum_{i=1}^{n}\sum_{j=1}^{b}K^{2}(s_{i},X_{j}),\quad\text{
and
}\quad\hat{T}_{2}^{k,l}(\mathbf{X})=\frac{N}{b}\sum_{j=1}^{b}\partial_{[s_{k}]_{l}}^{[\mathrm{l}]}K^{2}(s_{k},X_{j}).$

In what follows, we discuss the properties of some stochastic approximations
of the gradient of $R$ that can be defined from such estimators.

#### One-Sample Approximation.

Using a single random sample $\mathbf{X}=\\{X_{1},\cdots,X_{b}\\}$ of size
$b$, we can define the following stochastic approximation of the partial
derivative (5):

$\hat{\partial}_{[s_{k}]_{l}}R(\mathcal{S};\mathbf{X})=\frac{\hat{T}_{1}(\mathbf{X})^{2}}{\|\mathbf{K}_{\mathcal{S}}\|_{\mathrm{F}}^{4}}\Upsilon(\mathcal{S})-\frac{2\hat{T}_{1}(\mathbf{X})\hat{T}_{2}^{k,l}(\mathbf{X})}{\|\mathbf{K}_{\mathcal{S}}\|_{\mathrm{F}}^{2}}.$
(6)

An evaluation of $\hat{\partial}_{[s_{k}]_{l}}R(\mathcal{S};\mathbf{X})$ has
complexity $\mathcal{O}(n^{2}+nb)$, as opposed to $\mathcal{O}(n^{2}+nN)$ for
the corresponding exact partial derivative. However, due to the dependence
between $\hat{T}_{1}(\mathbf{X})$ and $\hat{T}_{2}^{k,l}(\mathbf{X})$, and to
the fact that $\hat{\partial}_{[s_{k}]_{l}}R(\mathcal{S};\mathbf{X})$ involves
the square of $\hat{T}_{1}(\mathbf{X})$, the stochastic partial derivative
$\hat{\partial}_{[s_{k}]_{l}}R(\mathcal{S};\mathbf{X})$ will generally be a
biased estimator of $\partial_{[s_{k}]_{l}}R(\mathcal{S})$.

#### Two-Sample Approximation.

To obtain an unbiased estimator of the partial derivative (5), instead of
considering a single random sample, we may define a stochastic approximation
based on two independent random samples
$\mathbf{X}=\\{X_{1},\cdots,X_{b_{\mathbf{X}}}\\}$ and
$\mathbf{Y}=\\{Y_{1},\cdots,Y_{b_{\mathbf{Y}}}\\}$, consisting of
$b_{\mathbf{X}}$ and $b_{\mathbf{Y}}\in\mathbb{N}$ copies of $X$ (i.e.
consisting of uniform random variables on $\mathcal{D}$), with
$b=b_{\mathbf{X}}+b_{\mathbf{Y}}$. The two-sample estimator of (5) is then
given by

$\hat{\partial}_{[s_{k}]_{l}}R(\mathcal{S};\mathbf{X},\mathbf{Y})=\frac{\hat{T}_{1}(\mathbf{X})\hat{T}_{1}(\mathbf{Y})}{\|\mathbf{K}_{\mathcal{S}}\|_{\mathrm{F}}^{4}}\Upsilon(\mathcal{S})-\frac{2\hat{T}_{1}(\mathbf{X})\hat{T}_{2}^{k,l}(\mathbf{Y})}{\|\mathbf{K}_{\mathcal{S}}\|_{\mathrm{F}}^{2}},$
(7)

and since
$\mathbb{E}\big{[}\hat{T}_{1}(\mathbf{X})\hat{T}_{1}(\mathbf{Y})\big{]}=T_{1}^{2}$
and
$\mathbb{E}\big{[}\hat{T}_{1}(\mathbf{X})\hat{T}_{2}^{k,l}(\mathbf{Y})\big{]}=T_{1}T_{2}^{k,l}$,
we have

$\mathbb{E}\Big{[}\hat{\partial}_{[s_{k}]_{l}}R(\mathcal{S};\mathbf{X},\mathbf{Y})\Big{]}=\partial_{[s_{k}]_{l}}R(\mathcal{S}).$

Although being unbiased, for a common batch size $b$, the variance of the two-
sample estimator (7) will generally be larger than the variance of the one-
sample estimator (6). In our numerical experiments, the larger variance of the
unbiased estimator (7) seems to actually slow down the descent when compared
to the descent obtained with the one-sample estimator (6).

######  Remark 3.1.

While considering two independent samples $\mathbf{X}$ and $\mathbf{Y}$, the
two terms $\hat{T}_{1}(\mathbf{X})\hat{T}_{1}(\mathbf{Y})$ and
$\hat{T}_{1}(\mathbf{X})\hat{T}_{2}^{k,l}(\mathbf{Y})$ appearing in (7) are
dependent. This dependence may complicate the analysis of the properties of
the resulting SGD; nevertheless, this issue might be overcome by considering
four independent samples instead of two. $\triangleleft$

##  4 Numerical Experiments

Throughout this section, the matrices $\mathbf{K}$ are defined from multisets
$\mathcal{D}=\\{x_{1},\cdots,x_{N}\\}\subset\mathbb{R}^{d}$ and from kernels
$K$ of the form $K(x,t)=e^{-\rho\|x-t\|^{2}}$, with $\rho>0$ and where $\|.\|$
is the Euclidean norm of $\mathbb{R}^{d}$ (Gaussian kernel). Except for the
synthetic example of Section 4.1, all the multisets $\mathcal{D}$ we consider
consist of the entries of data sets available on the UCI Machine Learning
Repository; see [6].

Our experiments are based on the following protocol: for a given
$n\in\mathbb{N}$, we consider an initial Nyström sample $\mathcal{S}^{(0)}$
consisting of $n$ points drawn uniformly at random, without replacement, from
$\mathcal{D}$. The initial sample $\mathcal{S}^{(0)}$ is regarded as an
element of $\mathscr{X}^{n}$, and used to initialise a GD or SGD, with fixed
stepsize $\gamma>0$, for the minimisation of $R$ over $\mathscr{X}^{n}$,
yielding, after $T\in\mathbb{N}$ iterations, a locally optimised Nyström
sample $\mathcal{S}^{(T)}$. The SGDs are performed with the one-sample
estimator (6) and are based on independent and identically distributed uniform
random variables on $\mathcal{D}$ (i.e. i.i.d. sampling), with batch size
$b\in\mathbb{N}$; see Section 3. We assess the accuracy of the Nyström
approximations of $\mathbf{K}$ induced by $\mathcal{S}^{(0)}$ and
$\mathcal{S}^{(T)}$ in terms of radial SKD and of the classical criteria
(C.1)-(C.3).

For a Nyström sample $\mathcal{S}\in\mathscr{X}^{n}$ of size $n\in\mathbb{N}$,
the matrix $\hat{\mathbf{K}}(\mathcal{S})$ is of rank at most $n$. Following
[10, 4], to further assess the efficiency of the approximation of $\mathbf{K}$
induced by $\mathcal{S}$, we introduce the approximation factors

$\mathcal{E}_{\mathrm{tr}}(\mathcal{S})=\frac{\|\mathbf{K}-\hat{\mathbf{K}}(\mathcal{S})\|_{*}}{\|\mathbf{K}-\mathbf{K}_{n}\|_{*}},\quad\mathcal{E}_{\mathrm{F}}(\mathcal{S})=\frac{\|\mathbf{K}-\hat{\mathbf{K}}(\mathcal{S})\|_{\mathrm{F}}}{\|\mathbf{K}-\mathbf{K}_{n}\|_{\mathrm{F}}},\quad\text{
and
}\mathcal{E}_{\mathrm{sp}}(\mathcal{S})=\frac{\|\mathbf{K}-\hat{\mathbf{K}}(\mathcal{S})\|_{2}}{\|\mathbf{K}-\mathbf{K}_{n}\|_{2}},$
(8)

where $\mathbf{K}_{n}$ denotes an optimal rank-$n$ approximation of
$\mathbf{K}$ (i.e. the approximation of $\mathbf{K}$ obtained by truncation of
a spectral expansion of $\mathbf{K}$ and based on $n$ of the largest
eigenvalues of $\mathbf{K}$). The closer
$\mathcal{E}_{\mathrm{tr}}(\mathcal{S})$,
$\mathcal{E}_{\mathrm{F}}(\mathcal{S})$ and
$\mathcal{E}_{\mathrm{sp}}(\mathcal{S})$ are to $1$, the more efficient the
approximation is.

###  4.1 Bi-Gaussian Example

We consider a kernel matrix $\mathbf{K}$ defined by a set $\mathcal{D}$
consisting of $N=2{,}000$ points in $[-1,1]^{2}\subset\mathbb{R}^{2}$ (i.e.
$d=2$); for the kernel parameter, we use $\rho=1$. A graphical representation
of the set $\mathcal{D}$ is given in Figure 1; it consists of $N$ independent
realisations of a bivariate random variable whose density is proportional to
the restriction of a bi-Gaussian density to the set $[-1,1]^{2}$ (the two
modes of the underlying distribution are located at $(-0.8,0.8)$ and
$(0.8,-0.8)$, and the covariance matrix of the each Gaussian density is
$\mathbb{I}_{2}/2$, with $\mathbb{I}_{2}$ the $2\times 2$ identity matrix).

Figure 1: Graphical representation of the path followed by the landmark points
of a Nystrom sample during the local minimisation of $R$ through GD, with
$n=50$, $\gamma=10^{-6}$ and $T=1{,}300$; the green squares are the landmark
points of the initial sample $\mathcal{S}^{(0)}$, the red dots are the
landmark points of the locally optimised sample $\mathcal{S}^{(T)}$, and the
purple lines correspond to the paths followed by each landmark point (left).
The corresponding decay of the radial SKD is also presented (right).

The initial samples $\mathcal{S}^{(0)}$ are optimised via GD with stepsize
$\gamma=10^{-6}$ and for a fixed number of iterations $T$. A graphical
representation of the paths followed by the landmark points during the
optimisation process is given in Figure 1 (for $n=50$ and $T=1{,}300$); we
observe that the landmark points exhibit a relatively complex dynamic, some of
them showing significant displacements from their initial positions. The
optimised landmark points concentrate around the regions where the density of
points in $\mathcal{D}$ is the largest, and inherit a space-filling-type
property in accordance with the stationarity of the kernel $K$.

To assess the improvement yielded by the optimisation process, for a given
number of landmark points $n\in\mathbb{N}$, we randomly draw an initial
Nyström sample $\mathcal{S}^{(0)}$ from $\mathcal{D}$ (uniform sampling
without replacement) and compute the corresponding locally optimised sample
$\mathcal{S}^{(T)}$ (GD with $\gamma=10^{-6}$ and $T=1{,}000$). We then
compare $R\big{(}\mathcal{S}^{(0)}\big{)}$ with
$R\big{(}\mathcal{S}^{(T)}\big{)}$, and compute the corresponding
approximation factors with respect to the trace, Frobenius and spectral norms,
see (8). We consider three different values of $n$, namely $n=20$, $50$ and
$80$, and each time perform $m=1{,}000$ repetitions of this experiment. Our
results are presented in Figure 2; we observe that, independently of $n$, the
local optimisation produces a significant improvement of the Nyström
approximation accuracy for all the criterion considered; the improvements are
particularly noticeable for the trace and Frobenius norms, and slightly less
for the spectral norm (which of the three, appears the coarsest measure of the
approximation accuracy). Remarkably, the efficiencies of the locally optimised
Nyström samples are relatively close to each other, in particular in terms of
trace and Frobenius norms, suggesting that a large proportion of the local
minima of the radial SKD induce approximations of comparable quality.

Figure 2: For the Bi-Gaussian example, comparison of the efficiency of the
Nyström approximations for the initial samples $\mathcal{S}^{(0)}$ and the
locally optimised samples $\mathcal{S}^{(T)}$ (optimisation through GD with
$\gamma=10^{-6}$ and $T=1{,}000$). Each row corresponds to a given value of
$n$; in each case $m=1{,}000$ repetitions are performed. The first column
corresponds to the radial SKD, and the following three correspond to the
approximation factors defined in (8).

###  4.2 Abalone Data Set

We now consider the $d=8$ attributes of the Abalone data set. After removing
two observations that are clear outliers, we are left with $N=4{,}175$
entries. Each of the $8$ features is standardised such that it has zero mean
and unit variance. We set $n=50$ and consider three different values of the
kernel paramater $\rho$, namely $\rho=0.25$, $1$, and $4$; this values are
chosen so that the eigenvalues of the kernel matrix $\mathbf{K}$ exhibit
sharp, moderate and shallower decays, respectively. For the Nyström sample
optimisation, we use SGD with i.i.d. sampling and batch size $b=50$,
$T=10{,}000$ and $\gamma=8\times 10^{-7}$; these values were chosen to obtain
relatively efficient optimisations for the whole range of values of $\rho$ we
consider. For each value of $\rho$, we perform $m=200$ repetitions. The
results are presented in Figure 3.

Figure 3: For the Abalone data set with $n=50$ and $\rho\in\\{0.25,1,4\\}$,
comparison of the efficiency of the Nyström approximations for the initial
Nyström samples $\mathcal{S}^{(0)}$ and the locally optimised samples
$\mathcal{S}^{(T)}$ (SGD with i.i.d sampling, $b=50$, $\gamma=8\times 10^{-7}$
and $T=10{,}000$). Each row corresponds to a given value of $\rho$; in each
case, $m=200$ repetitions are performed.

We observe that regardless of the values of $\rho$ and in comparison with the
initial Nyström samples, the efficiencies of the locally optimised samples in
terms of trace, Frobenius and spectral norms are significantly improved. As
observed in Section 4.1, the gains yielded by the local optimisations are more
evident in terms of trace and Frobenius norms, and the impact of the
initialisation appears limited.

###  4.3 MAGIC Data Set

We consider the $d=10$ attributes of the MAGIC Gamma Telescope data set. In
pre-processing, we remove the $115$ duplicated entries in the data set,
leaving us with $N=18{,}905$ data points; we then standardise each of the
$d=10$ features of the data set. For the kernel parameter, we use $\rho=0.2$.

In Figure 4, we present the results obtained after the local optimisation of
$m=200$ random initial Nyström samples of size $n=100$ and $200$. Each
optimisation was performed through SGD with i.i.d. sampling, batch size $b=50$
and stepsize $\gamma=5\times 10^{-8}$; as number of iterations, for $n=100$,
we used $T=3{,}000$, and $T=4{,}000$ for $n=200$. The optimisation parameters
were chosen to obtain relatively efficient but not fully completed descents,
as illustrated in Figure 4. Alongside the radial SKD, we only compute the
approximation factor corresponding to the trace norm (the trace norm is indeed
the least costly to evaluate of the three matrix norms we consider, see
Section 1.2). As in the previous experiments, we observe a significant
improvement of the initial Nyström samples obtained by local optimisation of
the radial SKD.

Figure 4: For the MAGIC data set, boxplots of the radial SKD $R$ and of the
approximation factor $\mathcal{E}_{\mathrm{tr}}$ before and after the local
optimisation via SGD of random Nyström samples of size $n=100$ and $200$; for
each value of $n$, $m=200$ repetitions are performed. The SGD is based on
i.i.d. sampling, with $b=50$ and $\gamma=5\times 10^{-8}$; for $n=100$, the
descent is stopped after $T=3{,}000$ iterations, and after $T=4{,}000$
iterations for $n=200$ (left). A graphical representation of the decay of the
radial SKD is also presented for $n=200$ (right).

###  4.4 MiniBooNE Data Set

In this last experiment, we consider the $d=50$ attributes of the MiniBooNE
particle identification data set. In pre-processing, we remove the $471$
entries in the data set with missing values, and $1$ entry appearing as a
clear outlier, leaving us with $N=129{,}592$ data points; we then standardise
each of the $d=50$ features of the data set. We use $\rho=0.04$ (kernel
parameter).

Figure 5: For the MiniBooNE data set, decay of the radial SKD during the
optimisation of a random initial Nyström sample of size $n=1{,}000$. The SGD
is based on i.i.d. sampling with batch size $b=200$ and stepsize
$\gamma=2\times 10^{-7}$, and the descent is stopped after $T=8{,}000$
iterations; the cost is evaluated every $100$ iterations.

We consider a random initial Nyström sample of size $n=1{,}000$, and optimise
it through SGD with i.i.d. sampling, batch size $b=200$, stepsize
$\gamma=2\times 10^{-7}$; the descent is stopped after $T=8{,}000$ iterations.
The resulting decay of the radial SKD is presented in Figure 5 (the cost is
evaluated every $100$ iterations), and the trace norm of the Nyström
approximation error for the initial and locally optimised samples are
reported. In terms of computation time, on our machine (endowed with an 3.5
GHz Dual-Core Intel Core i7 processor, and using a single-threaded C
implementation interfaced with R), for $n=1{,}000$, an evaluation of the
radial SKD (up to the constant $\|\mathbf{K}\|_{\mathrm{F}}^{2}$) takes
$6.8~{}\mathrm{s}$, while an evaluation of the term
$\|\mathbf{K}-\hat{\mathbf{K}}(S)\|_{*}$ takes $6{,}600~{}\mathrm{s}$;
performing the optimisation reported in Figure 5 without checking the decay of
the cost takes $1{,}350~{}\mathrm{s}$. This experiment illustrates the ability
of the considered framework to tackle relatively large problems.

##  5 Conclusion

We demonstrated the relevance of the radial-SKD-based framework for the local
optimisation, through SGD, of Nyström samples for SPSD kernel-matrix
approximation. We studied the Lipschitz continuity of the underlying gradient
and discussed its stochastic approximation. We performed numerical experiments
illustrating that local optimisation of the radial SKD yields significant
improvement of the Nyström approximation in terms of trace, Frobenius and
spectral norms.

In our experiments, we implemented SGD with i.i.d. sampling, fixed stepsize
and fixed number of iterations; although already bringing satisfactory
results, to improve the time efficiency of the approach, the optimisation
strategy could be accelerated by considering for instance adaptive stepsize,
parallelisation or momentum-type techniques (see [16] for an overview). The
initial Nyström samples $\mathcal{S}^{(0)}$ we considered were draw uniformly
at random without replacement; while our experiments suggest that the local
minima of the radial SKD often induce approximations of comparable quality,
the use of more efficient initialisation strategies may be investigated (see
e.g. [13, 11, 18, 3, 4]).

As a side note, when considering the trace norm, the Nyström sampling problem
is intrinsically related to the _integrated-mean-squared-error_ design
criterion in kernel regression (see e.g. [15, 8, 17]); consequently the
approach considered in this paper may be used for the design of experiments
for such models.

## Appendix

###### Proof of Theorem 2.1.

We consider a Nyström sample $\mathcal{S}\in\mathscr{X}^{n}$ and introduce

$c_{\mathcal{S}}=\frac{1}{\|\mathbf{K}_{\mathcal{S}}\|_{\mathrm{F}}^{2}}\sum_{i=1}^{N}\sum_{j=1}^{n}K^{2}(x_{i},s_{j}).$
(9)

In view of (5), the partial derivative of $R$ at $\mathcal{S}$ with respect to
the $l$-th coordinate of the $k$-th landmark point $s_{k}$ can be written as

$\partial_{[s_{k}]_{l}}R(\mathcal{S})=c_{\mathcal{S}}^{2}\bigg{(}\partial_{[s_{k}]_{l}}^{[\mathrm{d}]}K^{2}(s_{k},s_{k})+2\sum_{\begin{subarray}{c}j=1,\\\
j\neq
k\end{subarray}}^{n}\partial_{[s_{k}]_{l}}^{[\mathrm{l}]}K^{2}(s_{k},s_{j})\bigg{)}-2c_{\mathcal{S}}\sum_{i=1}^{N}\partial_{[s_{k}]_{l}}^{[\mathrm{l}]}K^{2}(s_{k},x_{i}).$
(10)

For $k$ and $k^{\prime}\in\\{1,\cdots,n\\}$ with $k\neq k^{\prime}$, and for
$l$ and $l^{\prime}\in\\{1,\cdots,d\\}$, the second-order partial derivatives
of $R$ at $\mathcal{S}$, with respect to the coordinates of the landmark
points in $\mathcal{S}$, verify

$\displaystyle\begin{split}\partial_{[s_{k}]_{l}}\partial_{[s_{k}]_{l^{\prime}}}R(\mathcal{S})&=c_{\mathcal{S}}^{2}\partial_{[s_{k}]_{l}}^{[\mathrm{d}]}\partial_{[s_{k}]_{l^{\prime}}}^{[\mathrm{d}]}K^{2}(s_{k},s_{k})+2c_{\mathcal{S}}(\partial_{[s_{k}]_{l^{\prime}}}c_{\mathcal{S}})\partial_{[s_{k}]_{l}}^{[\mathrm{d}]}K^{2}(s_{k},s_{k})\\\
&\quad+2c_{\mathcal{S}}^{2}\sum_{\begin{subarray}{c}j=1,\\\ j\neq
k\end{subarray}}^{n}\partial_{[s_{k}]_{l}}^{[\mathrm{l}]}\partial_{[s_{k}]_{l^{\prime}}}^{[\mathrm{l}]}K^{2}(s_{k},s_{j})+4c_{\mathcal{S}}(\partial_{[s_{k}]_{l^{\prime}}}c_{\mathcal{S}})\sum_{\begin{subarray}{c}j=1,\\\
j\neq
k\end{subarray}}^{n}\partial_{[s_{k}]_{l}}^{[\mathrm{l}]}K^{2}(s_{k},s_{j})\\\
&\quad-2c_{\mathcal{S}}\sum_{i=1}^{N}\partial_{[s_{k}]_{l}}^{[\mathrm{l}]}\partial_{[s_{k}]_{l^{\prime}}}^{[\mathrm{l}]}K^{2}(s_{k},x_{i})-2(\partial_{[s_{k}]_{l^{\prime}}}c_{\mathcal{S}})\sum_{i=1}^{N}\partial_{[s_{k}]_{l}}^{[\mathrm{l}]}K^{2}(s_{k},x_{i}),\text{
and}\end{split}$ (11)
$\displaystyle\begin{split}\partial_{[s_{k}]_{l}}\partial_{[s_{k^{\prime}}]_{l^{\prime}}}&R(\mathcal{S})=2c_{\mathcal{S}}(\partial_{[s_{k^{\prime}}]_{l^{\prime}}}c_{\mathcal{S}})\partial_{[s_{k}]_{l}}^{[\mathrm{d}]}K^{2}(s_{k},s_{k})+2c_{\mathcal{S}}^{2}\partial_{[s_{k}]_{l}}^{[\mathrm{l}]}\partial_{[s_{k^{\prime}}]_{l^{\prime}}}^{[\mathrm{r}]}K^{2}(s_{k},s_{k^{\prime}})\\\
&\quad+4c_{\mathcal{S}}(\partial_{[s_{k^{\prime}}]_{l^{\prime}}}c_{\mathcal{S}})\sum_{\begin{subarray}{c}j=1,\\\
j\neq
k\end{subarray}}^{n}\partial_{[s_{k}]_{l}}^{[\mathrm{l}]}K^{2}(s_{k},s_{j})-2(\partial_{[s_{k^{\prime}}]_{l^{\prime}}}c_{\mathcal{S}})\sum_{i=1}^{N}\partial_{[s_{k}]_{l}}^{[\mathrm{l}]}K^{2}(s_{k},x_{i}),\end{split}$
(12)

where the partial derivative of $c_{\mathcal{S}}$ with respect to the $l$-th
coordinate of the $k$-th landmark point $s_{k}$ is given by

$\partial_{[s_{k}]_{l}}c_{\mathcal{S}}=\frac{1}{\|\mathbf{K}_{\mathcal{S}}\|_{\mathrm{F}}^{2}}\bigg{(}\sum_{i=1}^{N}\partial_{[s_{k}]_{l}}^{[\mathrm{l}]}K^{2}(s_{k},x_{i})-c_{\mathcal{S}}\partial_{[s_{k}]_{l}}^{[\mathrm{d}]}K^{2}(s_{k},s_{k})-2c_{\mathcal{S}}\sum_{\begin{subarray}{c}j=1,\\\
j\neq
k\end{subarray}}^{n}\partial_{[s_{k}]_{l}}^{[\mathrm{l}]}K^{2}(s_{k},s_{j})\bigg{)}.$
(13)

From (C.1), we have

$\|\mathbf{K}_{\mathcal{S}}\|_{\mathrm{F}}^{2}=\sum_{i=1}^{n}\sum_{j=1}^{n}K^{2}(s_{i},s_{j})\geqslant\sum_{i=1}^{n}K^{2}(s_{i},s_{i})\geqslant
n\alpha.$ (14)

By the Schur product theorem, the squared kernel $K^{2}$ is SPSD; we denote by
$\mathcal{G}$ the RKHS of real-valued functions on $\mathscr{X}$ for which
$K^{2}$ is reproducing. For $x$ and $y\in\mathscr{X}$, we have
$K^{2}(x,y)=\langle k_{x}^{2},k_{y}^{2}\rangle_{\mathcal{G}}$, with
$\langle\cdot,\cdot\rangle_{\mathcal{G}}$ the inner product on $\mathcal{G}$,
and where $k_{x}^{2}\in\mathcal{G}$ is such that $k_{x}^{2}(t)=K^{2}(t,x)$,
for all $t\in\mathscr{X}$. From the Cauchy-Schwartz inequality, we have

$\displaystyle\sum_{i=1}^{N}\sum_{j=1}^{n}K^{2}(s_{j},x_{i})$
$\displaystyle=\sum_{i=1}^{N}\sum_{j=1}^{n}\langle
k_{s_{j}}^{2},k_{x_{i}}^{2}\rangle_{\mathcal{G}}=\bigg{\langle}\sum_{j=1}^{n}k_{s_{j}}^{2},\sum_{i=1}^{N}k_{x_{i}}^{2}\bigg{\rangle}_{\mathcal{G}}$
$\displaystyle\leqslant\bigg{\|}\sum_{j=1}^{n}k_{s_{j}}^{2}\bigg{\|}_{\mathcal{G}}\bigg{\|}\sum_{i=1}^{N}k_{x_{i}}^{2}\bigg{\|}_{\mathcal{G}}=\|\mathbf{K}_{\mathcal{S}}\|_{\mathrm{F}}\|\mathbf{K}\|_{\mathrm{F}}.$
(15)

By combining (9) with inequalities (14) and (15), we obtain

$0\leqslant
c_{\mathcal{S}}\leqslant\frac{\|\mathbf{K}\|_{\mathrm{F}}}{\|\mathbf{K}_{\mathcal{S}}\|_{\mathrm{F}}}\leqslant\frac{\|\mathbf{K}\|_{\mathrm{F}}}{\sqrt{n\alpha}}=C_{0}.$
(16)

Let $k\in\\{1,\ldots,n\\}$ and let $l\in\\{1,\ldots,d\\}$. From equation (13),
and using inequalities (14) and (16) together with (C.2), we obtain

$|\partial_{[s_{k}]_{l}}c_{\mathcal{S}}|\leqslant\frac{M_{1}}{n\alpha}[N+(2n-1)C_{0}]=C_{1}.$
(17)

In addition, let $k^{\prime}\in\\{1,\ldots,n\\}\setminus\\{k\\}$ and
$l^{\prime}\in\\{1,\ldots,d\\}$; from equations (11), (12), (16) and (17), and
conditions (C.2) and (C.3), we get

$\displaystyle|\partial_{[s_{k}]_{l}}\partial_{[s_{k}]_{l^{\prime}}}R(\mathcal{S})|$
$\displaystyle\quad\leqslant
C_{0}^{2}M_{2}+2C_{0}C_{1}M_{1}+2(n-1)C_{0}^{2}M_{2}+4(n-1)C_{0}C_{1}M_{1}+2C_{0}M_{2}N+2C_{1}M_{1}N$
$\displaystyle\quad=(2n-1)C_{0}^{2}M_{2}+(4n-2)C_{0}C_{1}M_{1}+2N(C_{0}M_{2}+C_{1}M_{1}),$
(18)

and

$\displaystyle|\partial_{[s_{k}]_{l}}\partial_{[s_{k^{\prime}}]_{l^{\prime}}}R(\mathcal{S})|$
$\displaystyle\leqslant
2C_{0}C_{1}M_{1}+2C_{0}^{2}M_{2}+4(n-1)C_{0}C_{1}M_{1}+2C_{1}M_{1}N$
$\displaystyle=2C_{0}^{2}M_{2}+(4n-2)C_{0}C_{1}M_{1}+2NC_{1}M_{1}.$ (19)

For $k,k^{\prime}\in\\{1,\ldots,n\\}$, we denote by
$\mathbf{B}^{k,k^{\prime}}$ the $d\times d$ matrix with $l,l^{\prime}$ entry
given by (11) if $k=k^{\prime}$, and by (12) otherwise. The Hessian
$\nabla^{2}R(\mathcal{S})$ can then be represented as a block-matrix, that is

$\nabla^{2}R(\mathcal{S})=\begin{bmatrix}\mathbf{B}^{1,1}&\cdots&\mathbf{B}^{1,n}\\\
\vdots&\ddots&\vdots\\\
\mathbf{B}^{n,1}&\cdots&\mathbf{B}^{n,n}\end{bmatrix}\in\mathbb{R}^{nd\times
nd}.$

The $d^{2}$ entries of the $n$ diagonal blocks of $\nabla^{2}R(\mathcal{S})$
are of the form (11), and the $d^{2}$ entries of the $n(n-1)$ off-diagonal
blocks of $\nabla^{2}R(\mathcal{S})$ are the form (12). From inequalities (18)
and (19), we obtain

$\|\nabla^{2}R(\mathcal{S})\|_{2}^{2}\leqslant\|\nabla^{2}R(\mathcal{S})\|_{\mathrm{F}}^{2}=\sum_{k=1}^{n}\sum_{l=1}^{d}\sum_{l^{\prime}=1}^{d}[\mathbf{B}^{k,k}]_{l,l^{\prime}}^{2}+\sum_{k=1}^{n}\sum_{\begin{subarray}{c}k^{\prime}=1,\\\
k^{\prime}\neq
k\end{subarray}}^{n}\sum_{l=1}^{d}\sum_{l^{\prime}=1}^{d}[\mathbf{B}^{k,k^{\prime}}]_{l,l^{\prime}}^{2}\leqslant
L^{2},$

with

$\displaystyle\begin{split}L&=\big{(}nd^{2}[(2n-1)C_{0}^{2}M_{2}+(4n-2)C_{0}C_{1}M_{1}+2N(C_{0}M_{2}+C_{1}M_{1})]^{2}\\\
&\quad\phantom{x}+4n(n-1)d^{2}[C_{0}^{2}M_{2}+(2n-1)C_{0}C_{1}M_{1}+NC_{1}M_{1}]^{2}\big{)}^{\frac{1}{2}}.\end{split}$

For all $\mathcal{S}\in\mathscr{X}^{n}$, the constant $L$ is an upper bound
for the spectral norm of the Hessian matrix $\nabla^{2}R(\mathcal{S})$, so the
gradient of $R$ is Lipschitz continuous over $\mathscr{X}^{n}$, with Lipschitz
constant $L$. ∎

## References

  * [1] Alain Berlinet and Christine Thomas-Agnan.  Reproducing Kernel Hilbert Spaces in Probability and Statistics.  Springer Science, 2004.
  * [2] Léon Bottou, Frank E Curtis, and Jorge Nocedal.  Optimization methods for large-scale machine learning.  Siam Review, 60(2):223–311, 2018.
  * [3] Difeng Cai, Edmond Chow, Lucas Erlandson, Yousef Saad, and Yuanzhe Xi.  SMASH: Structured matrix approximation by separation and hierarchy.  Numerical Linear Algebra with Applications, 25, 2018.
  * [4] Michal Derezinski, Rajiv Khanna, and Michael W. Mahoney.  Improved guarantees and a multiple-descent curve for Column Subset Selection and the Nyström method.  In Advances in Neural Information Processing Systems, 2020.
  * [5] Petros Drineas and Michael W. Mahoney.  On the Nyström method for approximating a Gram matrix for improved kernel-based learning.  Journal of Machine Learning Research, 6:2153–2175, 2005.
  * [6] Dheeru Dua and Casey Graff.  UCI machine learning repository, 2019.
  * [7] Bertrand Gauthier.  Nyström approximation and reproducing kernels: embeddings, projections and squared-kernel discrepancy.  Preprint, 2021.
  * [8] Bertrand Gauthier and Luc Pronzato.  Convex relaxation for IMSE optimal design in random-field models.  Computational Statistics and Data Analysis, 113:375–394, 2017.
  * [9] Bertrand Gauthier and Johan Suykens.  Optimal quadrature-sparsification for integral operator approximation.  SIAM Journal on Scientific Computing, 40:A3636–A3674, 2018.
  * [10] Alex Gittens and Michael W. Mahoney.  Revisiting the Nyström method for improved large-scale machine learning.  Journal of Machine Learning Research, 17:1–65, 2016.
  * [11] Sanjiv Kumar, Mehryar Mohri, and Ameet Talwalkar.  Sampling methods for the Nyström method.  Journal of Machine Learning Research, 13:981–1006, 2012.
  * [12] Jason D Lee, Max Simchowitz, Michael I Jordan, and Benjamin Recht.  Gradient descent only converges to minimizers.  In Conference on learning theory, pages 1246–1257. PMLR, 2016.
  * [13] Harald Niederreiter.  Random Number Generation and Quasi-Monte Carlo Methods.  SIAM, 1992.
  * [14] Vern I. Paulsen and Mrinal Raghupathi.  An Introduction to the Theory of Reproducing Kernel Hilbert Spaces.  Cambridge University Press, 2016.
  * [15] C.E. Rasmussen and C.K.I. Williams.  Gaussian Processes for Machine Learning.  MIT press, Cambridge, MA, 2006.
  * [16] Sebastian Ruder.  An overview of gradient descent optimization algorithms.  arXiv preprint arXiv:1609.04747, 2016.
  * [17] Thomas J. Santner, Brian J. Williams, and William I. Notz.  The Design and Analysis of Computer Experiments.  Springer, 2018.
  * [18] Shusen Wang, Zhihua Zhang, and Tong Zhang.  Towards more efficient SPSD matrix approximation and CUR matrix decomposition.  Journal of Machine Learning Research, 17:7329–7377, 2016.