A Large Dimensional Analysis of Least Squares Support Vector Machines

TL;DR

This paper provides a large-dimensional theoretical analysis of LS-SVMs using random matrix theory, showing how their decision functions behave under high-dimensional Gaussian mixture models and applying findings to real datasets.

Contribution

It offers a novel large-dimensional analysis of LS-SVMs, revealing how kernel functions influence their performance in high-dimensional settings.

Findings

Decision function approximates a normal distribution in high dimensions

Explicit dependence of mean and variance on kernel function

Application to MNIST datasets confirms theoretical predictions

Abstract

In this article, a large dimensional performance analysis of kernel least squares support vector machines (LS-SVMs) is provided under the assumption of a two-class Gaussian mixture model for the input data. Building upon recent advances in random matrix theory, we show, when the dimension of data $p$ and their number $n$ are both large, that the LS-SVM decision function can be well approximated by a normally distributed random variable, the mean and variance of which depend explicitly on a local behavior of the kernel function. This theoretical result is then applied to the MNIST and Fashion-MNIST datasets which, despite their non-Gaussianity, exhibit a convincingly close behavior. Most importantly, our analysis provides a deeper understanding of the mechanism into play in SVM-type methods and in particular of the impact on the choice of the kernel function as well as some of their theoretical limits in separating high dimensional Gaussian vectors.

Tables

Table 1. TABLE I: Empirical estimation of differences in means and covariances of MNIST/Fashion-MNIST data (numbers 8 8 8 and 9 9 9 )

	MNIST/Fashion-MNIST without noise	MNIST/Fashion-MNIST with $0\mathrm{dB}$ noise
${\Vert {𝝁}_2-{𝝁}_2\Vert }^2$	$251/483$	$96/197$
$\frac{1}{p}{\left(\mathrm{tr}\left({𝐂}_2-{𝐂}_1\right)\right)}^2$	$19/89$	$3/13$
$\frac{1}{p}\mathrm{tr}\left({\left({𝐂}_2-{𝐂}_1\right)}^2\right)$	$30/86$	$5/13$

Equations

$$\underset{\mathbf{w},b}{\mathrm{argmin}}$$

$$y_i={\mathbf{w}}^{\mathsf{T}}\mathbf{\phi }({\mathbf{x}}_i)+b+e_i,i=1,\dots ,n$$

$$\{\begin{array}{ll}\mathbf{\alpha }& ={\mathbf{S}}^{-1}\left({\mathbf{I}}_n-\frac{1_n{1}_n^{\mathsf{T}}{\mathbf{S}}^{-1}}{1_n^{\mathsf{T}}{\mathbf{S}}^{-1}{1}_n}\right)\mathbf{y}={\mathbf{S}}^{-1}\left(\mathbf{y}-b{1}_n\right)\\ b& =\frac{1_n^{\mathsf{T}}{\mathbf{S}}^{-1}\mathbf{y}}{1_n^{\mathsf{T}}{\mathbf{S}}^{-1}{1}_n}\end{array}$$

$$g(\mathbf{x})={\mathbf{\alpha }}^{\mathsf{T}}\mathbf{k}(\mathbf{x})+b$$

$$\mathbf{K}$$

$$\mathbf{k}(\mathbf{x})$$

$${\mathbf{x}}_i={\mathbf{\mu }}_a+\sqrt{p}{\mathbf{\omega }}_i$$

$${\left(\mathbf{x}-{\mathbf{\mu }}_2\right)}^{\mathsf{T}}{\mathbf{C}}_2^{-1}\left(\mathbf{x}-{\mathbf{\mu }}_2\right)-{\left(\mathbf{x}-{\mathbf{\mu }}_1\right)}^{\mathsf{T}}{\mathbf{C}}_1^{-1}\left(\mathbf{x}-{\mathbf{\mu }}_1\right)\lessgtr \log \frac{\det {\mathbf{C}}_1}{\det {\mathbf{C}}_2}$$

$$t(\mathbf{x})$$

$$-\frac{1}{p}\log \frac{\det {\mathbf{C}}_1}{\det {\mathbf{C}}_2}\lessgtr 0$$

$$t(\mathbf{x})-\hat{t}\stackrel{d}{\to }0.$$

$$\Vert {\mathbf{x}}_i-{\mathbf{x}}_j{\Vert }^2/p\stackrel{a.s.}{\to }\tau$$

$$\mathbf{P}$$

$$\mathbf{\psi }$$

$$\Vert \mathbf{K}-\hat{\mathbf{K}}\Vert \stackrel{a.s.}{\to }0$$

$${\mathbf{\omega }}_{\mathbf{x}}\triangleq (\mathbf{x}-{\mathbf{\mu }}_a)/\sqrt{p}\in {\mathbb{R}}^p,{\psi }_{\mathbf{x}}\triangleq \Vert {\mathbf{\omega }}_{\mathbf{x}}{\Vert }^2-\mathbb{E}\left[\Vert {\mathbf{\omega }}_{\mathbf{x}}{\Vert }^2\right]\in \mathbb{R}.$$

$$\hat{g}(\mathbf{x})=\{\begin{array}{ll}{c}_2-c_1+\gamma \left(\mathfrak{P}-2c_1{c}_2^2\mathfrak{D}\right),& \text{if}\mathbf{x}\in {\mathcal{C}}_1\\ c_2-c_1+\gamma \left(\mathfrak{P}+2c_1^2{c}_2\mathfrak{D}\right),& \text{if}\mathbf{x}\in {\mathcal{C}}_2\end{array}$$

$$\mathfrak{P}$$

$$+2c_1{c}_2{f}^{\prime \prime }(\tau ){\psi }_{\mathbf{x}}\frac{\mathrm{tr}\left({\mathbf{C}}_2-{\mathbf{C}}_1\right)}{p}$$

$$\mathfrak{D}$$

$$+\frac{2f^{\prime \prime }(\tau )}{p^2}\mathrm{tr}\left(({\mathbf{C}}_2-{\mathbf{C}}_1{)}^2\right).$$

$$G_a\sim \mathcal{N}({\mathrm{E}}_a,{\mathrm{Var}}_a)$$

$${\mathrm{E}}_a=\{\begin{array}{ll}{c}_2-c_1-2c_2\cdot c_1{c}_2\gamma \mathfrak{D},& a=1\\ c_2-c_1+2c_1\cdot c_1{c}_2\gamma \mathfrak{D},& a=2\end{array}$$

$${\mathrm{Var}}_a=8{\gamma }^2{c}_1^2{c}_2^2\left({\mathcal{V}}_1^a+{\mathcal{V}}_2^a+{\mathcal{V}}_3^a\right)$$

$${\mathcal{V}}_1^a$$

$${\mathcal{V}}_2^a$$

$${\mathcal{V}}_3^a$$

$$g(\mathbf{x})-(c_2-c_1)=2c_1{c}_2{g}^{\ast }(\mathbf{x}).$$

$$g(x)$$

$$g(\mathbf{x})-(c_2-c_1)$$

$$=2c_1{c}_2{\left(\frac{\mathbf{y}-(c_2-c_1)1_n}{2c_1{c}_2}\right)}^{\mathsf{T}}\mathbf{\varpi }$$

$$=2c_1{c}_2({\mathbf{y}}^{\ast }{)}^{\mathsf{T}}\mathbf{\varpi }=2c_1{c}_2{g}^{\ast }(\mathbf{x})$$

$$G_a^{\ast }\sim \mathcal{N}({\mathrm{E}}_a^{\ast },{\mathrm{Var}}_a^{\ast })$$

$${\mathrm{E}}_{\mathrm{a}}^{\ast }=\{\begin{array}{ll}-c_2\gamma \mathfrak{D},& a=1\\ +c_1\gamma \mathfrak{D},& a=2\end{array}$$

$$\mathrm{Va}{\mathrm{r}}_{\mathrm{a}}^{\ast }=2{\gamma }^2\left({\mathcal{V}}_1^a+{\mathcal{V}}_2^a+{\mathcal{V}}_3^a\right)$$

$$\mathrm{P}(g(\mathbf{x})>{\xi }_n|\mathbf{x}\in {\mathcal{C}}_1)-Q\left(\frac{{\xi }_n-{\mathrm{E}}_1}{\sqrt{{\mathrm{Var}}_1}}\right)\to 0$$

$$\mathrm{P}(g(\mathbf{x})<{\xi }_n|\mathbf{x}\in {\mathcal{C}}_2)-Q\left(\frac{{\mathrm{E}}_2-{\xi }_n}{\sqrt{{\mathrm{Var}}_2}}\right)\to 0$$

Full text

I Introduction

In the past two decades, due to their surprising classification capability and simple implementation, kernel support vector machine (SVM) [1] and its variants [2, 3, 4] have been used in a wide variety of classification applications, such as face detection [5, 6], handwritten digit recognition [7], and text categorization [8, 9]. In all aforementioned applications, the dimension of data $p$ and their number $n$ are large: in the hundreds and even thousands. The significance of working in this large $n,p$ regime is even more convincing in the Big Data paradigm today where handling data which are both numerous and large dimensional becomes increasingly common.

Firmly grounded in the framework of statistical learning theory [10], support vector machine has two main features: (i) in SVM, the training data $\mathbf{x}_{1},\ldots,\mathbf{x}_{n}\in\mathbb{R}^{p}$ are mapped into some feature space through a non-linear function $\bm{\varphi}$, which, thanks to the so-called “kernel trick” [11], needs not be computed explicitly, so that some kernel function $f$ is introduced in place of the inner product in the feature space: $f(\mathbf{x},\mathbf{y})=\bm{\varphi}(\mathbf{x})^{{\sf T}}\bm{\varphi}(\mathbf{y})$, and (ii) a standard (convex) optimization method is used to find the classifier that both minimizes the training error and yields a good generalization performance for unknown data.

As the training of SVMs involves a quadratic programming problem, the computation complexity of SVM training algorithms can be intensive when the number of training examples $n$ becomes large (at least quadratic with respect to $n$). It is thus difficult to deal with large scale problems with traditional SVMs. To cope with this limitation, least squares SVM (LS-SVM, also later referred to as kernel regularized least-squares estimator or kernel ridge regression [12, 13, 14]) was proposed in [2], providing a more computationally efficient implementation of the traditional SVMs, by taking equality optimization constraints instead of inequalities, which results in an explicit solution (from a set of linear equations) rather than an implicit one in SVMs. This article is mostly concerned with this particular type of SVMs.

Trained SVMs are strongly data-dependent: the data with generally unknown statistics are passed through a nonlinear kernel function $f$ and standard optimization methods are used to find the best classifier. All these features make the performance of SVM hardly traceable (at least within the classical finite $n,p$ regime). To understand the mechanism of SVMs, the notion of VC dimension was introduced to provide bounds on the generalization performance of SVM [10], while a probabilistic interpretation of LS-SVM was discussed in [15] through a Bayesian inference approach. In other related works, connections between LS-SVMs and SVMs were revealed in [16], and more relationships were shown between SVM-type and other learning methods, e.g., LS-SVMs and extreme learning machines (ELMs) [17]; SVMs and regularization networks (RNs) [18], etc. Theoretical analyses on the generalization performance of LS-SVM have been developed, under the conventional asymptotic statistics framework (i.e., assuming $n\to\infty$), to obtain optimal convergence rates in [13, 14]. Nonetheless, a proper adaptation to the large $n,p$ setting to address LS-SVM performance for large dimensional datasets (of growing interest today) is still missing.

Similar to classical analysis of asymptotic statistics where $n\to\infty$ while $p$ is fixed, where the diversity of the number of data provides convergence through laws of large numbers, working in the large $n,p$ regime by letting in addition $p\to\infty$ helps exploit the diversity offered by the size of each data vector, providing us with another dimension to guarantee the convergence of some key objects in our analysis, and thus makes the asymptotic analysis of the elusive kernel matrix $\mathbf{K}=\left\{f\left(\mathbf{x}_{i},\mathbf{x}_{j}\right)\right\}_{i,j=1}^{n}$ technically more accessible. Recent breakthroughs in random matrix theory have allowed one to overtake the theoretical difficulty posed by the nonlinearity of the aforementioned kernel function $f$ [19, 20] and thus make an in-depth analysis of LS-SVM possible in the large $n,p$ regime. These tools were notably used to assess the performance of the popular Ng-Weiss-Jordan kernel spectral clustering methods for large datasets [20], in the analysis of graphed-based semi-supervised learning [21] or for the development of novel kernel subspace clustering methods [22].

Similar to these works, in this article, we provide a performance analysis of LS-SVM, in the regime of $n,p\to\infty$ and $p/n\to\bar{c}_{0}\in(0,\infty)$, under the assumption of a two-class Gaussian mixture model of means $\bm{\mu}_{1},\bm{\mu}_{2}$ and covariance matrices $\mathbf{C}_{1},\mathbf{C}_{2}$ for the input data. The Gaussian assumption may seem artificial to the practitioners, but reveals first insights into how SVM-type methods deal with the information in means and covariances from a more quantitative point of view. Besides, the early investigations [20, 21] have revealed that the behavior of some machine learning methods under Gaussian or deterministic practical input datasets are a close match, despite the obvious non-Gaussianity of the latter.

Our main finding is that, as in [20], in the large $n,p$ regime and under suitable conditions on the input statistics, a non-trivial asymptotic classification error rate (i.e., neither 0 nor 1) can be obtained and the decision function of LS-SVM converges to a Gaussian random variable whose mean and variance depend on the statistics of the two different classes as well as on the behavior of the kernel function $f$ evaluated at $2\operatorname{{\rm tr}}(n_{1}\mathbf{C}_{1}+n_{2}\mathbf{C}_{2})/(np)$, with $n_{1}$ and $n_{2}$ the number of instances in each class. This brings novel insights into some key issues of SVM-type methods such as kernel function selection and parameter optimization (see for example [15, 23, 24, 25, 26, 27] and the references therein), as far as large dimensional data are concerned. More importantly, we confirm through simulations that our theoretical findings closely match the performance obtained on the MNIST [28] and the Fashion-MNIST datasets [29], which conveys a strong applicative motivation for this work.

In the remainder of the article, we provide a rigorous statement of our main results. The problem of LS-SVM is discussed in Section II and our model and main results presented in Section III, while all proofs are deferred to the appendices in the Supplementary Material. In Section IV, attention will be paid on some special cases that are more analytically tractable. Section V concludes the paper by summarizing the main results and outlining future research directions.

Reproducibility: Python 3 codes to reproduce the results in this article are available at https://github.com/Zhenyu-LIAO/RMT4LSSVM.

Notations: Boldface lowercase (uppercase) characters stand for vectors (matrices), and scalars non-boldface respectively. $\mathbf{1}_{n}$ is the column vector of ones of size $n$, $\mathbf{0}_{n}$ the column vector of zeros, and $\mathbf{I}_{n}$ the $n\times n$ identity matrix. The notation $(\cdot)^{{\sf T}}$ denotes the transpose operator. The norm $\|\cdot\|$ is the Euclidean norm for vectors and the operator norm for matrices. The notation ${\rm P}(\cdot)$ denotes the probability measure of a random variable. The notation $\overset{d}{\to}$ denotes convergence in distribution and $\overset{\rm{a.s.}}{\to}$ almost sure convergence, respectively. The operator $\mathcal{D}(\mathbf{v})=\mathcal{D}\{v_{a}\}_{a=1}^{k}$ is the diagonal matrix having $v_{a},\ldots,v_{k}$ as its ordered diagonal elements. We denote $\{v_{a}\}_{a=1}^{k}$ a column vector with $a$-th entry (or block entry) $v_{a}$ (which may be a vector), while $\{V_{ab}\}_{a,b=1}^{k}$ denotes a square matrix with entry (or block-entry) $(a,b)$ given by $V_{ab}$ (which may be a matrix).

II Problem statement

Least squares support vector machines (LS-SVMs) are a modification of the standard SVM introduced in [2] to overcome the drawbacks of SVM related to computational efficiency. The optimization problem has half the number of parameters and benefits from solving a linear system of equations instead of a quadratic programming problem as in standard SVM and is thus more practical for large dimensional learning tasks. In this article, we will focus on a binary classification problem using LS-SVM as described in the following paragraph.

Given a training set $\left\{(\mathbf{x}_{1},y_{1}),\ldots,(\mathbf{x}_{n},y_{n})\right\}$ of size $n$, where data $\mathbf{x}_{i}\in\mathbb{R}^{p}$ and labels $y_{i}\in\{-1,1\}$, the objective of LS-SVM is to devise a decision function $g(\mathbf{x})$ that ideally maps all $\mathbf{x}_{i}$ in the training set to $y_{i}$ and subsequently all unknown data $\mathbf{x}$ to their corresponding $y$ value. Here we denote $\mathbf{x}_{i}\in\mathcal{C}_{1}$ if $y_{i}=-1$ and $\mathbf{x}_{i}\in\mathcal{C}_{2}$ if $y_{i}=1$ and shall say that $\mathbf{x}_{i}$ belongs to class $\mathcal{C}_{1}$ or class $\mathcal{C}_{2}$, respectively. Due to the often nonlinear separability of these training data in the input space $\mathbb{R}^{p}$, in most cases, one associates the training data $\mathbf{x}_{i}$ to some feature space $\mathcal{H}$ through a nonlinear mapping $\bm{\varphi}:\ \mathbf{x}_{i}\mapsto\bm{\varphi}(\mathbf{x}_{i})\in\mathcal{H}$. Constrained optimization methods are then used to define a separating hyperplane in $\mathcal{H}$ with direction vector $\mathbf{w}$ and correspondingly to find a function $g(\mathbf{x})=\mathbf{w}^{{\sf T}}\bm{\varphi}(\mathbf{x})+b$ that minimizes the training errors $e_{i}=y_{i}-\left(\mathbf{w}^{{\sf T}}\bm{\varphi}(\mathbf{x}_{i})+b\right)$, and meanwhile yields good generalization performance by minimizing the norm of $\mathbf{w}$ [30]. More specifically, the LS-SVM approach consists in minimizing the squared errors $e_{i}^{2}$, thus resulting in111 We include the bias term $b$ as in the (classical) LS-SVM formulation [2], which may be different from kernel ridge regression in some literature [12, 15] where no bias term is used.

[TABLE]

where $\gamma>0$ is a penalty factor that weights the structural risk $\|\mathbf{w}\|^{2}$ against the empirical one $\frac{1}{n}\sum_{i=1}^{n}e_{i}^{2}$.

The problem can be solved by introducing Lagrange multipliers $\alpha_{i},i=1,\ldots,n$ with solution $\mathbf{{w}}=\sum_{i=1}^{n}\alpha_{i}\bm{\varphi}(\mathbf{x}_{i})$, where, letting $\mathbf{{y}}=[y_{1},\ldots,y_{n}]^{{\sf T}},\mathbf{\mathbb{\bm{\alpha}}}=[\alpha_{1},\ldots,\alpha_{n}]^{{\sf T}}$, we obtain

[TABLE]

with $\mathbf{S}={\mathbf{K}}+\frac{n}{\gamma}{\bf I}_{n}$ and $\mathbf{K}\triangleq\left\{\bm{\varphi}(\mathbf{x}_{i})^{{\sf T}}\bm{\varphi}(\mathbf{x}_{j})\right\}_{i,j=1}^{n}$ referred to as the kernel matrix [2].

Given $\bm{\alpha}$ and $b$, a new datum $\mathbf{x}$ is then classified into class $\mathcal{C}_{1}$ or $\mathcal{C}_{2}$ depending on the value of the following decision function

[TABLE]

where $\mathbf{k}(\mathbf{x})=\left\{\bm{\varphi}(\mathbf{x})^{{\sf T}}\bm{\varphi}(\mathbf{x}_{j})\right\}_{j=1}^{n}\in\mathbb{R}^{n}$. More precisely, ${\mathbf{x}}$ is associated to class $\mathcal{C}_{1}$ if $g({\bf x})$ takes a small value (below a certain threshold $\xi$) and to class $\mathcal{C}_{2}$ otherwise.222Since data from $\mathcal{C}_{1}$ are labeled $-1$ while data from $\mathcal{C}_{2}$ are labeled $1$.

With the “kernel trick” [11], as shown in (2) and (3) that, both in the “training” and “testing” steps, one only needs to evaluate the inner product $\bm{\varphi}(\mathbf{x}_{i})^{{\sf T}}\bm{\varphi}(\mathbf{x}_{j})$ or $\bm{\varphi}(\mathbf{x})^{{\sf T}}\bm{\varphi}(\mathbf{x}_{j})$, and never needs to know explicitly the mapping $\bm{\varphi}(\cdot)$. In the rest of this article, we assume that the kernel is translation invariant and focus on kernel functions $f:~{}\mathbb{R}^{+}\to\mathbb{R}^{+}$ that satisfy $\bm{\varphi}(\mathbf{x}_{i})^{{\sf T}}\bm{\varphi}(\mathbf{x}_{j})=f(\|\mathbf{x}_{i}-\mathbf{x}_{j}\|^{2}/p)$ and shall redefine $\mathbf{K}$ and $\mathbf{k}(\mathbf{x})$ for data point $\mathbf{x}$ as333As shall be seen later, the division by $p$ here is a convenient normalization in the large $n,p$ regime. For example, we have the (normalized) norm $\|\mathbf{x}_{i}\|/\sqrt{p}$ is of order $O(1)$ with high probability for large $n,p$. The motivation of studying “translation invariant” kernel is that, being one of the most popular types of kernel used in practice, it offers (additionally) technical tractability as a result of the “concentration” phenomenon of large dimensional Gaussian vector, as we shall see later for example in (5). Similar results can be obtained for “inner-product” kernel of the type $f(\mathbf{x}_{i}^{\sf T}\mathbf{x}_{j}/p)$ as presented in [19, 31].

[TABLE]

Some commonly used kernel functions are the Gaussian radial basis (RGB) kernel $f(x)=\exp\left(-\frac{x}{2\sigma^{2}}\right)$ with $\sigma>0$ and the polynomial kernel $f(x)=\sum_{i=0}^{d}{a_{i}x^{i}}$ with $d\geq 1$.

In the rest of this article, we will focus on the performance of LS-SVM, in the large $n,p$ regime, by studying the asymptotic behavior of the decision function $g(\mathbf{x})$ defined in (3), in a binary classification problem with some statistical properties of the data, the model of which will be specified in the next section.

III Main results

III-A Model and assumptions

Evaluating the performance of LS-SVM is made difficult by the heavily data-driven aspect of the method. In this article, we assume that all $\mathbf{x}_{i}$’s are extracted from a Gaussian mixture, thereby allowing for a thorough theoretical analysis.

Let $\mathbf{x}_{1},\ldots,\mathbf{x}_{n}\in\mathbb{R}^{p}$ be independent vectors belonging to two distribution classes $\mathcal{C}_{1},\mathcal{C}_{2}$, with $\mathbf{x}_{1},\ldots,\mathbf{x}_{n_{1}}\in\mathcal{C}_{1}$ and $\mathbf{x}_{n_{1}+1},\dots,\mathbf{x}_{n}\in\mathcal{C}_{2}$ (so that class $\mathcal{C}_{1}$ has cardinality $n_{1}$ and class $\mathcal{C}_{2}$ has cardinality $n-n_{1}=n_{2}$). We assume that $\mathbf{x}_{i}\in\mathcal{C}_{a}$ for $a\in\{1,2\}$ if

[TABLE]

for some $\bm{\mu}_{a}\in\mathbb{R}^{p}$ and $\bm{\omega}_{i}\sim\mathcal{N}(0,\mathbf{C}_{a}/p)$, with $\mathbf{C}_{a}\in\mathbb{R}^{p\times p}$ some positive definite matrix.

As the ${\bm{\mu}}_{a}$’s and ${\mathbf{C}}_{a}$’s scale with $p$, to avoid asymptotic trivial misclassification rates (i.e., neither [math] or $1$ in the limit of $n,p\to\infty$), we shall (as in [32, 20]) technically place ourselves under the following controlled growth rate assumption:

Assumption 1** (Growth Rate).**

As $n\to\infty$, for $a\in\{1,2\}$, the following conditions hold.

• •

Data scaling*: $\frac{p}{n}\triangleq c_{0}\to\bar{c}_{0}>0$.*

• •

Class scaling*: $\frac{n_{a}}{n}\triangleq c_{a}\to\bar{c}_{a}>0$.*

• •

Mean scaling*: $\|\bm{\mu}_{2}-\bm{\mu}_{1}\|=O(1)$.*

• •

Covariance scaling*: $\|\mathbf{C}_{a}\|=O(1)$ and $\operatorname{{\rm tr}}(\mathbf{C}_{2}-\mathbf{C}_{1})=O(\sqrt{p})$.*

• •

for $\mathbf{C}^{\circ}\triangleq\frac{n_{1}}{n}\mathbf{C}_{1}+\frac{n_{2}}{n}\mathbf{C}_{2}$, $\frac{2}{p}\operatorname{{\rm tr}}\mathbf{C}^{\circ}\to\tau>0$ as $n,p\to\infty$.

From a practical aspect, where $p$ and $n$ are fixed quantities, the dual condition $n\to\infty$ and $\frac{p}{n}\to\bar{c}_{0}>0$ must be understood as requesting that both $p$ and $n$ be large and such that the ratio $\frac{p}{n}$ is sufficiently distinct from [math] and $\infty$.444As a matter of fact, as our results will demonstrate, the case where $\frac{p}{n}\to\bar{c}_{0}=0$ is also valid as an extension by continuity through $\bar{c}_{0}\to 0$.

Aside from the last assumption, stated here mostly for technical convenience, it can be shown that the growth rate demanded in Assumption 1 is rate-optimal in the sense that an oracle Neyman–Pearson hypothesis testing procedure (with known $\bm{\mu}_{a}$ and $\mathbf{C}_{a}$) is (in general) ineffective at any smaller distance rates (so that the misclassification rate will constantly be $1$), as discussed in the following remark.

Remark 1** (Optimal Growth Rate).**

Assume that both $\|\mathbf{C}_{a}\|$ and $\|\mathbf{C}_{a}^{-1}\|$ are of order $O(1)$ and let $\mathbf{x}$ be a vector belonging to class $\mathcal{C}_{1}$, i.e., $\mathbf{x}\sim\mathcal{N}(\bm{\mu}_{1},\mathbf{C}_{1})$. Then, for perfectly known means $\bm{\mu}_{1},\bm{\mu}_{2}$ and covariances $\mathbf{C}_{1},\mathbf{C}_{2}$, the Neyman–Pearson test for $\mathbf{x}$ to belong to $\mathcal{C}_{1}$ consists in the following comparison,

[TABLE]

which is further equivalent to

[TABLE]

where we denote $\Delta\bm{\mu}\triangleq\bm{\mu}_{1}-\bm{\mu}_{2}$, $\bm{\omega}\triangleq\frac{1}{\sqrt{p}}(\mathbf{x}-\bm{\mu}_{1})$ and thus $\bm{\omega}\sim\mathcal{N}(0,\mathbf{C}_{1}/p)$. To explore the difference in means $\Delta\bm{\mu}$ we take $\mathbf{C}_{1}=\mathbf{C}_{2}=\mathbf{C}$ and by Lyapunov’s CLT [33, Theorem 27.3] we have, as $p\to\infty$,

[TABLE]

where $\hat{t}\sim\mathcal{N}\left(\frac{1}{p}\Delta\bm{\mu}^{{\sf T}}\mathbf{C}^{-1}\Delta\bm{\mu},\frac{2}{p}\Delta\bm{\mu}^{{\sf T}}\mathbf{C}^{-1}\Delta\bm{\mu}\right)$.

For a non-trivial classification rate, the mean of $\hat{t}$ must scale with $p$ at least at the same rate as its standard deviation and thus, since $\|\mathbf{C}_{a}^{-1}\|=O(1)$, this implies that $\|\bm{\Delta\mu}\|$ be at least of order $O(1)$. Similar analysis can be performed to obtain the rate $\|\mathbf{C}_{1}-\mathbf{C}_{2}\|=O(1/\sqrt{p})$ and consequently $\operatorname{{\rm tr}}(\mathbf{C}_{2}-\mathbf{C}_{1})=O(\sqrt{p})$. We refer the readers to [32] for more discussions in this respect.

A key observation, also made in [20], is that, as a consequence of Assumption 1, for all pairs $i\neq j$,

[TABLE]

and the convergence is even uniform across all $i\neq j$. This remark is the crux of all subsequent results (note that, surprisingly at first, it states that all data are essentially at the same distance from one another, irrespective of classes, and that the matrix $\mathbf{K}$ defined in (4) has all its entries essentially equal “in the limit” due to the the high dimensional nature of the data; this can be seen as a manifestation of the “curse of dimensionality” with respect to the Euclidean distance in high-dimensional space).

The function $f$ defining the kernel matrix $\mathbf{K}$ in (4) shall be requested to satisfy the following assumption:

Assumption 2** (Kernel Function).**

The function $f$ is a three-times differentiable function in a neighborhood of $\tau$.

The objective of this article is to assess the performance of LS-SVM, under the setting of Assumptions 1 and 2, by studying the asymptotic behavior of the decision function $g(\mathbf{x})$ defined in (3). Following the work of [19] and [20], under our basic settings, the convergence in (5) makes it possible to linearize the kernel matrix $\mathbf{K}$ around the matrix $f(\tau)\mathbf{1}_{n}\mathbf{1}_{n}^{{\sf T}}$, and thus the intractable nonlinear kernel matrix $\mathbf{K}$ can be asymptotically linearized in the large $n,p$ regime. As such, since the decision function $g(\mathbf{x})$ is explicitly defined as a function of $\mathbf{K}$ (through $\bm{\alpha}$ and $b$ as defined in (2)), one can work out an asymptotic linearization of $g(\mathbf{x})$ as a function of the kernel function $f$ and the statistics of the data. This analysis, presented in detail in Appendix A of the Supplementary Material, allows one to reveal the relationship between the performance of LS-SVM and the kernel function $f$ as well as the given learning task, for Gaussian input data as $n,p\to\infty$, as presented in the following subsection.

III-B Asymptotic behavior of the decision function $g(\mathbf{x})$

Before going into our main results, a few notations need to be introduced. In the remainder of the article, we shall use the following deterministic and random elements notations:

[TABLE]

Under Assumptions 1 and 2, following up [20], one can approximate the kernel matrix $\mathbf{K}$ by $\hat{\mathbf{K}}$ in such a way that

[TABLE]

with $\hat{\mathbf{K}}=-2f^{\prime}(\tau)(\mathbf{M}+\mathbf{V}\mathbf{V}^{{\sf T}})+\left(f(0)-f(\tau)+\tau f^{\prime}(\tau)\right)\mathbf{I}_{n}$ for some matrices $\mathbf{M}$ and $\mathbf{V}$, where $\mathbf{M}$ is a standard random matrix model (of operator norm $O(1)$) and $\mathbf{V}\mathbf{V}^{{\sf T}}$ a small rank matrix (of operator norm $O(n)$), which depends both on $\mathbf{P},\bm{\Omega},\bm{\psi}$ and on the class statistics $\bm{\mu}_{1},\bm{\mu}_{2}$ and $\mathbf{C}_{1},\mathbf{C}_{2}$. The same analysis is applied to the vector $\mathbf{k}(\mathbf{x})$ by similarly defining the following random variables for a new datum $\mathbf{x}\in\mathcal{C}_{a}$, $a\in\{1,2\}$:

[TABLE]

Based on the (operator norm) approximation $\mathbf{K}\approx\hat{\mathbf{K}}$, a Taylor expansion is then performed on $\mathbf{S}^{-1}=\left({\mathbf{K}}+n\mathbf{I}_{n}/\gamma\right)^{-1}$ to obtain an (asymptotic) approximation of $\mathbf{S}^{-1}$, and subsequently on $\bm{\alpha}$ and $b$ which depend explicitly on $\mathbf{S}^{-1}$. At last, plugging these results into (3), one finds the main technical result of this article as follows.

Theorem 1** (Asymptotic Approximation).**

Let Assumptions 1 and 2 hold, and $g(\mathbf{x})$ be defined by (3). Then, as $n,p\to\infty$, $n(g({\bf x})-\hat{g}(\mathbf{x}))\overset{\rm{a.s.}}{\to}0$, where

[TABLE]

with

[TABLE]

Leaving the proof to Appendix A in the Supplementary Material, Theorem 1 tells us that the decision function $g(\mathbf{x})$ has an asymptotic equivalent $\hat{g}(\mathbf{x})$ that consists of three parts:

the deterministic term $c_{2}-c_{1}$ of order $O(1)$ that depends on the number of instances in each class of the training set, which essentially comes from the term $\mathbf{1}_{n}^{{\sf T}}\mathbf{y}/n$ in $b$; 2. 2.

the “noisy” term $\mathfrak{P}$ of order $O(n^{-1})$ which is a function of the zero mean random variables $\bm{\omega}_{\mathbf{x}}$ and $\psi_{\mathbf{x}}$, thus in particular $\mathbb{E}[\mathfrak{P}]=0$; 3. 3.

the “informative” term containing $\mathfrak{D}$, also of order $O(n^{-1})$, which features the deterministic differences between the two classes.

From Theorem 1, under the basic settings of Assumption 1, for Gaussian data ${\bf x}\in\mathcal{C}_{a}$, $a\in\{1,2\}$, we can show that $\hat{g}(\mathbf{x})$ (and therefore $g(\mathbf{x})$) converges to a random Gaussian variable the mean and variance of which are given in the following theorem. The proof is deferred to Appendix B.

Theorem 2** (Gaussian Approximation).**

Under the setting of Theorem 1, $n(g({\bf x})-G_{a})\overset{d}{\to}0$, where

[TABLE]

with

[TABLE]

and

[TABLE]

Theorem 2 is our main practical result as it allows one to evaluate the large $n,p$ performance of LS-SVM for Gaussian data. While dwelling on the implications of Theorem 1 and 2, several remarks and discussions are in order.

Remark 2** (Dominant Bias).**

From Theorem 1, under the key Assumption 1, both the random noise $\mathfrak{P}$ and the deterministic “informative” term $\mathfrak{D}$ are of order $O(n^{-1})$, which means that the decision function $g(\mathbf{x})=c_{2}-c_{1}+O(n^{-1})$. This result somehow contradicts the classical decision criterion proposed in [2], based on the sign of $g(\mathbf{x})$, i.e., $\mathbf{x}$ is associated to class $\mathcal{C}_{1}$ if $g({\mathbf{x}})<0$ and to class $\mathcal{C}_{2}$ otherwise. When $c_{1}\neq c_{2}$, this would lead to an asymptotic classification of all new data $\mathbf{x}$’s in the same class as $n\to\infty$. Practically speaking, this means for $n,p$ large that the decision function $g(\mathbf{x})$ of a new datum $\mathbf{x}$ lies (sufficiently) away from [math] ([math] being the classically considered threshold), so that the sign of $g(\mathbf{x})$ is constantly positive (in the case of $\bar{c}_{2}>\bar{c}_{1}$) or negative (in the case of $\bar{c}_{2}<\bar{c}_{1}$). As such, all new data will be trivially classified into the same class. Instead, a first result of Theorem 1 is that the decision threshold $\xi$ should be taken as $\xi=\xi_{n}=c_{2}-c_{1}+O(n^{-1})$ for imbalanced classification problem.

The conclusion of Remark 2 was in fact already known since the work of [15] who reached the same conclusion through a Bayesian inference analysis, for all finite $n,p$. From their Bayesian perspective, the term $c_{2}-c_{1}$ appears in the “bias term” $b$ under the form of prior class probabilities ${\rm P}(y=-1)$, ${\rm P}(y=1)$ and allows for adjusting classification problems with different prior class probabilities in the training and test sets. This idea of a (static) bias term correction has also been applied in [34] in order to improve the validation set performance. Here we confirm the problem of imbalanced datasets in Remark 2 by Figure 1 with $c_{1}=1/4$ and $c_{2}=3/4$, where the histograms of $g(\mathbf{x})$ for $\mathbf{x}\in\mathcal{C}_{1}$ and $\mathcal{C}_{2}$ center somewhere close to $c_{2}-c_{1}=0.5$, thus resulting in a trivial classification by assigning all new data to $\mathcal{C}_{2}$ if one takes $\xi=0$ because ${\rm P}(g(\mathbf{x})<\xi\mid\mathbf{x}\in\mathcal{C}_{1})\to 0$ and ${\rm P}(g(\mathbf{x})>\xi\mid\mathbf{x}\in\mathcal{C}_{2})\to 1$ as $n,p\to\infty$ (the convergence being in fact an equality for finite $n,p$ in this particular figure).

An alternative to alleviate this imbalance issue is to normalize the label vector $\mathbf{y}$. From the proof of Theorem 1 in Appendix A we see the term $c_{2}-c_{1}$ is due to the fact that in $b$ one has $\mathbf{1}_{n}^{{\sf T}}\mathbf{y}/n=c_{2}-c_{1}\neq 0$. Thus, one may normalize the labels $y_{i}$ as $y^{*}_{i}=-1/c_{1}$ if $\mathbf{x}_{i}\in\mathcal{C}_{1}$ and $y^{*}_{i}=1/c_{2}$ if $\mathbf{x}_{i}\in\mathcal{C}_{2}$, so that the relation $\mathbf{1}_{n}^{{\sf T}}\mathbf{y}^{*}=0$ is satisfied. This formulation is also referred to as the Fisher’s targets: $\{-n/n_{1},n/n_{2}\}$ in the context of kernel fisher discriminant analysis [35, 36]. With the aforementioned normalized labels $\mathbf{y}^{*}$, we have the following lemma that reveals the connection between the corresponding decision function $g^{*}(\mathbf{x})$ and $g(\mathbf{x})$.

Lemma 1**.**

Let $g(\mathbf{x})$ be defined by (3) and $g^{*}(\mathbf{x})$ be defined as $g^{*}(\mathbf{x})=(\bm{\alpha}^{*})^{{\sf T}}\mathbf{k}(\mathbf{x})+b^{*}$, with $(\bm{\alpha}^{*},b^{*})$ given by (2) for $\mathbf{y}^{*}$ in the place of $\mathbf{y}$, where $y^{*}_{i}=-1/c_{1}$ if $\mathbf{x}_{i}\in\mathcal{C}_{1}$ and $y^{*}_{i}=1/c_{2}$ if $\mathbf{x}_{i}\in\mathcal{C}_{2}$. Then,

[TABLE]

Proof:

From (2) and (3) we get

[TABLE]

with $\bm{\varpi}=\left(\mathbf{S}^{-1}-\frac{\mathbf{S}^{-1}\mathbf{1}_{n}\mathbf{1}_{n}^{{\sf T}}\mathbf{S}^{-1}}{\mathbf{1}_{n}^{{\sf T}}\mathbf{S}^{-1}\mathbf{1}_{n}}\right)\mathbf{k}(\mathbf{x})+\frac{\mathbf{S}^{-1}\mathbf{1}_{n}}{\mathbf{1}_{n}^{{\sf T}}\mathbf{S}^{-1}\mathbf{1}_{n}}$. Besides, note that $\mathbf{1}_{n}^{{\sf T}}\bm{\varpi}=1$. We thus have

[TABLE]

which concludes the proof. ∎

As a consequence of Lemma 1, instead of Theorem 2 for standard labels $\mathbf{y}$, one would have the following corollary for the corresponding Gaussian approximation of $g^{*}(\mathbf{x})$ when normalized labels $\mathbf{y}^{*}$ are used.

Corollary 1** (Gaussian Approximation of $g^{*}(\mathbf{x})$).**

Under the setting of Theorem 1, and with $g^{*}(\mathbf{x})$ defined in Lemma 1, $n(g^{*}({\bf x})-G^{*}_{a})\overset{d}{\to}0$, where

[TABLE]

with

[TABLE]

and $\mathfrak{D}$ is defined by (8), $\mathcal{V}_{1}^{a},\mathcal{V}_{2}^{a}$ and $\mathcal{V}_{3}^{a}$ as in Theorem 2.

Figure 2 illustrates this result in the same settings as Figure 1. Compared to Figure 1, one can observe that in Figure 2 both histograms are now centered close to [math] (at distance $O(n^{-1})$ from zero) instead of $c_{2}-c_{1}=1/2$. Still, even in the case where normalized labels $\mathbf{y}^{*}$ are used as observed in Figure 2 (where the histograms cross at about $-0.004\approx 1/n$), taking $\xi=0$ as a decision threshold may not be an appropriate choice, as ${\rm E}^{*}_{1}\neq-{\rm E}^{*}_{2}$.

Remark 3** (Insignificance of $\gamma$).**

As a direct result of Theorem 1 and Remark 2, note in (6) that $\hat{g}(\mathbf{x})-\left(c_{2}-c_{1}\right)$ is proportional to the hyperparameter $\gamma$, which indicates that, rather surprisingly, the tuning of $\gamma$ is (asymptotically) of no importance when $n,p\to\infty$ since it does not alter the classification statistics when one uses the sign of $g(\mathbf{x})-\left(c_{2}-c_{1}\right)$ for the decision.555This remark is only valid only under Assumption 1 and $\gamma=O(1)$, i.e., $\gamma$ is considered to remain a constant as $n,p\to\infty$. Recall that this is in sharp contrast with [13] where $\gamma=O(\sqrt{n})$ (or $O(n)$, depending on the problem) is claimed optimal in the large $n$ only regime. From Remark 1 on the growth rate optimality reached by LS-SVM, we see here that $\gamma=O(1)$ is rate-optimal under the present large $n,p$ setting; yet we believe that more elaborate kernels (such as those explored in [37]) may allow for improved performances (not in the rate but in the constants), possibly for different scales of $\gamma$. This intuition will be explored in future investigations.

Letting $Q(x)=\frac{1}{2\pi}\int_{x}^{\infty}\exp\left(-t^{2}/2\right)dt$, from Theorem 2 and Corollary 1, we now have the following immediate corollary for the (asymptotic) classification error rate.

Corollary 2** (Asymptotic Error Rate).**

Under the setting of Theorem 1, for a threshold $\xi_{n}$ possibly depending on $n$, as $n\to\infty$,

[TABLE]

with ${\rm E}_{a}$ and ${\rm Var}_{a}$ given in Theorem 2.

Obviously, Corollary 2 is only meaningful when $\xi_{n}=c_{2}-c_{1}+O(n^{-1})$ as recalled earlier. Besides, it is clear from Lemma 1 and Corollary 1 that ${\rm P}(g(\mathbf{x})>\xi_{n}~{}|~{}\mathbf{x}\in\mathcal{C}_{a})={\rm P}(g^{*}(\mathbf{x})>\xi_{n}-(c_{2}-c_{1})~{}|~{}\mathbf{x}\in\mathcal{C}_{a})$, so that Corollary 2 extends naturally to $g^{*}(\mathbf{x})$ when normalized labels $\mathbf{y}^{*}$ are applied.

Corollary 2 allows one to compute the asymptotic misclassification rate as a function of $\rm{E}_{a},\rm{Var}_{a}$ and the threshold $\xi_{n}$. Combined with Theorem 2, one may note the significance of a proper choice of the kernel function $f$. For instance, if $f^{\prime}(\tau)=0$, the term $\bm{\mu}_{2}-\bm{\mu}_{1}$ vanishes from the mean and variance of $G_{a}$, meaning that the classification of LS-SVM will not rely (at least asymptotically and under Assumption 1) on the differences in means of the two classes. Figure 3 corroborates this finding with the same theoretical Gaussian approximations $G_{1}$ and $G_{2}$ in subfigures (a) and (b). When $\|\bm{\mu}_{2}-\bm{\mu}_{1}\|^{2}$ varies from [math] in (a) to $18$ in (b), the distribution of $g(\mathbf{x})$, and in particular, the overlap between two classes, remain almost the same in (a) and (b).

More traceable special cases and discussions on the choice of kernel function $f$ will be given in the next section.

IV Special cases and further discussions

IV-A More discussions on the kernel function $f$

Following the discussion at the end of Section III, if $f^{\prime}(\tau)=0$, the information about the statistical means of the two different classes is lost and will not help perform the classification. Nonetheless, we find that, rather surprisingly, if one further assumes $\operatorname{{\rm tr}}\mathbf{C}_{1}=\operatorname{{\rm tr}}\mathbf{C}_{2}+o(\sqrt{p})$ (which is beyond the minimum “distance” rate in Assumption 1), using a kernel $f$ that satisfies $f^{\prime}(\tau)=0$ results in ${\rm Var}_{a}=0$ while ${\rm E}_{a}$ may remain non-zero, thereby ensuring a vanishing misclassification rate (as long as $f^{\prime\prime}(\tau)\neq 0$). Intuitively speaking, the kernels with $f^{\prime}(\tau)=0$ play an important role in extracting the information of “shape” of both classes, making the classification extremely accurate even in cases that are deemed impossible to classify according to Remark 1. This phenomenon was also remarked in [20] and deeply investigated in [22]. Figure 4 substantiates this finding for $\bm{\mu}_{1}=\bm{\mu}_{2}$, ${\bf C}_{1}={\bf I}_{p}$ and $\{{\bf C}_{2}\}_{i,j}=.4^{\mid i-j\mid}$, for which $\operatorname{{\rm tr}}\mathbf{C}_{1}=\operatorname{{\rm tr}}\mathbf{C}_{2}=p$. We observe a rapid drop of the classification error as $f^{\prime}(\tau)$ gets close to [math].

Remark 4** (Condition on Kernel Function $f$).**

From Theorem 2 and Corollary 1, one observes that $|{\rm E}_{1}-{\rm E}_{2}|$ is always proportional to the “informative” term $\mathfrak{D}$ and should, for fixed ${\rm Var}_{a}$, be made as large as possible to avoid the overlap of $g(\mathbf{x})$ for $\mathbf{x}$ from different classes. Since ${\rm Var}_{a}$ does not depend on the signs of $f^{\prime}(\tau)$ and $f^{\prime\prime}(\tau)$, it is easily deduced that, to achieve optimal classification performance, one needs to choose the kernel function $f$ such that $f(\tau)>0,f^{\prime}(\tau)<0$ and $f^{\prime\prime}(\tau)>0$.

Incidentally, the condition in Remark 4 is naturally satisfied for Gaussian kernel $f(x)=\exp\left(-x/(2\sigma^{2})\right)$ for any $\sigma$, meaning that, even without specific tuning of the kernel parameter $\sigma$ through cross validation or other techniques, LS-SVM is expected to perform rather well with a Gaussian kernel (as shown in Figure 5), which is not always the case for polynomial kernels. This especially entails, for a second-order polynomial kernel given by $f(x)=a_{2}x^{2}+a_{1}x+a_{0}$, that attention should be paid to meeting the aforementioned condition when tuning the kernel parameters $a_{2},a_{1}$ and $a_{0}$. Figure 6 attests of this remark with Gaussian input data. A rapid increase in classification error rate can be observed both in theory and in practice as soon as the condition $f^{\prime}(\tau)<0,f^{\prime\prime}(\tau)>0$ is no longer satisfied.

Note also from both Figure 4 and Figure 5 that, when $n,p$ are doubled (from $2048,512$ to $4\,096,1\,024$ in Figure 4 and from $256,512$ to $512,1\,024$ in Figure 5), the empirical error becomes closer to the theoretical one, which confirms the asymptotic result as $n,p\to\infty$.

Clearly, for practical use, one needs to know in advance the value of $\tau$ before training so that the kernel $f$ can be properly chosen during the training step. The estimation of $\tau$ is possible, in the large $n,p$ regime, with the following lemma.

Lemma 2**.**

Under Assumptions 1 and 2, as $n\to\infty$,

[TABLE]

with $\bar{\mathbf{x}}\triangleq\frac{1}{n}\sum_{i=1}^{n}{\mathbf{x}_{i}}$.

Proof:

Since

[TABLE]

with $\kappa=\frac{4}{n\sqrt{p}}(\bm{\mu}_{2}-\bm{\mu}_{1})^{{\sf T}}\left(-c_{2}\sum_{\mathbf{x}_{i}\in\mathcal{C}_{1}}\bm{\omega}_{i}+c_{1}\sum_{\mathbf{x}_{j}\in\mathcal{C}_{2}}\bm{\omega}_{j}\right)$ and $\bar{\bm{\omega}}=\frac{1}{n}\sum_{i=1}^{n}\bm{\omega}_{i}$.

According to Assumption 1 we have $\frac{2c_{1}c_{2}}{p}\|\bm{\mu}_{2}-\bm{\mu}_{1}\|^{2}=O(n^{-1})$. The term $\kappa$ is a linear combination of independent zero-mean Gaussian variables and thus $\kappa\sim\mathcal{N}(0,{\rm Var}[\kappa])$ with ${\rm Var}[\kappa]=\frac{16c_{1}c_{2}}{np^{2}}(\bm{\mu}_{2}-\bm{\mu}_{1})^{{\sf T}}\left(c_{2}\mathbf{C}_{1}+c_{1}\mathbf{C}_{2}\right)(\bm{\mu}_{2}-\bm{\mu}_{1})=O(n^{-3})$. We thus deduce from Chebyshev’s inequality and Borel-Cantelli lemma that $\kappa\overset{\rm{a.s.}}{\to}0$.

We then work on the last term $\frac{2}{n}\sum_{i=1}^{n}\|\bm{\omega}_{i}-\bar{\bm{\omega}}\|^{2}$ as

[TABLE]

Since $\bar{\bm{\omega}}\sim\mathcal{N}(\mathbf{0},\mathbf{C}^{\circ}/np)$, we deduce that $\|\bar{\bm{\omega}}\|^{2}\overset{\rm{a.s.}}{\to}0$. Ultimately by the strong law of large numbers, the term $\frac{2}{n}\sum_{i=1}^{n}\|\bm{\omega}_{i}\|^{2}\overset{\rm{a.s.}}{\to}\tau$, which concludes the proof. ∎

IV-B Some limiting cases

IV-B1 Dominant information in means

When $\|\bm{\mu}_{2}-\bm{\mu}_{1}\|^{2}$ is largely dominant over $(\operatorname{{\rm tr}}({\bf C}_{2}-{\bf C}_{1}))^{2}/p$ and $\operatorname{{\rm tr}}((\mathbf{C}_{2}-\mathbf{C}_{1})^{2})/p$, from Theorem 2, both $\rm{E}_{a}-(c_{2}-c_{1})$ and $\sqrt{\rm{Var}_{a}}$ are (approximately) proportional to $f^{\prime}(\tau)$, which eventually makes the choice of the kernel irrelevant (as long as $f^{\prime}(\tau)\neq 0$). This result also holds true for $\rm{E}_{a}^{*}$ and $\sqrt{\rm{Var}_{a}^{*}}$ when normalized labels $\mathbf{y}^{*}$ are applied, as a result of Lemma 1.

IV-B2 $c_{0}$ large or small

Note that, different from both $\mathcal{V}_{1}$ and $\mathcal{V}_{2}$, $\mathcal{V}_{3}$ is a function of $c_{0}$ as it can be rewritten as

[TABLE]

which indicates that the variance of $g(\mathbf{x})$ grows as $c_{0}$ becomes large. This result is easily understood since, with $p$ fixed, a small $c_{0}$ means a larger $n$, and with more training samples, one may “gain” more information of the two different classes, which reduces the “uncertainty” of the classifier. When $n\to\infty$ with a fixed $p$, we have $c_{0}\to 0$ and the LS-SVM is considered “well-trained” and its performance can be described with Theorem 2 by taking $\mathcal{V}_{3}=0$. However, it is worthy noting that the misclassification rate may not be [math] even in this case, since $\mathcal{V}_{1}$ and $\mathcal{V}_{2}$ may differ from [math], which indicated the theoretical limitation of LS-SVM in separating high dimensional Gaussian vectors. On the contrary, when $c_{0}\to\infty$, with few training data, LS-SVM does not sample sufficiently the high dimensional space of the $\mathbf{x}$’s, thus resulting in a classifier with arbitrarily large variance (for fixed means). Moreover, since the term $\mathcal{V}_{3}^{a}$ is proportional to $n^{-1}$, we see that for $f^{\prime}(\tau)$ away from zero and fixed large $p$, as $n$ grows large, the two Gaussians $G_{1}$ and $G_{2}$ in Theorem 2 separate from each other at a rate of $n^{-\frac{1}{2}}$, the overlapping section of the Gaussian tails then provides the misclassification rate via Corollary 2. Figure 7 confirms this result with $p$ fixed to $256$ while $n$ varies from $8$ to $8\,192$.

IV-B3 $c_{1}\to 0$

As revealed in Remark 2, the ratio $c_{1}/c_{2}$ plays a significant role in the performance of classification. A natural question arises: what happens when one class is strongly dominant over the other? Take the case of $c_{1}\to 0,c_{2}\to 1$. From Corollary 1, one has $\rm{E}_{1}^{*}\to-\gamma\mathfrak{D}$, $\rm{E}_{2}^{*}\to 0$ and $\mathcal{V}_{3}^{a}\to\infty$ because of $c_{1}\to 0$ in the denominator, which then makes the ratio $\frac{\rm{E}_{a}^{*}}{\sqrt{\rm{Var}_{a}^{*}}}$ (and thus $\frac{\rm{E}_{a}-(c_{2}-c_{1})}{\sqrt{\rm{Var}_{a}}}$) go to zero, resulting in a poorly-performing LS-SVM. The same occurs when $c_{1}\to 1$ and $c_{2}\to 0$. Figure 8 collaborates this remark with $c_{1}=1/32$ in subfigure (a) and $1/2$ in (b). Note that in subfigure (a), even with a smartly chosen threshold $\xi$, LS-SVM is impossible to perform as well as in the case $c_{1}=c_{2}$, as a result of the significant overlap between the two histograms.

IV-C Applying to real-world datasets

When the classification performance of real-world datasets is concerned, our theory may be limited by: i) the fact that it is an asymptotic result and allows for an estimation error of order $O(n^{-\frac{1}{2}})$ between theory and practice and ii) the strong Gaussian assumption for the input data.

However, when applied to real-world datasets, here to the popular MNIST [28] and Fashion-MNIST [29] datasets, our asymptotic results, which are theoretically only applicable for Gaussian data, show an unexpectedly similar behavior. Here we consider a two-class classification problem with a training set of $n=256$ vectorized images of size $p=784$ randomly selected from the MNIST and Fashion-MNIST datasets (numbers $8$ and $9$ in both cases as an example). Then a test set of $n_{\rm test}=256$ is used to evaluate the classification performance. Means and covariances are empirically obtained from the full set of $11\,800$ MNIST images ($5\,851$ images of number $8$ and $5\,949$ of number $9$) and of $11\,800$ Fashion-MNIST images ($5\,851$ images of number $8$ and $5\,949$ of number $9$), respectively. Despite the obvious non-Gaussianity as well as the clearly different nature of the input data (from the two datasets), the distribution of $g(\mathbf{x})$ is still surprisingly close to its Gaussian approximation computed from Theorem 2, as shown in Figure 9 and 10 for MNIST and Fashion-MNIST, respectively. In both cases we plot the results from (a) raw images as well as (b) when Gaussian white noise is artificially added to the image vectors.

In Figure 11 we plot the misclassification rate as a function of the decision threshold $\xi$ for MNIST and Fashion-MNIST data (number $8$ and $9$). We observe that although derived from a Gaussian mixture model, the conclusion from Remark 2, Lemma 1 and Corollary 1 that the decision threshold should approximately be $c_{2}-c_{1}$ rather than [math] approximately holds true in both cases.

In Figure 12 and 13 we evaluated the performance of LS-SVM on the MNIST and Fashion-MNIST datasets (with and without noise) as a function of the kernel parameter $\sigma$ of Gaussian kernel $f(x)=\exp(-x/2\sigma^{2})$. Surprisingly, compared to Figure 5, we face the situation where there is little difference in the performance of LS-SVM as soon as $\sigma^{2}$ is away from [math], which likely comes from the fact that the difference in means $\bm{\mu}_{2}-\bm{\mu}_{1}$ is so large that it becomes predominant over the influence of covariances as mentioned in the first paragraph of Section IV-B. This argument is numerically sustained by Table I. The gap between theory and practice observed as $\sigma^{2}\to 0$ is likely a result of the finite $n,p$ (as in Figure 5) rather than of the Gaussian assumption of the input data, since we observe a similar behavior even when Gaussian white noise is added.

V Concluding remarks

In this work, through a performance analysis of LS-SVM for large dimensional data, we reveal the significance of balanced dataset with $c_{1}=c_{2}$, as well as the interplay between the pivotal kernel function $f$ and the statistical structure of the data. The normalized labels $y^{*}_{i}\in\{-1/c_{1},1/c_{2}\}$ are proposed to mitigate the damage of $c_{2}-c_{1}$ in the decision function. We prove the irrelevance of $\gamma$ when it is considered to remain constant in the large $n,p$ regime; however, this argument is not guaranteed to hold true when $\gamma$ scales with $n,p$. Our theoretical results, even though built upon the assumption of Gaussian data, provide similar results when tested on real-world large dimensional datasets, which offers a possible application despite the strong Gaussian assumption in the general context of large scale supervised learning.

The major difference of the present work compared to other theoretical analyses (for example [13]) is that, by studying the rather simple problem of a two-class Gaussian mixture separation with comparably large instance number and data dimension, together with sufficiently smooth kernel function $f$ and regularization parameter $\gamma$ of order $O(1)$, we deduce explicit results for the output of LS-SVM which surprisingly coincide with observations on some large dimensional real-world datasets (including MNIST and beyond) and therefore allowing for novel insights into the behavior of LS-SVM for large dimensional datasets. Of interest to future work is the remark that, unlike in the work of [13] where, in the large $n$ alone asymptotics, $\gamma$ is best scaled large with $n$, in the present large $p$, large $n$ setting, where we demonstrate rate-optimality of LS-SVM for $\gamma=O(1)$. This apparent paradox could be deciphered through the analysis of more advanced (normalized inner product) kernels of the type $f(\mathbf{x}_{i}^{\sf T}\mathbf{x}_{j}/\sqrt{p})$, studied notably in [37], for which we believe that other scalings for $\gamma$ would be optimal; it is also importantly believed that such kernels could lead to improved performances (not in rate, as those are already optimal in the present setting, but possibly in absolute performance). These technically more involved considerations are left for future investigations.

The extension of the present work to the asymptotic performance analysis of the classical SVM requires more efforts since, there, the decision function $g(\mathbf{x})$ depends implicitly (through the solution to a quadratic programming problem) rather than explicitly on the underlying kernel matrix $\mathbf{K}$. Additional technical tools are thus required to cope with this dependence structure.

The link between LS-SVM and extreme learning machine (ELM) was brought to light in [17] and the performance analysis of ELM in large dimension has been investigated in the recent article [38]. Together with these works, we have the possibility to identify the tight but subtle relation between the kernel function and the activation function in the context of some simple structured neural networks. This is notably of interest when the datasets are so large that computing $\mathbf{K}$ and the decision function $g(\mathbf{x})$ becomes prohibitive, a problem largely alleviated by neural networks with controllable number of neurons. This link also generally opens up a possible direction of research into the complex neural networks realm.

Appendix A Proof of Theorem 1

Our key interest here is on the decision function of LS-SVM: $g(\mathbf{x})=\bm{\alpha}^{{\sf T}}\mathbf{k}(\mathbf{x})+b$ with $(\bm{\alpha},b)$ given by

[TABLE]

and $\mathbf{S}^{-1}=\left({\mathbf{K}}+\frac{n}{\gamma}\mathbf{I}_{n}\right)^{-1}$.

Before going into the detailed proof, as we will frequently deal with random variables evolving as $n,p$ grow large, we shall use the extension of the $O(\cdot)$ notation introduced in [20]: for a random variable $x\equiv x_{n}$ and $u_{n}\geq 0$, we write $x=O(u_{n})$ if for any $\eta>0$ and $D>0$, we have $n^{D}{\rm P}(x\geq n^{\eta}u_{n})\to 0$. Note that under Assumption 1 it is equivalent to use either $O(u_{n})$ or $O(u_{p})$ since $n,p$ scales linearly. In the following we shall use constantly $O(u_{n})$ for simplicity.

When multidimensional objects are concerned, $\mathbf{v}=O(u_{n})$ means the maximum entry of a vector (or a diagonal matrix) $\mathbf{v}$ in absolute value is of order $O(u_{n})$ and $\mathbf{M}=O(u_{n})$ means that the operator norm of $\mathbf{M}$ is of order $O(u_{n})$. We refer the reader to [20] for more discussions on these practical definitions.

Under the growth rate settings of Assumption 1, from [20], the approximation of the kernel matrix $\mathbf{K}$ is given by

[TABLE]

with $\beta=f(0)-f(\tau)+\tau f^{\prime}(\tau)$ and $\mathbf{A}=\mathbf{A}_{n}+\mathbf{A}_{\sqrt{n}}+\mathbf{A}_{1}$, $\mathbf{A}_{n}=-\frac{f(\tau)}{2f^{\prime}(\tau)}\mathbf{1}_{n}\mathbf{1}_{n}^{{\sf T}}$ and $\mathbf{A}_{\sqrt{n}}$, $\mathbf{A}_{1}$ given by (18) and (A) at the top of next page, where we denote

[TABLE]

We start with the term $\mathbf{S}^{-1}$. The terms of leading order in $\mathbf{K}$, i.e.,$-2f^{\prime}(\tau)\mathbf{A}_{n}$ and $\frac{n}{\gamma}\mathbf{I}_{n}$ are both of operator norm $O(n)$. Therefore a Taylor expansion can be performed as

[TABLE]

with $\mathbf{L}=\left(f(\tau)\frac{\mathbf{1}_{n}\mathbf{1}_{n}^{{\sf T}}}{n}+\frac{\mathbf{I}_{n}}{\gamma}\right)^{-1}$ of order $O(1)$ and $\mathbf{Q}=\frac{2f^{\prime}(\tau)}{n^{2}}\left(\mathbf{A}_{1}+\mathbf{P}\bm{\Omega}^{{\sf T}}\bm{\Omega}\mathbf{P}+\frac{2f^{\prime}(\tau)}{n}\mathbf{A}_{\sqrt{n}}\mathbf{L}\mathbf{A}_{\sqrt{n}}\right)$.

With the Sherman-Morrison formula we are able to compute explicitly $\mathbf{L}$ as

[TABLE]

Writing $\mathbf{L}$ as a linear combination of $\mathbf{I}_{n}$ and $\mathbf{P}$ is useful when computing $\mathbf{L}\mathbf{1}_{n}$ or $\mathbf{1}_{n}^{{\sf T}}\mathbf{L}$, because by the definition of $\mathbf{P}=\mathbf{I}_{n}-\frac{\mathbf{1}_{n}\mathbf{1}_{n}^{{\sf T}}}{n}$, we have $\mathbf{1}_{n}^{{\sf T}}\mathbf{P}=\mathbf{P}\mathbf{1}_{n}=\mathbf{0}$.

We shall start with the term $\mathbf{1}_{n}^{{\sf T}}\mathbf{S}^{-1}$, since it is the basis of several other terms appearing in $\bm{\alpha}$ and $b$,

[TABLE]

since $\mathbf{1}_{n}^{{\sf T}}\mathbf{L}=\frac{\gamma}{1+\gamma f(\tau)}\mathbf{1}_{n}^{{\sf T}}$.

With $\mathbf{1}_{n}^{{\sf T}}\mathbf{S}^{-1}$ at hand, we next obtain,

[TABLE]

[TABLE]

[TABLE]

The inverse of $\mathbf{1}_{n}^{{\sf T}}\mathbf{S}^{-1}\mathbf{1}_{n}$ can consequently be computed using a Taylor expansion around its leading order, allowing an error term of $O(n^{-\frac{3}{2}})$ as

[TABLE]

Combing (14) with (15) we deduce

[TABLE]

and similarly the following approximation of $b$ as

[TABLE]

where

[TABLE]

which gives the asymptotic approximation of $b$.

Moving to $\bm{\alpha}$, note from (13) that $\mathbf{L}-\frac{\gamma}{1+\gamma f(\tau)}\frac{\mathbf{1}_{n}\mathbf{1}_{n}^{{\sf T}}}{n}=\gamma\mathbf{P}$, and we can thus rewrite:

[TABLE]

At this point, for $\bm{\alpha}=\mathbf{S}^{-1}\left(\mathbf{I}_{n}-\frac{\mathbf{1}_{n}\mathbf{1}_{n}^{{\sf T}}\mathbf{S}^{-1}}{\mathbf{1}_{n}^{{\sf T}}\mathbf{S}^{-1}\mathbf{1}_{n}}\right)\mathbf{y}$, we have

[TABLE]

Here again, we use $\mathbf{1}_{n}^{{\sf T}}\mathbf{L}=\frac{\gamma}{1+\gamma f(\tau)}\mathbf{1}_{n}^{{\sf T}}$ and $\mathbf{L}-\frac{\gamma}{1+\gamma f(\tau)}\frac{\mathbf{1}_{n}\mathbf{1}_{n}^{{\sf T}}}{n}=\gamma\mathbf{P}$, to eventually get

[TABLE]

Note here the absence of a term of order $O(n^{-3/2})$ in the expression of $\bm{\alpha}$ since $\mathbf{P}\mathbf{A}_{\sqrt{n}}\mathbf{P}=0$ from (18).

We shall now work on the vector $\mathbf{k}(\mathbf{x})$ for a new datum $\mathbf{x}$, following the same analysis as in [20] for the kernel matrix $\mathbf{K}$, assuming that $\mathbf{x}\sim\mathcal{N}(\bm{\mu}_{a},\mathbf{C}_{a})$ and recalling the random variables definitions,

[TABLE]

we show that the $j$-th entry of $\mathbf{k}(\mathbf{x})$ can be written as

[TABLE]

Combining (21) and (22), we deduce

[TABLE]

with $\tilde{\mathbf{k}}(\mathbf{x})$ given in (A).

At this point, note that the term of order $O(n^{-\frac{1}{2}})$ in the final object $g(\mathbf{x})=\bm{\alpha}^{{\sf T}}\mathbf{k}(\mathbf{x})+b$ disappears because in both (17) and (23) the term of order $O(n^{-1/2})$ is $\frac{2\gamma}{\sqrt{p}}c_{1}c_{2}f^{\prime}(\tau)(t_{2}-t_{1})$ but of opposite signs. Also, we see that the leading term $c_{2}-c_{1}$ in $b$ will remain in $g(\mathbf{x})$ as stated in Remark 2.

The development of $\mathbf{y}^{{\sf T}}\mathbf{P}\tilde{\mathbf{k}}(\mathbf{x})$ induces many simplifications, since i) $\mathbf{P}\mathbf{1}_{n}=\mathbf{0}$ and ii) random variables as $\bm{\omega}_{\mathbf{x}}$ and $\bm{\psi}$ in $\tilde{\mathbf{k}}(\mathbf{x})$, once multiplied by $\mathbf{y}^{{\sf T}}\mathbf{P}$, thanks to probabilistic averaging of independent zero-mean terms, are of smaller order and thus become negligible. We thus get

[TABLE]

This result, together with (23), completes the analysis of the term $\bm{\alpha}^{{\sf T}}\mathbf{k}(\mathbf{x})$. Combining (23)-(24) with (17) we conclude the proof of Theorem 1.

Appendix B Proof of Theorem 2

This section is dedicated to the proof of the central limit theorem for

[TABLE]

with the shortcut $c_{\mathbf{x}}=-2c_{1}c_{2}^{2}$ for $\mathbf{x}\in\mathcal{C}_{1}$ and $c_{\mathbf{x}}=2c_{1}^{2}c_{2}$ for $\mathbf{x}\in\mathcal{C}_{2}$, and $\mathfrak{P},\mathfrak{D}$ as defined in (7) and (8).

Our objective is to show that for $a\in\{1,2\}$, $n(\hat{g}(\mathbf{x})-G_{a})\overset{d}{\to}0$ with

[TABLE]

where ${\rm E}_{a}$ and ${\rm Var}_{a}$ are given in Theorem 2. We recall that $\mathbf{x}=\bm{\mu}_{a}+\sqrt{p}\bm{\omega}_{\mathbf{x}}$ with $\bm{\omega}_{\mathbf{x}}\sim\mathcal{N}(0,\mathbf{C}_{a}/p)$.

Letting $\mathbf{z}_{\mathbf{x}}$ such that $\bm{\omega}_{\mathbf{x}}=\mathbf{C}_{a}^{1/2}\mathbf{z}_{\mathbf{x}}/\sqrt{p}$, we have $\mathbf{z}_{\mathbf{x}}\sim\mathcal{N}(\mathbf{0},\mathbf{I}_{n})$ and we can rewrite $\hat{g}(\mathbf{x})$ in the following quadratic form (of $\mathbf{z}_{\mathbf{x}}$) as

[TABLE]

with

[TABLE]

Since $\mathbf{z}_{\mathbf{x}}$ is (standard) Gaussian and has the same distribution as $\mathbf{U}\mathbf{z}_{\mathbf{x}}$ for any orthogonal matrix $\mathbf{U}$ (i.e., such that $\mathbf{U}^{{\sf T}}\mathbf{U}=\mathbf{U}\mathbf{U}^{{\sf T}}=\mathbf{I}_{n}$), we choose $\mathbf{U}$ that diagonalize $\mathbf{A}$ such that $\mathbf{A}=\mathbf{U}\bm{\Lambda}\mathbf{U}^{{\sf T}}$, with $\bm{\Lambda}$ diagonal so that $\hat{g}(\mathbf{x})$ and $\tilde{g}(\mathbf{x})$ have the same distribution where

[TABLE]

and $\tilde{}=\mathbf{U}^{{\sf T}}$, $\lambda_{i}$ the diagonal elements of $\bm{\Lambda}$ and $z_{i}$ the elements of $\mathbf{z}_{\mathbf{x}}$.

Conditioning on $\bm{\Omega}$, we thus result in the sum of independent but not identically distributed random variables $r_{i}=z_{i}^{2}\lambda_{i}+z_{i}\tilde{b}_{i}+\frac{c}{n}$. We then resort to the Lyapunov CLT [33, Theorem 27.3].

We begin by estimating the expectation and the variance

[TABLE]

of $r_{i}$, so that

[TABLE]

We shall rewrite $\bm{\Omega}$ into two blocks as:

[TABLE]

where $\mathbf{Z}_{1}\in\mathbb{R}^{p\times n_{1}}$ and $\mathbf{Z}_{2}\in\mathbb{R}^{p\times n_{2}}$ with i.i.d. Gaussian entries with zero mean and unit variance. Then

[TABLE]

and with $\mathbf{P}\mathbf{y}=\mathbf{y}-(c_{2}-c_{1})\mathbf{1}_{n}$, we deduce

[TABLE]

Since $\mathbf{Z}_{i}\mathbf{1}_{n_{i}}\sim\mathcal{N}(\mathbf{0},n_{i}\mathbf{I}_{n_{i}})$, by applying the trace lemma [39, Lemma B.26] we get

[TABLE]

Consider now the events

[TABLE]

for any fixed $\epsilon$ with $\rho=\frac{4nc_{1}^{2}c_{2}^{2}}{p^{2}}\left(\frac{\operatorname{{\rm tr}}\mathbf{C}_{1}\mathbf{C}_{a}}{c_{1}}+\frac{\operatorname{{\rm tr}}\mathbf{C}_{2}\mathbf{C}_{a}}{c_{2}}\right)$ and write

[TABLE]

We start with the variable $\tilde{g}(\mathbf{x})|E$ and check that Lyapunov’s condition for $\bar{r}_{i}=r_{i}-\mathbb{E}[r_{i}]$, conditioning on $E$,

[TABLE]

holds by rewriting

[TABLE]

since both $\lambda_{i}$ and $\tilde{b}_{i}$ are of order $O(n^{-3/2})$.

As a consequence of the above, we have the CLT for the random variable $\tilde{g}(\mathbf{x})|E$, thus

[TABLE]

Next, we see that the second term in (26) goes to zero because \big{|}\mathbb{E}\big{[}\exp\left(iun\frac{\tilde{g}(\mathbf{x})-{\rm E}_{a}}{s_{n}}\right)\big{|}\bar{E}\big{]}\big{|}\leq 1 and ${\rm P}(\bar{E})\to 0$ from (25) and we eventually deduce

[TABLE]

With the help of Lévy’s continuity theorem, we thus prove the CLT of the variable $n\frac{\tilde{g}(\mathbf{x})-{\rm E}_{a}}{s_{n}}$. Since $s_{n}^{2}\to{\rm Var}_{a}$, with Slutsky’s theorem, we have the CLT for $n\frac{\tilde{g}(\mathbf{x})-{\rm E}_{a}}{\sqrt{{\rm Var}_{a}}}$ (thus for $n\frac{\hat{g}(\mathbf{x})-{\rm E}_{a}}{\sqrt{{\rm Var}_{a}}}$), and eventually for $n\frac{{g}(\mathbf{x})-{\rm E}_{a}}{\sqrt{{\rm Var}_{a}}}$ by Theorem 1 which completes the proof.