Revisiting Memory Efficient Kernel Approximation:
An Indefinite Learning Perspective

Simon Heilig Maximilian Münch Frank-Michael Schleif

University of Bamberg UAS Würzburg-Schweinfurt, University of Groningen UAS
Würzburg-Schweinfurt

###### Abstract

Matrix approximations are a key element in large-scale algebraic machine
learning approaches. The recently proposed method MEKA (Si et al.,, 2017)
effectively employs two common assumptions in Hilbert spaces: the low-rank
property of an inner product matrix obtained from a shift-invariant kernel
function and a data compactness hypothesis by means of an inherent block-
cluster structure. In this work, we extend MEKA to be applicable not only for
shift-invariant kernels but also for non-stationary kernels like polynomial
kernels and an extreme learning kernel. We also address in detail how to
handle non-positive semi-definite kernel functions within MEKA, either caused
by the approximation itself or by the intentional use of general kernel
functions. We present a Lanczos-based estimation of a spectrum shift to
develop a stable positive semi-definite MEKA approximation, also usable in
classical convex optimization frameworks. Furthermore, we support our findings
with theoretical considerations and a variety of experiments on synthetic and
real-world data.

##  1 INTRODUCTION

Typical machine learning algorithms contain matrix operations that do not
scale to scenarios with multiple thousand samples. Especially in the field of
kernel machines, algorithms like the well-known support vector machine (SVM),
developed by Cortes and Vapnik, (1995), contain a square kernel matrix
$\mathbf{K}\in\mathbb{R}^{n\times n}$ as the critical part, where $n$ is the
number of data points. Similar issues also exist in Gaussian Process Models
(Rasmussen and Williams,, 2006). Those matrices cause both runtime and space
complexity issues for practical applications. For this purpose, matrix
approximation has been widely studied to speed up these algorithms (Halko et
al.,, 2011; Si et al.,, 2017; Oglic and Gärtner,, 2019; Gisbrecht and
Schleif,, 2015).

An alternative to matrix approximation is to approximate a particular kernel
function using a mapping approach as in Le et al., (2013) or by an explicit
feature expansion with, e.g., random fourier features as proposed by Rahimi et
al., (2007). This, however, performed inferior in comparison to the original
MEKA. Additionally, one would like to be able to use a wide variety of kernel
functions, preferable arbitrary similarity measures, to obtain an optimal and
memory storage efficient prediction model.
As recently shown, symmetric matrices of arbitrary similarity measures can be
effectively approximated by the Nyström-method (Oglic and Gärtner,, 2019;
Gisbrecht and Schleif,, 2015).
In order to approximate a matrix using the Nyström approximation, a low-rank
structure is required to guarantee small error. The matrix obtained by the
classical Nyström method can be stored with $O(n\hat{k})$ space, where $n$ is
the number of data points and $\hat{k}$ is the total rank of the approximation
(Si et al.,, 2017). The objective of the original MEKA approach proposed by Si
et al., (2017) is to obtain a much better space complexity compared to
existing methods. After observing that clustering in the input space reveals a
block-wise structure of the matrix, Si et al., (2017) developed a novel
technique to exploit this structure without reconstructing the original kernel
matrix at any step. By this, a memory storage of $O(nk+(ck)^{2})$ for the
approximated matrix is achieved, where $k$ is the rank per clustering-block
and $c$ the number of clusters and the overall rank is $ck$ (note:
$O(n\hat{k})=O(nck)$).

The kernel theory (Hofmann et al.,, 2008) is based on positive semi-definite
(psd) and symmetric matrices, but there is recent research that argues for
also considering indefinite kernels (Huang et al.,, 2018; Mehrkanoon et al.,,
2018; Schleif and Tino,, 2015; Schleif et al.,, 2018). Such kernels emerge
from domain-specific object comparisons like protein sequence alignments or
like in case of the truncated Manhattan kernel (TL1) Huang et al., (2018) are
motivated by local learning properties. Since these kernels do not fulfill the
psd assumption of the kernel theory, new techniques were proposed to handle
such non-psd matrices (Loosli et al.,, 2015; Liu et al.,, 2018; Oglic and
Gärtner,, 2019; Gisbrecht and Schleif,, 2015) or to convert the matrices to be
psd (Schleif and Tino,, 2015).
As briefly noticed in Si et al., (2017), MEKA can lead to non-psd matrix
approximations and hence this paper aims to discuss different aspects of this
problem and to show how these limitations can be overcome.

  * •

To what extend does indefiniteness occur in MEKA and what is the effect?

  * •

How to implement a correction strategy while preserving memory efficiency?

  * •

How can MEKA be extended to handle a more extensive range of kernel functions,
especially indefinite ones?

To answer these questions, the paper is organized as follows: Section 2
introduces the relevant notation and background. Section 3 describes MEKA with
its core steps, studies the occurrence of invalid approximations, and extends
MEKA to non-stationary and indefinite kernels. Subsequently in Section 4, the
shift correction is realized by computing the shifting parameter using the
Lanczos iteration in a highly efficient manner. Our proposal is evaluated with
different classification tasks solved by an SVM using different kernel
functions.

##  2 PRELIMINARIES

This section introduces basic notation and concepts of matrix approximation
and kernel methods.

###  2.1 Matrix Approximation

Given a matrix $\mathbf{A}\in\mathbb{R}^{n\times n}$, the objective of matrix
approximation algorithms is to achieve a low approximation error at a given
memory constraint, specified by the rank of the resulting matrix
$\mathbf{\tilde{A}}$:

$\min{\|\mathbf{A}-\mathbf{\tilde{A}}\|_{F}}$ (1)

where $\mathbf{\tilde{A}}$ is the approximated matrix with rank $k$ and
$\|\cdot\|_{F}$ the Frobenius norm (Kumar and Schneider,, 2017). This concept
can be seen as a data compression technique, since an intrinsic low-rank
matrix contains many linear dependent rows or columns, reflecting redundant
information. Note that the optimal solution for Equation (1) is the truncated
singular value decomposition (SVD), which is prohibitive in the large-scale
setting. In order to maintain efficiency with huge amounts of data,
probability-based decomposition methods such as randomised CUR decomposition
have been developed (Drineas et al.,, 2008). The idea is to select certain
rows of $\mathbf{A}$ in matrix $\mathbf{R}$ and columns in matrix $\mathbf{C}$
to capture the relevant subspace, with $\mathbf{U}$ as a link matrix and
eventually solve an optimization problem:

$\min\;||\mathbf{A}-\mathbf{C}\mathbf{U}\mathbf{R}||^{2}_{F}.$ (2)

In particular for kernel matrices $\mathbf{K}$ (square and symmetric),
according to Williams and Seeger, (2001) and Sun et al., (2015), it is
possible to approximate the eigenspectrum with the Nyström expansion, where
only $l$ landmarks (selected columns or rows) and a pseudo-inverse of a
rank-$k$ approximation of the $l\times l$ intersection matrix $\mathbf{W}$ are
required to build a rank-$k$ approximation of the original matrix:

$\mathbf{\tilde{K}}=\mathbf{K}_{(\cdot,\mathcal{I})}\mathbf{W}_{k}^{\dagger}\mathbf{K}_{(\cdot,\mathcal{I})}^{T},$
(3)

where $\mathcal{I}$ is a set of landmark indicators consisting of $l$ indices
and $\mathbf{W}_{k}^{\dagger}$ as the pseudo-inverse of the rank-$k$
approximation obtained by the truncated SVD. As it is later on required in the
MEKA formulation to decompose a matrix not in three but in two matrices, one
can split the $\mathbf{W}_{k}^{\dagger}$ matrix based on the eigenvalue
decomposition (EVD):
$\mathbf{W}_{k}=\mathbf{U}_{W,k}\mathbf{\Lambda}_{W,k}\mathbf{U}_{W,k}^{T}$.
Assuming only non-negative eigenvalues, Equation (3) is written as follows:

$\displaystyle\begin{split}\mathbf{\tilde{K}}&=\underbrace{\mathbf{K}_{(\cdot,\mathcal{I})}\mathbf{U}_{W,k}\left(\mathbf{\Lambda}_{W,k}^{\dagger}\right)^{\frac{1}{2}}}_{\mathbf{H}}\underbrace{\left(\mathbf{\Lambda}_{W,k}^{\dagger}\right)^{\frac{1}{2}}\mathbf{U}_{W,k}^{T}\mathbf{K}_{(\cdot,\mathcal{I})}^{T}}_{\mathbf{H}^{T}},\\\
&=\mathbf{H}\mathbf{H}^{T}.\end{split}$ (4)

###  2.2 Kernel Theory

Given a collection of $n$ data points
$X=\\{\mathbf{x}_{1},\mathbf{x}_{2},\cdots,\mathbf{x}_{n}\\},\text{ where
}\mathbf{x}_{i}\in\mathcal{X}$, the kernel matrix
$\mathbf{K}_{i,j}=k(\mathbf{x}_{i},\mathbf{x}_{j})\in\mathbb{R}^{n\times n}$
is created with a kernel function
$k:\mathcal{X}\times\mathcal{X}\to\mathbb{R}$, where
$k(\mathbf{x},\mathbf{y})=\langle\phi(\mathbf{x}),\phi(\mathbf{y})\rangle_{\mathcal{H}}$
with $\phi:\mathcal{X}\to\mathcal{H}$ as a mapping into an arbitrarily high
dimensional Hilbert space $\mathcal{H}$. Following the line of Si et al.,
(2017), we assume $\mathcal{X}=\mathbb{R}^{d}$ as $d$-dimensional input space.

As a matter of fact, any function of the form of $k$ which gives rise to a psd
kernel matrix reflects such an implicit inner product in a possibly infinite-
dimensional Hilbert space (Hofmann et al.,, 2008). This technique allows a
flexible way of defining a data representation leading to a powerful and more
important linear classifier in such a Hilbert space. By employing the inner
product of this space implicitly, without the need of the explicit mapping
$\phi$ to be present, a non-linear classification in the input space
$\mathcal{X}$ is achieved, also known as the kernel trick (Hofmann et al.,,
2008).

If a kernel function leads to a non-psd kernel matrix, a Kreĭn space is
induced instead (Duin and Pekalska,, 2005; Schleif and Tino,, 2015). A Kreĭn
space $\mathcal{K}$ is spanned by two orthogonal Hilbert spaces
$\mathcal{H}_{+}$ and $\mathcal{H}_{-}$ so that
$\langle\mathbf{x},\mathbf{y}\rangle_{\mathcal{K}}=\langle\mathbf{x}_{+},\mathbf{y}_{+}\rangle_{\mathcal{H}_{+}}-\langle\mathbf{x}_{-},\mathbf{y}_{-}\rangle_{\mathcal{H}_{-}}$
for $\mathbf{x},\mathbf{y}\in\mathcal{K}$ and $\mathbf{x}_{+},\mathbf{y}_{+}$
as the projection onto $\mathcal{H}_{+}$ and $\mathbf{x}_{-},\mathbf{y}_{-}$
as the projection onto $\mathcal{H}_{-}$.
Consequently, the norm induced by the indefinite inner product lacks of metric
properties, which is a problem for branch-and-bound methods (Morrison et al.,,
2016) and the objective function of a SVM becomes non-convex (Schleif and
Tino,, 2015). For kernel approaches, this problem can be tackled by various
strategies, e.g., learning a proxy matrix which is constrained to be psd (Luss
and d’Aspremont,, 2009), by correcting the eigenspectrum of the matrix to
remove the negative eigenvalues (Schleif and Tino,, 2015; Münch et al.,, 2020)
or by developing methods in the Kreĭn space, e.g., the indefinite SVM from
Loosli et al., (2015); Schleif et al., (2018); Schleif and Tiño, (2017) or
Oglic and Gärtner, (2019); Gisbrecht and Schleif, (2015) and by a shift
correction based SVM from Loosli, (2019). These approaches are in general
costly and do not scale to large data as considered by MEKA, but only to a few
thousand samples at most.

A kernel function is shift-invariant if and only if
$k(\mathbf{x},\mathbf{y})=f(\eta(\mathbf{x}-\mathbf{y}))$ for a function
$f:\mathcal{X}\to\mathbb{R}$ and a scaling parameter $\eta$. As a consequence,
kernel functions based on distances, such as the Gaussian rbf kernel:

$k_{rbf}(\mathbf{x},\mathbf{y})=\exp\left(-\gamma\|\mathbf{x}-\mathbf{y}\|_{2}^{2}\right),$

are shift-invariant. However, for inner product based (non-stationary)
functions, a restriction to the unit sphere has to be made (Pennington et
al.,, 2015).

######  Proposition 1 (Pennington et al., (2015)).

If the input data is restricted to a unit sphere via the norm $\|\cdot\|_{2}$,
then any function $f:\mathcal{X}\times\mathcal{X}\to\mathbb{R}$ that is based
on a computation of an inner product of the input data is shift-invariant.

###### Proof.

Given $\mathbf{x},\mathbf{y}\in\mathcal{X}$ in $d$ dimensions, so that
$\mathbf{x}=(x_{1},x_{2},\dots,x_{d-1},x_{d})$ and
$\mathbf{y}=(y_{1},y_{2},\dots,y_{d-1},y_{d})$.

$\displaystyle\|\mathbf{x}-\mathbf{y}\|_{2}^{2}=\langle\mathbf{x}-\mathbf{y},\mathbf{x}-\mathbf{y}\rangle_{2}$
$\displaystyle=\sum_{i=1}^{d}{(x_{i}-y_{i})^{2}}=\sum_{i=1}^{d}{(x_{i}^{2}-2x_{i}y_{i}+y_{i}^{2})}$
$\displaystyle=\underbrace{\sum_{i=1}^{d}{x_{i}^{2}}}_{=1}+\sum_{i=1}^{d}{(-2x_{i}y_{i})}+\underbrace{\sum_{i=1}^{d}{y_{i}^{2}}}_{=1}\;\;,\|\mathbf{x}\|_{2}=\|\mathbf{y}\|_{2}=1$
$\displaystyle=2-2\sum_{i=1}^{d}{x_{i}y_{i}}=2-2\langle\mathbf{x},\mathbf{y}\rangle_{2}\qed$

##### Remark 1.

In case there is no real vector space given and only a kernel matrix is
present, a matrix normalization is required. Since the diagonal elements
represent the kernel function evaluation of
$k(\mathbf{x}_{i},\mathbf{x}_{i}),\;i=1\dots n$, restricting the input space
on an unit sphere is equivalent to a constant diagonal. This can be achieved
by:
$\langle\mathbf{x}_{i},\mathbf{x}_{j}\rangle=\frac{\langle\mathbf{x}_{i},\mathbf{x}_{j}\rangle}{\sqrt{\langle\mathbf{x}_{i},\mathbf{x}_{i}\rangle}\sqrt{\langle\mathbf{x}_{j},\mathbf{x}_{j}\rangle}}$,
where
$\langle\mathbf{x}_{i},\mathbf{x}_{i}\rangle>0,\langle\mathbf{x}_{j},\mathbf{x}_{j}\rangle>0,\;i,j=1\dots
n$. Moreover, by this transformation, the self-similarities of any kernel
become a constant.

##### Remark 2.

As noted in Pennington et al., (2015), there are kernel functions that become
non-psd in consequence of this normalization since it has been proven that
specific kernel functions are not a Fourier transform of a finite non-negative
Borel measure anymore. As an example, the polynomial kernel:

$k_{poly}(\mathbf{x},\mathbf{y})=\left(1-\frac{\|\mathbf{x}-\mathbf{y}\|^{2}_{2}}{a^{2}}\right)^{p}=\alpha(q+\langle\mathbf{x},\mathbf{y}\rangle)^{p},$

with $q=(a^{2}/2)-1$ and $\alpha=(a^{2}/2)^{p}$ becomes non-psd for $a>2$ and
$p>=1$. Hence, non-psd kernel functions occur more frequently than expected
and need to be addressed (Pennington et al.,, 2015; Liu et al.,, 2021).

##### Remark 3.

The projection of the input space onto a sphere can lead to overlapping of
clusters, which were on a similar projection ray but had a distinctive
separation beforehand. Further, we will assume that this effect is negligible,
since it is a data specific use case and multiple works follow the line of
Pennington et al., (2015) such as Liu et al., (2021), which extends also to
dot-product based kernel functions. A learned embedding on such a manifold
while preserving the differences is able to overcome this problem, as
mentioned in (Wilson and Hancock,, 2010), but it needs to fit in the large-
scale setting of this work.

##  3 ANALYSIS OF MEKA

(a)

(b)

(c)

(d)

Figure 1: Illustration of the block-wise structure in the kernel matrix
associated to a Gaussian rbf kernel, applied on the pendigit data set. The
darker the higher the value of the kernel function. (a) the unordered kernel
matrix for $\gamma=0.05$. (b) rearranged matrix with respect to a c-means
clustering. (c) $\gamma=1\cdot 10^{-6}$ leads to the low-rank case. (d)
$\gamma=1.2$ leads to a block-diagonal kernel matrix and finally approaches an
identity matrix.

After introducing the basic foundations of kernels and matrix approximations
thereof, a brief outline of the MEKA algorithm is presented below. Further, we
enable MEKA to be applicable for non-stationary and indefinite kernels and the
occurrence of both positive and negative eigenvalues is studied on real-world
and synthetic data.

###  3.1 Memory Efficient Kernel Approximation

The key aspect of MEKA is the combination of a clustering in the input space
exploiting the effect of the scaling parameter of a shift-invariant kernel
function. As Si et al., (2017) observed, a kernel matrix arranged with respect
to a clustering changes from low-rank to block-diagonal when the scaling
parameter $\gamma$ increases. Figure 1 (a) shows a kernel matrix created
without any ordered structure. A rearranged version in Figure 1 (b) clearly
shows a block-diagonal structure. The influence of the scaling parameter can
be seen in Figure 1 (c) - (d), where a small scaling parameter leads to a low-
rank structure, whereas a large scaling parameter causes a block-diagonal
matrix. Exploiting this structure, a highly memory-efficient representation of
the matrix is generated by MEKA, which is capable of handling both the low-
rank and the block-wise scenario. The algorithm is structured as follows:

##### Part 1.

Applying the c-means algorithm on a random subset of $X$ to compute $c$
cluster centers and assign the remaining data points to the nearest cluster.
This results in $c$ cluster centroids
$\mathbf{c}_{1},\mathbf{c}_{2},\dots,\mathbf{c}_{c}$ with an associated set of
indices per cluster $\mathcal{C}_{1},\mathcal{C}_{2},\dots,\mathcal{C}_{c}$
leading to $c$ disjoint sets $V_{1},V_{2},\dots,V_{c}$, where
$V_{i}=\\{\mathbf{x}_{j}\in\mathbb{R}^{d}\;|\;j\in\mathcal{C}_{i}\\}$ and
$V_{i}\cap V_{j}=\emptyset$ for $i\neq j$, hence $X=\bigcup_{i=1}^{c}V_{i}$.

##### Part 2.

Each block on the diagonal
$\mathbf{K}^{i}=k(\mathbf{x}_{p},\mathbf{x}_{q})\;\forall\mathbf{x}_{p},\mathbf{x}_{q}\in
V_{i}$ has a low-rank structure as shown by Si et al., (2017), so a Nyström
approximation is conducted to build a low-rank approximation
$\mathbf{K}^{i}\approx\mathbf{Q}^{i}\left(\mathbf{Q}^{i}\right)^{T}$ according
to Equation (4). The block-diagonal of $\mathbf{K}$ is represented as
$\mathbf{K}_{blk\,diag}\approx\mathbf{Q}\mathbf{Q}^{T},\;\text{with
}\mathbf{Q}=\bigoplus_{i=1}^{c}\mathbf{Q}^{i}$, using the notation of Si et
al., (2017), this provides a block-kernel approximation with zeros off-
diagonal blocks. Si et al., (2017) proposed several strategies to determine a
proper rank per block: 1) Using the same rank $k$ in each Nyström
approximation. This results in a fast algorithm but could lead to a high
approximation error. In general, the overall target rank is $ck$ of the MEKA
approximation. 2) Distribute the rank according to singular value
contribution. Each block gets a rank $k_{i}$ such that
$\sum_{i=1}^{c}k_{i}=ck$ depending on the singular value contribution to the
first $ck$ singular values of $\mathbf{K}_{blk\,diag}$. This strategy
decreases MEKA’s speed but results in a better approximation error. The
overall target rank is $ck$. 3) Using a rank $k_{i}$ proportional to the ratio
of the cluster size and the number of all samples $n$, namely
$k_{i}=\frac{|V_{i}|}{n}k$ so that $\sum_{i=1}^{c}k_{i}=k$. This results in an
overall rank of $k$. In accordance with Si et al., (2017), we follow this
strategy.

##### Part 3.

Since many relevant off-diagonal blocks are present in the low-rank cases, Si
et al., (2017) suggested a highly efficient least squares approximation of
such blocks
$\mathbf{K}^{i,j}=k(\mathbf{x}_{p},\mathbf{x}_{q})\;\forall\mathbf{x}_{p}\in
V_{i},\forall\mathbf{x}_{q}\in V_{j}$. They proposed to subsample the original
block $\mathbf{K}^{i,j}$ with random row and column indices sets
$\mathcal{I}_{i},\mathcal{I}_{j}$ to get
$\mathbf{\bar{K}}^{i,j}=\mathbf{K}_{\mathcal{I}_{i},\mathcal{I}_{i}}$ and to
solve:

$\underset{\mathbf{L}^{i,j}}{\operatorname{arg}\,\operatorname{min}}\;\left\|\mathbf{K}^{i,j}-\mathbf{Q}^{i}\mathbf{L}^{i,j}\left(\mathbf{Q}^{j}\right)^{T}\right\|_{F},\;\text{where
}i\neq j,$

with a closed form solution:

$\mathbf{L}^{i,j}=((\mathbf{Q}^{i}_{\mathcal{I}_{i}})^{T}\mathbf{Q}^{i}_{\mathcal{I}_{i}})^{-1}(\mathbf{Q}^{i}_{\mathcal{I}_{i}})^{T}\mathbf{\bar{K}}^{i,j}\mathbf{Q}^{j}_{\mathcal{I}_{j}}((\mathbf{Q}^{j}_{\mathcal{I}_{j}})^{T}\mathbf{Q}^{j}_{\mathcal{I}_{j}})^{-1}$
(5)

to finally approximate $\mathbf{K}$ by
$\mathbf{\tilde{K}}=\mathbf{Q}\mathbf{L}\mathbf{Q}^{T}$. Algorithm 1 provides
a detailed procedure of MEKA.

Algorithm 1 Memory Efficient Kernel Approximation (MEKA), adapted from Si et
al., (2017)

Data points in $X$, $\mathcal{X}=\mathbb{R}^{\times d}$, kernel parameter
$\gamma$, target rank $k$, number of clusters $c$ and percentage of blocks to
be truncated

Matrices $\mathbf{Q}\in\mathbb{R}^{n\times
ck},\mathbf{L}\in\mathbb{R}^{ck\times ck}$

Run c-means to obtain $V_{1},V_{2},\dots,V_{c}$$\triangleright$ Part 1

for $i\leftarrow 1\;\text{to}\;c$ do$\triangleright$ Part 2

Approximate $\mathbf{K}^{i}$ with rank $k_{i}$ by Equation (4)

$\mathbf{Q}\leftarrow\mathbf{Q}\oplus\mathbf{Q}^{i}$

$\mathbf{L}\leftarrow\mathbf{L}\oplus\mathbf{I}\in\mathbb{R}^{k\times k}$

end for

for all $(i,j)(i\neq j)$ do $\triangleright$ Part 3

if Clusters too dissimilar then

$\mathbf{L}^{i,j}\leftarrow\mathbf{0}$

Continue

end if

Obtain $\mathbf{L}^{i,j}$ by Equation (5)

end for

(a)

(b)

(c)

Figure 2: MEKA experiments with pendigit data set with $c=10$ and Gaussian rbf
kernel. (a) illustrates the relative approximation error according to the
target rank and the scaling parameter. Note that after $\gamma>0.1$, the
matrix tends approach the full rank, hence, the approximation error worsens
and is limited to 100%. (b) shows the value of the largest negative
eigenvalue. (c) presents the total number of negative eigenvalues.

###  3.2 Non-Stationary and Indefinite Kernels

Si et al., (2017) assumed the input data to originate from a real vector space
and accessed the Hilbert space via a kernel function. Furthermore, their
theorems and error-bounds are based on the concept of shift-invariant
functions. In order to ensure a wider applicability of MEKA, an extension to
non-rbf like kernels is necessary. By applying Proposition 1 to the input
space, non-stationary kernel functions become shift-invariant. While MEKA
considers the effect of the Gaussian rbf scaling parameter in the light of a
shift-invariant kernel function, Proposition 1 allows a significant increase
of freedom to select a kernel function 111Further visualisations of non-
stationary matrices with a clear block-wise structure are given in the
supplements.. However, as mentioned above, kernels such as the polynomial,
arc-cosine or neural tangent kernel can become indefinite by applying this
normalization (Remark 1) to the input space or the kernel function (Pennington
et al.,, 2015; Liu et al.,, 2021). This leads to the question of how to
approximate indefinite kernels with MEKA.
Therefore, the crucial step of MEKA is the approximation of each diagonal
block computed by Nyström. This can be accomplished due to Oglic and Gärtner,
(2019), who showed that Nyström methods are able to handle both psd and
indefinite kernel matrices within a bounded error.

###  3.3 Occurrence of Negative Eigenvalues

|  |
---|---|---
(a) |  (b) |  (c)
Figure 3: MEKA experiments with gesture data set and extreme learning kernel.
(a) shows the relative approximation error along the different target ranks
and numbers of clusters. (b) presents the value of the smallest eigenvalue and
(c) shows the total number of negative eigenvalues.

As Si et al., (2017) mentioned, the subsampling strategy of the off-diagonal
block approximation (Part 3 of MEKA) potentially results in a non-psd
approximation. In the following, we critically discuss the importance of
MEKA’s hyperparameters with respect to indefiniteness. Following Si et al.,
(2017), we first focus on the Gaussian rbf kernel, and then we apply
Proposition 1 to show the approximation capabilities on inner product based
kernels, in particular, the extreme learning kernel (ELM):

$k_{elm}(\mathbf{x},\mathbf{y})=\frac{2}{\pi}\arcsin\left(\frac{1+\langle\mathbf{x},\mathbf{y}\rangle}{\sqrt{(\frac{1}{2\sigma^{2}_{w}}+1+\langle\mathbf{x},\mathbf{x}\rangle)(\frac{1}{2\sigma^{2}_{w}}+1+\langle\mathbf{y},\mathbf{y}\rangle)}}\right)$

(6)

from Frénay and Verleysen, (2011), as a parameter-insensitive Gaussian rbf
kernel alternative independent of $\sigma_{w}$.

Table 1: Overview of data set characteristics. Data set | ID |  Samples $n$ |  Attributes $d$
---|---|---|---
spambase | 1 | 4,601 | 57
artificial 1 | 2 | 7,500 | 2
cpusmall | 3 | 8,192 | 12
gesture | 4 | 9,873 | 32
artificial 2 | 5 | 10,000 | 15
pendigit | 6 | 10,992 | 16

In Table 1, the utilized data sets are listed. We retrieved the gesture data
set from Madeo et al., (2013) and spambase from the UCI
repository222https://archive.ics.uci.edu/ml/index.php and cpusmall and
pendigit are obtained from the LIBSVM
repository333https://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/. We also
included two artificial data sets created with scikit-learn to demonstrate
performance with overlapping clusters. Within the experiments, we collect the
following statistics: relative approximation error
$\|\mathbf{K}-\mathbf{\tilde{K}}\|_{F}/\|\mathbf{K}\|_{F}$, the smallest
eigenvalue and the number of negative eigenvalues. Note that for
reproducibility, a random generator seed has been set and the calculations
were executed with Matlab R2020a on a machine with
Intel${}^{\text{\textregistered}}$ Xeon W-2175 CPU @ 2.50GHz $\times$ 28\. The
code of the extended MEKA implementation and the experimental setups of all
experiments in this work can be found at:
https://github.com/simonheilig/indefinite-meka.

Figure 2 presents an analysis of MEKA with respect to the scaling parameter
$\gamma$ of the Gaussian rbf kernel and the target rank of the approximation.
To follow the line of Si et al., (2017), the rank distribution strategy
relative to the cluster size is utilized in this experiment. It is apparent
that a substantial negative eigenspace is present for $10^{-6}\leq\gamma\leq
1$. For $\gamma>1$, the matrix is close to full rank, resulting in an
increasing approximation error. Figure 2 (c) shows a direct correlection
between the target rank of the approximation and the number of negative
eigenvalues.
Note that the results for the remaining data sets are given in the
supplements, because they show similar behavior except for small changes.

Figure 3 shows the results for the gesture data set and the extreme learning
kernel. Since the kernel performance is invariant to its parameter (see
Equation (6)), the development along with the target rank and the number of
clusters is presented. Figure 3 (a) shows the connection between the relative
approximation error and the target rank of the approximation.
Obviously, an increasing target rank results in a lower approximation error.
Furthermore, an increasing number of clusters leads to higher approximation
errors, since more off-diagonal blocks are present in those cases. Figure 3
(b) reflects a similar behavior, as more off-diagonal blocks require many
subsampling steps in part 3 of MEKA. Finally, Figure 3 (c) shows a positive
correlation between the target rank and the number of negative eigenvalues.

##  4 REVISION

From the analyses above, it becomes evident that a psd-correction is required
to apply MEKA to kernel methods in a variety of optimization frameworks. This
is necessary due to psd violations in the original MEKA approximation, but
also if more general kernel functions are used, which are non-psd. Si et al.,
(2017) mentioned briefly that psd violations may occur but claimed it as a
rare event. They proposed to clip the link matrix $\mathbf{L}$, which leads to
the closest psd version of the approximation on the cone of all psd matrices
under the Frobenius norm (Schleif and Tino,, 2015). Since this matrix is
$ck\times ck$, it is prohibitive to obtain an eigenvalue decomposition if
there is a large number of clusters, but if the negative eigenvalues are due
to noise, clip is the method of choice. As shown previously, indefinite
contributions may occur regularly, but we cannot yet ensure whether this is
exclusively due to noise or not. The second approach, suggested by Si et al.,
(2017), is to apply an appropriate shift on the diagonal of the entire
approximation without specifying how to determine a shift parameter while
maintaining the time and memory efficiency. In particular, if the indefinite
part contains relevant information, a shift correction is a better choice
(Duin and Pekalska,, 2010; Münch et al.,, 2020).
Given a computed shift parameter $\lambda_{shift}$ and the smallest eigenvalue
$\lambda_{min}$ of $\mathbf{\tilde{K}}$, the factorized MEKA approximation can
be efficiently corrected by
$\mathbf{\tilde{K}}=\mathbf{Q}\mathbf{L}\mathbf{Q}^{T}+\lambda_{shift}\mathbf{I}$,
where $\lambda_{shift}\geq|\lambda_{min}|$ and $\mathbf{I}$ as identity
matrix, since the shift operation on the diagonal of the similarity matrix
causes an analogous shift of the eigenspectrum of the respective matrix (Laub
et al.,, 2006; Schleif and Tino,, 2015; Loosli,, 2019; Münch et al.,, 2020).
Note that it is prohibitive to determine this shift parameter by an eigenvalue
decomposition, as it violates the runtime and memory complexity constraint.
The out-of-sample extension is discussed in the supplements.

A well-known method to calculate only a few of the eigenvalues, in particular
the extreme values, is given by the Lanczos iteration in a very efficient
manner (Cullum and Willoughby,, 2002). The essential characteristic of this
approach is that it works with a matrix-times-vector multiplication, since it
would be prohibitive to build up kernel matrices entirely in the large-scale
setting. The factorization obtained with MEKA enables a matrix-times-vector
multiplication $\mathbf{y}=\mathbf{\tilde{K}}\mathbf{x}$ within a time
complexity proportional to the memory requirement $O(nk+(ck)^{2})$ (Si et
al.,, 2017).
See Table 2 for an overview of the complexities of the correction strategies,
where $n$ is the original number of rows or columns, $m$ the remaining number
of eigenvalues and eigenvectors after truncation, $j$ the number of Lanczos
iterations (usually $j\approx 2\sqrt{n}$ is enough (Parlett et al.,, 1982))
and $T$ as the complexity of matrix-times-vector multiplication (for MEKA
$T=O(j(nk+(ck)^{2})+j\ln{j})$) (Mahadevan,, 2008; Coakley and Rokhlin,, 2013).
The advantage of a Lanczos-based shift is to allow any matrix representation,
as long as the matrix-times-vector multiplication can be guaranteed. By
utilizing the MEKA factorization of the kernel matrix, it is unnecessary to
keep the entire matrix in memory, thus maintaining memory efficiency.

Table 2: Overview of cost complexities of different correction methods and their storage requirements (n.d. = not defined).  Type | Clip | Shift | Lanczos-Shift
---|---|---|---
Runtime | $O(n^{3})$ | $O(n^{3})$ | $O(jT+j\ln{j})$
Storage | $O(nm+m^{2})$ | $O(n^{2})$ | n.d.

###  4.1 Lanczos Iteration

The Lanczos iteration aims to minimize the following optimization problem
(also known as the Rayleigh quotient): $\lambda_{min}=\min_{\mathbf{x}\neq
0}\frac{\langle\mathbf{x},\,\mathbf{A}\mathbf{x}\rangle}{\langle\mathbf{x},\,\mathbf{x}\rangle},$
with a symmetric matrix $\mathbf{A}\in\mathbb{R}^{n\times n}$ and its smallest
eigenvalue $\lambda_{min}$ (Cullum and Willoughby,, 2002). Note that
situations occur where the algorithm needs to be restarted or a
reorthogonalization is required. However, there exist approaches without the
need for such reorthogonalization, by identifying the eigenvalues which lead
to the linear dependencies in the Lanczos vectors and ignoring them (Cullum
and Willoughby,, 2002). This is particularly important in our scenario to keep
a good runtime performance. The convergence rate is heavily related to the
overall stiffness of the matrix, but usually, the extreme Ritz values are
obtained after $i\approx 2\sqrt{n}$ iterations, leading to good approximations
of the corresponding eigenvalues of $\mathbf{A}$ (Parlett et al.,, 1982).

###  4.2 Error Bound Analysis

Table 3: Classification results using a SVM with L1QP solver on the original and unapproximated matrix (M0), approximated matrix before (M1) and after correction (M2) by a 10-fold CV accuracy and standard deviation in % (n.c. $=$ not converged). Typical memory costs of the approximation in (Mem). Experiments for M0∗ are conducted on a shift corrected matrix.  | Shift Invariant | Non-Stationary | Indefinite
---|---|---|---
| RBF | ELM | Poly | TL1
ID | Mem | M0 | M1 | M2 | M0 | M1 | M2 | M0 | M1 | M2 |  M0* |  | M2
1 | 4.14 % | 93.63 | n.c. | 88.85 | 95.09 | n.c. | 81.81 | 91.48 | n.c. | 73.66 | 94.41 |  | 85.57
($1.27$) | ($4.32$) | ($0.81$) | ($5.17$) | ($1.92$) | ($7.40$) | ($1.29$) |  | ($3.74$)
2 | 1.11 % | 89.19 | 85.48 | 89.23 | 88.60 | 82.20 | 88.63 | 82.43 | n.c. | 82.49 | 78.96 |  | 89.31
($1.34$) | ($11.93$) | ($1.42$) | ($1.32$) | ($17.63$) | ($1.37$) | ($1.15$) | ($1.03$) | ($13.32$) |  | ($0.50$)
3 | 0.30 % | 89.42 | n.c. | 86.47 | 89.01 | n.c. | 86.93 | 82.46 | n.c. | 77.27 | 89.22 |  | 88.53
($0.91$) | ($1.32$) | ($1.15$) | ($1.29$) | ($1.81$) | ($1.76$) | ($1.26$) |  | ($1.48$)
4 | 0.18 % | 58.99 | n.c. | 44.13 | 68.02 | n.c. | 45.50 | 64.18 | n.c. | 38.57 | 54.78 |  | 45.56
($0.95$) | ($1.13$) | ($2.45$) | ($1.23$) | ($1.75$) | ($1.53$) | ($1.51$) | ($0.95$)
5 | 1.75 % | 88.79 | 72.77 | 88.76 | 88.74 | n.c. | 85.71 | 88.79 | 88.79 | 88.79 | 88.66 |  | 88.61
($1.08$) | ($24.32$) | ($0.98$) | ($1.35$) | ($1.87$) | ($0.95$) | ($0.95$) | ($0.95$) | ($0.83$) |  | ($1.00$)
6 | 2.18 % | 99.58 | 21.04 | 87.79 | 99.52 | 10.37 | 83.26 | 99.52 | 39.23 | 98.01 | 99.18 |  | 96.53
($0.10$) | ($12.88$) | ($2.54$) | ($0.18$) | ($0.00$) | ($1.96$) | ($0.29$) | ($32.40$) | ($0.56$) | ($0.26$) |  | ($0.71$)

Given the error bound of Theorem 3 of Si et al., (2017):

$\displaystyle\|\mathbf{K}-\mathbf{\tilde{K}}\|_{F}\leq\|\mathbf{K}-\mathbf{K}_{k}\|_{F}+\left(\frac{64k}{l}\right)^{\frac{1}{4}}n\mathbf{K}_{max}(1+\theta)^{\frac{1}{2}}$
$\displaystyle\qquad\qquad\qquad+2\|\Delta\|_{F}$

where
$\theta=\sqrt{\frac{n-l}{n-0.5}\frac{1}{\beta(l,n)}\log{\frac{1}{\delta}}}d_{max}^{\mathbf{K}}/\mathbf{K}^{\frac{1}{2}}_{max}$;
$\beta(l,n)=1-\frac{1}{2\max\\{l,n-l\\}}$;
$\mathbf{K}_{max}=\max_{i}\mathbf{K}_{(i,i)}$; and $d_{max}^{\mathbf{K}}$
represents the distance
$max_{ij}\sqrt{\mathbf{K}_{(i,i)}+\mathbf{K}_{(j,j)}-2\mathbf{K}_{(i,j)}}$.
The error bound for the corrected MEKA (assuming that
$\mathbf{K},\mathbf{\tilde{K}}$ are non-psd) is:

$\displaystyle\|(\mathbf{K}-\mathbf{\tilde{K}})+\lambda_{shift}\mathbf{I}\|_{F}$
$\displaystyle=\|(\mathbf{K}^{+}-\mathbf{\tilde{K}}^{+})+(-\mathbf{K}^{-}+\mathbf{\tilde{K}}^{-})+\lambda_{shift}\mathbf{I}\|_{F}$
$\displaystyle\leq\|\mathbf{K}^{+}-\mathbf{\tilde{K}}^{+}\|_{F}+\|-\mathbf{K}^{-}+\mathbf{\tilde{K}}^{-}\|_{F}+\|\lambda_{shift}\mathbf{I}\|_{F}$
$\displaystyle=\|\mathbf{K}^{+}-\mathbf{\tilde{K}}^{+}\|_{F}+\|\mathbf{K}^{-}-\mathbf{\tilde{K}}^{-}\|_{F}+\|\lambda_{shift}\mathbf{I}\|_{F}$
$\displaystyle\leq\|\mathbf{K}^{+}-\mathbf{K}^{+}_{k}\|_{F}+\left(\frac{64k}{l}\right)^{\frac{1}{4}}n\mathbf{K}^{+}_{max}(1+\theta)^{\frac{1}{2}}+2\|\Delta_{+}\|_{F}$
$\displaystyle+\|\mathbf{K}^{-}-\mathbf{K}^{-}_{k}\|_{F}+\left(\frac{64k}{l}\right)^{\frac{1}{4}}n\mathbf{K}^{-}_{max}(1+\theta)^{\frac{1}{2}}+2\|\Delta_{-}\|_{F}$
$\displaystyle+\sqrt{n}|\lambda_{shift}|$

where $\mathbf{K}=\mathbf{K}^{+}-\mathbf{K}^{-}$ and
$\mathbf{\tilde{K}}=\mathbf{\tilde{K}}^{+}-\mathbf{\tilde{K}}^{-}$ are
provided by the definition of Kreĭn spaces, which can be decomposed into two
orthogonal Hilbert spaces. Apparently, the shift parameter linearly affects
the error bound and it is important to estimate a reliable $\lambda_{shift}$
as proposed in Section 4.1. The value of the smallest eigenvalue is bounded by
the row sums of the matrix, as given by Gerschgorin disks (Gershgorin,, 1931).

###  4.3 Classification Experiments

To study the effect of the Lanczos-shift, we compare our proposed correction
on a kernel SVM with a classical QP solver. Four kernel functions are
considered: the Gaussian rbf kernel, polynomial kernel, extreme learning
kernel, and an instance of an indefinite kernel, the shift invariant truncated
Manhattan kernel (TL1):

$k_{tl1}(\mathbf{x},\mathbf{y})=\max\\{\rho-\|\mathbf{x}-\mathbf{y}\|_{1},0\\}$

The TL1 kernel is successfully used for local training and is less sensitive
to its hyperparameter (Huang et al.,, 2018; Mehrkanoon et al.,, 2018).
According to Huang et al., (2018), $\rho=0.7d$ is superior to some
hyperparameter optimized Gaussian rbf kernel results. The other kernel
parameters are also optimized using a grid search
$\gamma\in\\{1,2,\ldots,10,15\\}$ and $a\in\\{2,3\\}$, $p\in\\{2,4,8,10\\}$.
The $c,k$ parameters of MEKA are determined to have the lowest approximation
error with minimal memory usage and the soft-margin parameter $C$ of the SVM
is chosen to minimize the classification loss of the original kernel SVM using
a 5-fold cross-validation (CV) in a grid of $C\in\\{10^{-3},1,10,100,1000\\}$.
For comparison purposes, the rank distribution strategy with the same rank per
cluster is used. Further details on preprocessing, typical hyperparameters and
a brief summary of utilized kernels is given in the supplements. For the TL1
kernel, the original kernel matrix is non-psd and thus the parameter
optimizations are executed on a shift corrected matrix fulfilling the psd
constraint.
In Table 3, the classification accuracy of a 10-fold CV (with a nested 5-fold
CV for parameter estimation) of our Lanczos-shift correction is compared to
the original kernel SVM and the original MEKA approximation. Also, the typical
memory consumption of the factorization provided by MEKA is given relative to
the original matrix. The approximation quality can be found in the
supplementary files.

The results show that the frequent negative eigenvalues (see Sec. 3.3), result
in common optimization problems, occurring with indefinite kernels. The shift
correction effectively solved this issue. Comparing the psd results M0 (exact)
for {RBF,ELM,Poly} with M2 (approximated), we see that the proposal has often
a similar accuracy. In the rare cases where classical MEKA (M1) converged, the
proposed solution (M2) is consistently better (one case on par). For the
indefinite TL1 we always used a shift correction such that no M1 exists. Here
we see competitive results of the exact and corrected M0∗ compared with the
approximated and corrected M2. As expected the memory consumption using M1 or
M2 can be kept very small, with only a few percent of the original space and
minimum information loss for M2, whereas M1 (original MEKA) causes convergence
problems due to the non-psd effect. In particular, the results show that MEKA
is now also effectively applicable to more generic kernel functions using the
proposal with the underlying normalization procedure.

##  5 CONCLUSION

We showed that MEKA can be effectively extended to a variety of common kernel
functions. In particular, we showed that a spherical projection extends the
applicability to non-stationary kernels. Moreover, we analyzed the properties
of the approximated matrices and showed that MEKA frequently generates non-psd
matrices and a correction approach is inevitable. Accordingly, we suggested an
efficient shift correction, based on the Lanczos iteration, to fulfill the psd
assumption of kernel methods, while maintaining the memory efficiency. This
allows us now to approximate indefinite kernel functions with MEKA. Our
evaluation demonstrated the importance of the psd correction and that the
revised MEKA approach is not only more flexible, but remains competitive in
accuracy and memory consumption. In particular, the technique overcomes
limitations of the original MEKA and enables the valid usage of a wide range
of kernel functions in large-scale data analysis as demanded in recent work
(Lee et al.,, 2018; Loosli,, 2019; Liu et al.,, 2020). As the revised MEKA
technique is now applicable to approximate a variety of kernel functions and
leads to consistent approximations, there are other properties of the original
kernel function which are not yet preserved. For example, the approximated
Gaussian rbf kernel may contain negative similarities, although being psd.
This is not meaningful and can be handled by truncations followed by the
suggested shift correction. Here it may be desirable to enforce additional
constraints in the least squares approximation of the off-diagonal blocks.
However, this can not be easily achieved in closed form and will be subject of
further research.

#### Acknowledgements

SH is supported by the ESF (WiT-HuB 4/2014-2020), project KI-trifft-KMU,
StMBW-W-IX.4-6-190065. MM is supported by the Bavarian HighTech agenda and the
Würzburg Center for Artificial Intelligence and Robotics (CAIRO).

## References

  * Coakley and Rokhlin, (2013) Coakley, E. S. and Rokhlin, V. (2013).  A fast divide-and-conquer algorithm for computing the spectra of real symmetric tridiagonal matrices.  Applied and Computational Harmonic Analysis, 34(3):379–414.
  * Cortes and Vapnik, (1995) Cortes, C. and Vapnik, V. (1995).  Support-vector networks.  Machine learning, 20(3):273–297.
  * Cullum and Willoughby, (2002) Cullum, J. and Willoughby, R. (2002).  Lanczos algorithms for large symmetric eigenvalue computations: Vol. I: Theory.  SIAM, Philadelphia (PA).
  * Drineas et al., (2008) Drineas, P., Mahoney, M., and Muthukrishnan, S. (2008).  Relative-error CUR matrix decompositions.  SIAM Journal on Matrix Analysis and Applications, 30(2):844–881.
  * Duin and Pekalska, (2005) Duin, R. and Pekalska, E. (2005).  Dissimilarity representation for pattern recognition, foundations and applications, volume 64.  World Scientific, Singapore.
  * Duin and Pekalska, (2010) Duin, R. P. W. and Pekalska, E. (2010).  Non-Euclidean Dissimilarities: Causes and Informativeness.  In Hancock, E. R., Wilson, R. C., Windeatt, T., Ulusoy, I., and Escolano, F., editors, Structural, Syntactic, and Statistical Pattern Recognition, Joint IAPR International Workshop, SSPR&SPR 2010, Cesme, Izmir, Turkey, August 18-20, 2010. Proceedings, volume 6218 of Lecture Notes in Computer Science, pages 324–333. Springer.
  * Frénay and Verleysen, (2011) Frénay, B. and Verleysen, M. (2011).  Parameter-insensitive kernel in extreme learning for non-linear support vector regression.  Neurocomputing, 74(16):2526–2531.
  * Gershgorin, (1931) Gershgorin, S. (1931).  Über die Abgrenzung der Eigenwerte einer Matrix.  Bull. Acad. Sci. URSS, 1931(6):749–754.
  * Gisbrecht and Schleif, (2015) Gisbrecht, A. and Schleif, F. (2015).  Metric and non-metric proximity transformations at linear costs.  Neurocomputing, 167:643–657.
  * Halko et al., (2011) Halko, N., Martinsson, P.-G., and Tropp, J. (2011).  Finding structure with randomness: probabilistic algorithms for constructing approximate matrix decompositions.  SIAM review, 53(2):217–288.
  * Hofmann et al., (2008) Hofmann, T., Schölkopf, B., and Smola, A. (2008).  Kernel methods in machine learning.  The Annals of Statistics, pages 1171–1220.
  * Huang et al., (2018) Huang, X., Suykens, J. A. K., Wang, S., Hornegger, J., and Maier, A. (2018).  Classification With Truncated $\ell_{1}$ Distance Kernel.  IEEE Transactions on Neural Networks and Learning Systems, 29(5):2025–2030.
  * Kumar and Schneider, (2017) Kumar, K. and Schneider, J. (2017).  Literature survey on low rank approximation of matrices.  Linear and Multilinear Algebra, 65(11):2212–2244.
  * Laub et al., (2006) Laub, J., Roth, V., Buhmann, J. M., and Müller, K. (2006).  On the information and representation of non-Euclidean pairwise data.  Pattern Recognition, 39(10):1815–1826.
  * Le et al., (2013) Le, Q., Sarlós, T., and Smola, A. (2013).  Fastfood-approximating kernel expansions in loglinear time.  In Proceedings of International Conference on Machine Learning, pages 244–252.
  * Lee et al., (2018) Lee, J., Bahri, Y., Novak, R., Schoenholz, S., Pennington, J., and Sohl-Dickstein, J. (2018).  Deep neural networks as Gaussian processes.  In 6th International Conference on Learning Representations, ICLR 2018 - Conference Track Proceedings.
  * Liu et al., (2021) Liu, F., Huang, X., Chen, Y., and Suykens, J. (2021).  Fast Learning in Reproducing Kernel Krein Spaces via Signed Measures.  In International Conference on Artificial Intelligence and Statistics, pages 388–396. PMLR.
  * Liu et al., (2018) Liu, F., Huang, X., Gong, C., Yang, J., and Suykens, J. A. (2018).  Indefinite kernel logistic regression with concave-inexact-convex procedure.  IEEE transactions on neural networks and learning systems, 30(3):765–776.
  * Liu et al., (2020) Liu, H., Ong, Y., Shen, X., and Cai, J. (2020).  When Gaussian Process Meets Big Data: A Review of Scalable GPs.  IEEE Transactions on Neural Networks and Learning Systems, 31(11):4405–4423.
  * Loosli, (2019) Loosli, G. (2019).  TrIK-SVM : an alternative decomposition for kernel methods in Krein spaces.  In 27th European Symposium on Artificial Neural Networks, ESANN 2019, Bruges, Belgium, April 24-26, 2019.
  * Loosli et al., (2015) Loosli, G., Canu, S., and Ong, C. S. (2015).  Learning SVM in kreĭn spaces.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 38(6):1204–1216.
  * Luss and d’Aspremont, (2009) Luss, R. and d’Aspremont, A. (2009).  Support vector machine classification with indefinite kernels.  Mathematical Programming Computation, 1(2):97–118.
  * Madeo et al., (2013) Madeo, R., Lima, C., and Peres, S. (2013).  Gesture unit segmentation using support vector machines: segmenting gestures from rest positions.  In Proceedings of the 28th Annual ACM Symposium on Applied Computing, pages 46–52.
  * Mahadevan, (2008) Mahadevan, S. (2008).  Fast spectral learning using lanczos eigenspace projections.  In AAAI, pages 1472–1475.
  * Mehrkanoon et al., (2018) Mehrkanoon, S., Huang, X., and Suykens, J. A. (2018).  Indefinite kernel spectral learning.  Pattern Recognition, 78:144–153.
  * Münch et al., (2020) Münch, M., Raab, C., Biehl, M., and Schleif, F.-M. (2020).  Data-Driven Supervised Learning for Life Science Data.  Frontiers in Applied Mathematics and Statistics, 6:56.
  * Morrison et al., (2016) Morrison, D. R., Jacobson, S. H., Sauppe, J. J., and Sewell, E. C. (2016).  Branch-and-bound algorithms: A survey of recent advances in searching, branching, and pruning.  Discrete Optimization, 19:79–102.
  * Münch et al., (2020) Münch, M., Raab, C., Biehl, M., and Schleif, F.-M. (2020).  Structure Preserving Encoding of Non-euclidean Similarity Data.  In ICPRAM, pages 43–51.
  * Oglic and Gärtner, (2019) Oglic, D. and Gärtner, T. (2019).  Scalable learning in reproducing kernel krein spaces.  In International Conference on Machine Learning, pages 4912–4921. PMLR.
  * Parlett et al., (1982) Parlett, B. N., Simon, H., and Stringer, L. (1982).  On estimating the largest eigenvalue with the Lanczos algorithm.  Mathematics of computation, 38(157):153–165.
  * Pennington et al., (2015) Pennington, J., Yu, F. X., and Kumar, S. (2015).  Spherical random features for polynomial kernels.  In Advances in Neural Information Processing Systems, pages 1837–1845. MIT Press.
  * Rahimi et al., (2007) Rahimi, A., Recht, B., et al. (2007).  Random features for large-scale kernel machines.  In Proceedings in Neural Information Processing Systems. MIT Press.
  * Rasmussen and Williams, (2006) Rasmussen, C. E. and Williams, C. K. I. (2006).  Gaussian processes for machine learning.  Adaptive computation and machine learning. MIT Press.
  * Schleif and Tiño, (2017) Schleif, F. and Tiño, P. (2017).  Indefinite core vector machine.  Pattern Recognit., 71:187–195.
  * Schleif et al., (2018) Schleif, F.-M., Gisbrecht, A., and Tino, P. (2018).  Supervised low rank indefinite kernel approximation using minimum enclosing balls.  Neurocomputing, 318:213–226.
  * Schleif and Tino, (2015) Schleif, F.-M. and Tino, P. (2015).  Indefinite proximity learning: a review.  Neural Computation, 27(10):2039–2096.
  * Si et al., (2017) Si, S., Hsieh, C.-J., and Dhillon, I. (2017).  Memory efficient kernel approximation.  The Journal of Machine Learning Research, 18(1):682–713.
  * Sun et al., (2015) Sun, S., Zhao, J., and Zhu, J. (2015).  A review of Nyström methods for large-scale machine learning.  Information Fusion, 26:36–48.
  * Williams and Seeger, (2001) Williams, C. and Seeger, M. (2001).  Using the Nyström method to speed up kernel machines.  In Advances in Neural Information Processing Systems, pages 682–688. MIT Press.
  * Wilson and Hancock, (2010) Wilson, R. C. and Hancock, E. R. (2010).  Spherical embedding and classification.  In Joint IAPR International Workshops on Statistical Techniques in Pattern Recognition (SPR) and Structural and Syntactic Pattern Recognition (SSPR), pages 589–599. Springer.