Least Squares Minimum Class Variance Support Vector Machines

Abstract

In this paper, we propose a Support Vector Machine (SVM)-type algorithm, which is statistically faster among other common algorithms in the family of SVM algorithms. The new algorithm uses distributional information of each class and, therefore, combines the benefits of using the class variance in the optimization with the least squares approach, which gives an analytic solution to the minimization problem and, therefore, is computationally efficient. We demonstrate an important property of the algorithm which allows us to address the inversion of a singular matrix in the solution. We also demonstrate through real data experiments that we improve on the computational time without losing any of the accuracy when compared to previously proposed algorithms.

1. Introduction

Support Vector Machines (SVMs) have been used in a number of disciplines since their introduction by [1] for classification and regression. Although the name is used to describe the classic algorithm proposed by [1], SVMs have been extended in many different directions and a lot of authors are now using the name to refer to a family of methods that are based on the original idea by [1]. The basic purpose of SVM algorithms in binary classification is to find an optimal hyperplane which separates the two classes of datapoints with the maximum margin when the data are separable (hard margin). In cases where the classes are not separable, a soft margin approach is used which finds the optimal hyperplane by maximizing the margin and minimizing the sum of the misclassification distances of the misclassified points.

There are three features that have made SVM algorithms popular since their introduction. The first one is the use of nonlinear kernels, which map the observations from the current space into a higher dimensional feature space to achieve linear separability of the points using the kernel trick, that is without the need to know the exact mapping to the feature space. The second important aspect of SVMs is the fact that they target the minimization of structural risk (minimizing the risk of the misclassification of unseen observations) rather than the minimization of empirical risk (minimizing the risk for points in the sample). Finally, the optimization problem is solved relatively efficiently using quadratic programming.

SVMs are continuing to be a popular option for researchers looking for classification methodology to apply to their datasets. Therefore, there is a constant need for new approaches to be developed within the SVM framework in classification to address the many challenges that the new era of massive and high-dimensional datasets brings to researchers. A very small sample of new methods being proposed in the SVM literature includes [2], which proposes a new fuzzy approach to Twin SVMs; [3], which presents an improved version of the SVM with the radial basis function; and [4], who proposed a twin SVM algorithm with the pinball loss. At the same time, there are some recent works which demonstrate the usefulness in applying SVM variants for classification in other sciences. See, for example, [5] who applied SVMs to Twitter data, the work by the authors in [6] who applied SVMs to remote sensing data, and [7] who applied them in seismic data.

In this paper, we propose a computationally efficient SVM-type algorithm which uses distributional information of the classes. To achieve a computationally fast algorithm, we replace the hinge loss in the classic SVM algorithms with the least squares approach, in a similar way as it was done in the Least Squares SVM (LSSVM) by [8]. To introduce the distributional information of the classes, we propose the use of the within-class variance in a similar way as it was proposed in the Minimum Class Variance SVM (MCVSVM) by [9], who reformulated the optimization problem in Fisher’s linear discriminant analysis (LDA) to achieve this.

In Section 2, we revisit the algorithms in the literature which are important for our development and in Section 3, we propose our new algorithm, presenting both the linear and nonlinear approaches to the algorithm. We demonstrate a very powerful property of the algorithm, which can overcome the issue of finding an inverse by using principal projections, which is needed for the solution in Section 4. We discuss some real data analysis in Section 5 and we close with a discussion section.

2. Literature Review on SVMs

In this section, we review some of the algorithms in the SVM family that were useful for the development of our idea. We start with the classic SVM by [1] and then we present the Least Squares SVM approach by [8]. Finally, we discuss the Minimum Class Variance SVM (MCVSVM) by [9]. All the algorithms were initially developed in the simple case where data are separable, and then extended in the soft margin case where the data are not linearly separable. In this paper, we talk about the most general approach, that is, the soft margin approach.

2.1. Support Vector Machines (SVMs)

The classic SVM algorithm was proposed by [1]. In the most general case, to find the optimal separating hyperplane, it was proposed to solve the following optimization problem:

$$\begin{array}{c}\hfill {\displaystyle min\frac{1}{2}{\mathit{\psi }}^{\mathsf{T}}\mathit{\psi }+\lambda \sum _{i=1}^n{\xi }_i}\end{array}$$

under the constraints:

$$\begin{array}{c}\hfill y_i({\mathit{\psi }}^{\mathsf{T}}{\mathit{x}}_i-t)\ge 1-{\xi }_i,{\xi }_i\ge 0,\end{array}$$

where $(\mathit{\psi },t)\in {\mathbb{R}}^p\times \mathbb{R}$ is the pair that characterizes the hyperplane that has equation ${\mathit{\psi }}^{\mathsf{T}}\mathit{X}-t=0$, $\lambda$ is a scalar known as the cost or the misclassification penalty, and the ${\xi }_i$’s are slack variables which denote misclassification distances. If a point is correctly classified, the slack variable ${\xi }_i$ associated with it is set to 0, and if the point is incorrectly classified, the ${\xi }_i$ denotes the distance of the variable to the to the hyperplane.

To find the solution to the optimization problem above, one uses the Lagrangian multipliers and tries to minimize the following Lagrangian equation:

$$\begin{array}{cc}\hfill L(\mathit{\psi },t,\mathit{\xi },\mathit{\alpha },\mathit{\beta })=& {\displaystyle \frac{1}{2}{\mathit{\psi }}^{\mathsf{T}}\mathit{\psi }+\lambda \sum _{i=1}^n{\xi }_i-\sum _{i=1}^n{\mathit{\alpha }}_i(1-{\xi }_i-y_i({\mathit{\psi }}^{\mathsf{T}}{\mathit{x}}_i-t))-\sum _{i=1}^n{\beta }_i{\xi }_i}\hfill \\ \hfill =& \frac{1}{2}{\mathit{\psi }}^{\mathsf{T}}\mathit{\psi }+\lambda {\mathbf{1}}^{\mathsf{T}}\mathit{\xi }-{\mathit{\alpha }}^{\mathsf{T}}(\mathbf{1}-\mathit{\xi }-\mathit{y}\odot ({\mathit{\psi }}^{\mathsf{T}}\mathit{x}-\mathit{t}))-{\mathit{\beta }}^{\mathsf{T}}\mathit{\xi }\hfill \end{array}$$

where $\mathit{y}={(y_1,\dots ,y_n)}^{\mathsf{T}}\in {\mathbb{R}}^n$, $\mathit{x}={({\mathit{x}}_1,\dots ,{\mathit{x}}_n)}^{\mathsf{T}}\in {\mathbb{R}}^{p\times n}$, $\mathbf{1}={(1,\dots ,1)}^{\mathsf{T}}\in {\mathbb{R}}^n$, $\mathit{t}={(t,\dots ,t)}^{\mathsf{T}}\in {\mathbb{R}}^n$, $\mathit{\alpha }={({\alpha }_1,\dots ,{\alpha }_n)}^{\mathsf{T}}\in {\mathbb{R}}^n$, $\mathit{\beta }={({\beta }_1,\dots ,{\beta }_n)}^{\mathsf{T}}\in {\mathbb{R}}^n$, and $\mathit{\xi }={({\xi }_1,\dots ,{\xi }_n)}^{\mathsf{T}}\in {\mathbb{R}}^n$. Using the derivatives, one finds the Karush–Kuhn–Tucker (KKT) equations:

$$\begin{array}{cc}\hfill \frac{\partial L}{\partial \mathit{\psi }}=& \mathit{\psi }-\mathit{x}(\mathit{\alpha }\odot \mathit{y})=0\Rightarrow \mathit{\psi }=\mathit{x}(\mathit{\alpha }\odot \mathit{y}),\hfill \\ \hfill \frac{\partial L}{\partial t}=& {\mathit{\alpha }}^{\mathsf{T}}\mathit{y}=0\Rightarrow {\mathit{\alpha }}^{\mathsf{T}}\mathit{y}=0,\hfill \\ \hfill \frac{\partial L}{\partial \xi }=& \lambda -\mathit{\alpha }-\mathit{\beta }=0.\hfill \end{array}$$

By replacing the result of the first derivative in the Lagrangian equation, one gets:

$$\begin{array}{c}\hfill L(\mathit{\psi },t,\mathit{\xi },\mathit{\alpha },\mathit{\beta })=\frac{1}{2}{\left(\mathit{x}(\mathit{\alpha }\odot \mathit{y})\right)}^{\mathsf{T}}\mathit{x}(\mathit{\alpha }\odot \mathit{y})+\lambda {\mathbf{1}}^{\mathsf{T}}\mathit{\xi }-{\mathit{\alpha }}^{\mathsf{T}}(\mathbf{1}-\mathit{\xi }-\mathit{y}\odot ({\left(\mathit{x}(\mathit{\alpha }\odot \mathit{y})\right)}^{\mathsf{T}}\mathit{x}-\mathit{t}))-{\mathit{\beta }}^{\mathsf{T}}\mathit{\xi }\\ \hfill =\frac{1}{2}{(\mathit{\alpha }\odot \mathit{y})}^{\mathsf{T}}{\mathit{x}}^{\mathsf{T}}\mathit{x}(\mathit{\alpha }\odot \mathit{y})+{\mathit{\alpha }}^{\mathsf{T}}\mathbf{1}-{(\lambda \mathbf{1}-\mathit{\alpha }-\mathit{\beta })}^{\mathsf{T}}\mathit{\xi }+{(\mathit{\alpha }\odot \mathit{y})}^{\mathsf{T}}{\mathit{x}}^{\mathsf{T}}\mathit{x}(\mathit{\alpha }\odot \mathit{y})+t{\mathit{\alpha }}^{\mathsf{T}}\mathit{y}.\end{array}$$

Now, the last term is equal to 0 from the result in the second KKT equation above, and the third term is 0 by the third KKT equation. Hence, the above Lagrangian equation reduces to:

$$\begin{array}{c}\hfill L(\mathit{\alpha },\mathit{\beta })={\mathit{\alpha }}^{\mathsf{T}}\mathbf{1}-\frac{1}{2}{(\mathit{\alpha }\odot \mathit{y})}^{\mathsf{T}}{\mathit{x}}^{\mathsf{T}}\mathit{x}(\mathit{\alpha }\odot \mathit{y})\end{array}$$

subject to the constraint $\mathbf{0}<\mathit{\alpha }<\lambda \mathbf{1}$. This is known as the dual problem. One can use quadratic programming optimization to solve the dual problem to obtain $\mathit{\alpha }$, which is essential in estimating the normal vector of the optimal hyperplane according to the solution in the first KKT equation above.

Two of the most important features of the classic SVM algorithm is the use of the hinge loss in the objective function to be minimized and the fact that, in constructing the hyperplane, only points that are incorrectly classified are used, as well as points that are closer to the hyperplane. This gives a form of sparsity to the SVM, as not all points are needed to construct the hyperplane. At the same time, it is computationally expensive as it requires the solution of a quadratic programming optimization problem.

2.2. Least Squares Support Vector Machines (LSSVMs)

A least squares approach was proposed by [8], which essentially changes the geometry of the problem from the way [1] framed it. It changes the hinge loss to the square loss, which allows one to take a least squares approach and find the solution analytically. At the same time, they changed the constraints from inequalities to equalities, allowing for the ${\xi }_i$’s to be either positive or negative, which implies that all the points are needed to find he optimal hyperplane, removing the sparsity of the classic SVM algorithm.

In the Least Squares SVM (LSSVM), one tries to find the optimal hyperplane by minimizing the following objective function:

$$\begin{array}{c}\hfill {\displaystyle min\frac{1}{2}{\mathit{\psi }}^{\mathsf{T}}\mathit{\psi }+\lambda \sum _{i=1}^n{\xi }_i^2}\end{array}$$

under the constraints:

$$\begin{array}{c}\hfill y_i({\mathit{\psi }}^{\mathsf{T}}{\mathit{x}}_i-t)=1-{\xi }_i.\end{array}$$

By replacing the equality constraint in the objective function, one gets that is needed to solve the objective function:

$$\begin{array}{c}\hfill {\displaystyle min\frac{1}{2}{\mathit{\psi }}^{\mathsf{T}}\mathit{\psi }+\lambda \sum _{i=1}^n{\left[1-y_i({\mathit{\psi }}^{\mathsf{T}}{\mathit{x}}_i-t)\right]}^2.}\end{array}$$

To find the values of the pair $(\mathit{\psi },t)$ which minimizes this objective function, one needs to take the derivative with respect to both parameters. Here, we rewrite the objective function by denoting $\mathit{r}={({\mathit{\psi }}^{\mathsf{T}},t)}^{\mathsf{T}}$ and reexpress the above as:

$$\begin{array}{c}\hfill {\displaystyle {\mathit{r}}^{\mathsf{T}}{\mathit{I}}_{p,1}^{\ast }\mathit{r}+\lambda \sum _{i=1}^n{\left[1-y_i\left({\mathit{r}}^{\mathsf{T}}{\mathit{x}}_i^{\ast }\right)\right]}^2}\end{array}$$

where ${\mathit{x}}_i^{\ast }=({\mathit{x}}_i^{\mathsf{T}},-1)$ and ${\mathit{I}}_{p,q}^{\ast }$ is the $(p+q)\times (p+q)$ diagonal matrix which has 1 on the first p diagonal elements and 0 on the last q diagonal elements. In matrix form, we can write this as:

$$\begin{array}{c}\hfill {\mathit{r}}^{\mathsf{T}}{\mathit{I}}_{p,1}^{\ast }\mathit{r}+\lambda {({\mathbf{1}}_n-{\mathit{D}}_{\mathit{Y}}{\mathit{X}}^{\ast }\mathit{r})}^{\mathsf{T}}({\mathbf{1}}_n-{\mathit{D}}_{\mathit{Y}}{\mathit{X}}^{\ast }\mathit{r})\end{array}$$

where ${\mathit{X}}^{\ast }=(\mathit{X},-{\mathbf{1}}_n)$ is the $n\times (p+1)$ matrix which contains the variables $\mathit{X}$ and an extra column of $-1$’s, and ${\mathit{D}}_{\mathit{Y}}$ is the diagonal matrix that has the vector $\mathit{Y}=(Y_1,\dots ,Y_n)$ on the diagonal. Taking the derivative we have:

$$\begin{array}{c}\hfill 2{\mathit{I}}_{p,1}^{\ast }\mathit{r}-2\lambda {\left({\mathit{X}}^{\ast }\right)}^{\mathsf{T}}{\mathit{D}}_{\mathit{Y}}({\mathbf{1}}_n-{\mathit{D}}_{\mathit{Y}}{\mathit{X}}^{\ast }\mathit{r})\end{array}$$

which, if we set it to equal 0, gives the solution:

$$\begin{array}{c}\hfill \mathit{r}={\left(\frac{{\mathit{I}}_{p,1}^{\ast }}{\lambda }+{\left({\mathit{X}}^{\ast }\right)}^{\mathsf{T}}{\mathit{X}}^{\ast }\right)}^{-1}{\left({\mathit{X}}^{\ast }\right)}^{\mathsf{T}}{\mathit{D}}_{\mathit{Y}}{\mathbf{1}}_n\end{array}$$

where we use the fact that ${\mathit{D}}_{\mathit{Y}}{\mathit{D}}_{\mathit{Y}}={\mathit{I}}_n$, the $n\times n$ identity matrix, to simplify the notation. We are mostly interested in $\mathit{\psi }={\left[\mathit{r}\right]}_p$ where ${[·]}_p$ denotes the first p entries in a vector or the first p rows of a matrix, depending on the type of argument being used.

As one can see from the developments in this section, LSSVM does not need quadratic optimization as it has an analytic solution and it is therefore much faster to find the equation of the optimal hyperplane.

2.3. Minimum Class Variances Support Vector Machines (MCVSVMs)

A different extension of the SVM algorithm was proposed by [9], which is called Minimum Class Variance SVM (MCVSVM). As its name suggest, this algorithm is focused on finding the optimal hyperplane, not only by maximizing the margin, but also by minimizing the class variance when projected on the normal vector. This method was inspired by Fisher’s linear discriminant analysis [10], as it uses information from the distribution of the classes to achieve better classification results. It is also interesting that the authors showed that their approach can be used in a large p small n setting, despite using the inverse of the pooled class covariance matrix (which is not invertible in large p small n settings) in their solution.

In the Minimum Class Variance SVM (MCVSVM), one tries to find the optimal hyperplane by minimizing the following objective function:

$$\begin{array}{c}\hfill {\displaystyle min\frac{1}{2}{\mathit{\psi }}^{\mathsf{T}}{\mathbf{\Sigma }}_w\mathit{\psi }+\lambda \sum _{i=1}^n{\xi }_i}\end{array}$$

under the constraints:

$$\begin{array}{c}\hfill y_i({\mathit{\psi }}^{\mathsf{T}}{x}_i-t)\ge 1-{\xi }_i,{\xi }_i\ge 0\end{array}$$

where ${\mathbf{\Sigma }}_w$ is the pooled covariance matrix as a weighted average of the covariance matrices of the classes.

Using the KKT equations, the solution of the hyperplane is found using:

$$\begin{array}{c}\hfill \mathit{\psi }={\mathbf{\Sigma }}_w^{-1}\mathit{x}(\mathit{\alpha }\odot \mathit{y})\end{array}$$

3. New Method: LS-MCVSVM

As we said in the introduction, our objective in this section is to introduce the Least Squares extension of the MCVSVM algorithm. This will allow us to utilize the advantages of both algorithms in an effort to create a much broader algorithm than the already existing one. First of all, the new algorithm is computationally very fast as it does not need quadratic optimization due to the use of the least squares approach to the MCVSVM algorithm. It is also an algorithm that utilizes the variability within each class due to the the use of the MCVSVM.

Therefore, the new algorithm minimizes the following objective function:

$$\begin{array}{c}\hfill {\displaystyle min\frac{1}{2}{\mathit{\psi }}^{\mathsf{T}}{\mathbf{\Sigma }}_w\mathit{\psi }+\lambda \sum _{i=1}^n{\xi }_i^2}\end{array}$$

(1)

under the equality constraints:

$$\begin{array}{c}\hfill y_i({\mathit{\psi }}^{\mathsf{T}}{\mathit{x}}_i-t)=1-{\xi }_i,\end{array}$$

(2)

Using a similar approach as the LSSVM which we described in the review in the previous section, we replace the equality constraint in the objective function to get the new objective functions:

$$\begin{array}{c}\hfill {\displaystyle min\frac{1}{2}{\mathit{\psi }}^{\mathsf{T}}{\mathbf{\Sigma }}_w\mathit{\psi }+\lambda \sum _{i=1}^n{\left[1-y_i({\mathit{\psi }}^{\mathsf{T}}{\mathit{x}}_i-t)\right]}^2}\end{array}$$

We set $\mathit{r}={({\mathit{\psi }}^{\mathsf{T}},t)}^{\mathsf{T}}$ and we rewrite the optimization function as:

$$\begin{array}{c}\hfill {\displaystyle {\mathit{r}}^{\mathsf{T}}\mathit{r}+\lambda \sum _{i=1}^n{\left[1-y_i\left({\mathit{r}}^{\mathsf{T}}{\mathit{x}}_i^{\ast }\right)\right]}^2}\end{array}$$

where ${\mathbf{\Sigma }}_w^{\ast }=\mathrm{diag}({\mathbf{\Sigma }}_w,0)$ a $(p+1)\times (p+1)$ matrix and ${\mathit{x}}_i^{\ast }=({\mathit{x}}_i^{\mathsf{T}},-1)$. We can also write this in matrix form as:

$$\begin{array}{c}\hfill {\mathit{r}}^{\mathsf{T}}{\mathbf{\Sigma }}_w^{\ast }\mathit{r}+\lambda {({\mathbf{1}}_n-{\mathit{D}}_{\mathit{Y}}{\mathit{X}}^{\ast }\mathit{r})}^{\mathsf{T}}({\mathbf{1}}_n-{\mathit{D}}_{\mathit{Y}}{\mathit{X}}^{\ast }\mathit{r})\end{array}$$

where ${\mathit{X}}^{\ast }=(\mathit{X},-{\mathbf{1}}_n)$ is the $n\times (p+1)$ matrix which contains the variables $\mathit{X}$ and an extra column of $-1$’s, and ${\mathit{D}}_{\mathit{Y}}$ is the diagonal matrix that has the vector $\mathit{Y}=(Y_1,\dots ,Y_n)$ on the diagonal. Now, if one takes the derivative and set it equal to 0, the solution is as follows:

$$\begin{array}{c}\hfill \mathit{r}={\left(\frac{{\mathbf{\Sigma }}_w^{\ast }}{\lambda }+{\left({\mathit{X}}^{\ast }\right)}^{\mathsf{T}}{\mathit{X}}^{\ast }\right)}^{-1}{\left({\mathit{X}}^{\ast }\right)}^{\mathsf{T}}{\mathit{D}}_{\mathit{Y}}{\mathbf{1}}_n\end{array}$$

where, as before, $\mathit{r}={({\mathit{\psi }}^{\mathsf{T}},t)}^{\mathsf{T}}$ and ${\mathbf{\Sigma }}_w^{\ast }$ is a $(p+1)\times (p+1)$ matrix which has ${\mathbf{\Sigma }}_w$ in the first $p\times p$ submatrix and everything else is completed with zeroes. We omitted the details of the development as it is very similar to the one described in the LSSVM above.

It is also important to note that there are similar developments in the nonlinear setting. Let $\varphi$ be a functions such that $\varphi :{\mathbb{R}}^p\to {\mathbb{R}}^q$ where $q>>p$ is the dimension of the feature space where the points are mapped to be separated linearly. Then, we can define the within sample variance in the ${\mathbf{\Sigma }}_w^{\mathbf{\Phi }}$ in the feature space as:

$$\begin{array}{c}\hfill {\mathbf{\Sigma }}_w^{\mathbf{\Phi }}=\sum _{\mathit{x}\in C-}(\varphi \left(\mathit{x}\right)-{\mathit{\mu }}_{C-}^{\mathbf{\Phi }}){(\varphi \left(\mathit{x}\right)-{\mathit{\mu }}_{C-}^{\mathbf{\Phi }})}^{\mathsf{T}}+\sum _{\mathit{x}\in C+}(\varphi \left(\mathit{x}\right)-{\mathit{\mu }}_{C+}^{\mathbf{\Phi }}){(\varphi \left(\mathit{x}\right)-{\mathit{\mu }}_{C+}^{\mathbf{\Phi }})}^{\mathsf{T}}\end{array}$$

where $C-,C+$ denote the points in each class and ${\mathit{\mu }}_{C-}^{\mathbf{\Phi }},{\mathit{\mu }}_{C+}^{\mathbf{\Phi }}$ the means of the predictor vectors when transformed by $\varphi$ to the feature space.

This means that the optimization problem we are solving involves the minimization of the following objective function:

$$\begin{array}{c}\hfill {\displaystyle min\frac{1}{2}{\mathit{\psi }}_{\mathbf{\Phi }}^{\mathsf{T}}{\mathbf{\Sigma }}_w^{\mathbf{\Phi }}{\mathit{\psi }}_{\mathbf{\Phi }}+\lambda \sum _{i=1}^n{\xi }_i^2}\end{array}$$

under the equality constraints:

$$\begin{array}{c}\hfill y_i({\mathit{\psi }}_{\mathbf{\Phi }}^{\mathsf{T}}\varphi \left({\mathit{X}}_i\right)-t_{\mathbf{\Phi }})=1-{\xi }_i,\end{array}$$

which will give us the solution:

$$\begin{array}{c}\hfill {\mathit{r}}^{\mathbf{\Phi }}={\left(\frac{{\left({\mathbf{\Sigma }}_w^{\mathbf{\Phi }}\right)}^{\ast }}{\lambda }+{\varphi }^{\ast }{\left(\mathit{X}\right)}^{\mathsf{T}}{\varphi }^{\ast }\left(\mathit{X}\right)\right)}^{-1}{\left({\varphi }^{\ast }\left(\mathit{X}\right)\right)}^{\mathsf{T}}{\mathit{D}}_{\mathit{Y}}{\mathbf{1}}_n\end{array}$$

where ${\mathit{r}}^{\mathbf{\Phi }}={({\mathit{\psi }}_{\mathbf{\Phi }}^{\mathsf{T}},t_{\mathbf{\Phi }})}^{\mathsf{T}}$, ${\left({\mathbf{\Sigma }}_w^{\mathbf{\Phi }}\right)}^{\ast }=\mathrm{diag}({\mathbf{\Sigma }}_w^{\mathbf{\Phi }},0)$, and ${\varphi }^{\ast }\left(\mathit{X}\right)={(\varphi {\left(\mathit{X}\right)}^{\mathsf{T}},1)}^{\mathsf{T}}$.

4. Addressing Singularity Using Principal Projections

As one can see in the discussion in the previous section, in order to solve the optimization problem and find $\mathit{r}$, we need the inverse matrix of:

$$\begin{array}{c}\hfill \mathit{A}=\frac{{\mathbf{\Sigma }}_w^{\ast }}{\lambda }+{\left({\mathit{X}}^{\ast }\right)}^{\mathsf{T}}{\mathit{X}}^{\ast }\end{array}$$

which may not be invertible. In this section, we will try to address the possible singularity of this matrix and demonstrate how one can overcome this difficulty using principal projections.

We first assume that the eigenvalues of $\mathit{A}$ form an orthonormal basis. Then, we can define the space $\mathcal{A}$ spanned by the eigenvectors corresponding to the nonzero eigenvalues of $\mathit{A}$ and the space ${\mathcal{A}}^{\perp }$ spanned by the eigenvectors corresponding to the zero eigenvalues of $\mathit{A}$. Therefore, we can write each vector in a $(p+1)$-dimensional space as $\mathit{r}=\mathit{\varphi }+\mathit{\zeta }$, where $\mathit{\varphi }\in \mathcal{A}$ and $\mathit{\zeta }\in {\mathcal{A}}^{\perp }$.

We also note that the optimization problem in (1) alongside the constraint in (2) can be rewritten as:

$$\begin{array}{c}\hfill {\displaystyle min\frac{1}{2}{\mathit{r}}^{\mathsf{T}}{\mathbf{\Sigma }}_w^{\ast }\mathit{r}+\lambda \sum _{i=1}^n{\xi }_i^2}\end{array}$$

under the equality constraints:

$$\begin{array}{c}\hfill y_i\left({\mathit{r}}^{\mathsf{T}}{\mathit{x}}_i^{\ast }\right)=1-{\xi }_i,\end{array}$$

and when we replace ${\xi }_i$ from the constraint to the optimization, the above simplifies to:

$$\begin{array}{c}\hfill {\displaystyle min\frac{1}{2}{\mathit{r}}^{\mathsf{T}}{\mathbf{\Sigma }}_w^{\ast }\mathit{r}+\lambda \sum _{i=1}^n{(1-y_i\left({\mathit{r}}^{\mathsf{T}}{\mathit{x}}_i^{\ast }\right))}^2}\end{array}$$

which, in matrix form, looks like:

$$\begin{array}{c}\hfill min\frac{1}{2}{\mathit{r}}^{\mathsf{T}}\left({\mathbf{\Sigma }}_w^{\ast }+\lambda {\left({\mathit{X}}^{\ast }\right)}^{\mathsf{T}}\left({\mathit{X}}^{\ast }\right)\right)\mathit{r}+\lambda n+2\lambda {\mathit{r}}^{\mathsf{T}}{\left({\mathit{X}}^{\ast }\right)}^{\mathsf{T}}{\mathit{D}}_Y\mathbf{1}\end{array}$$

which can also be rewritten as:

$$\begin{array}{c}\hfill min\frac{1}{2}{\mathit{r}}^{\mathsf{T}}\left(\frac{{\mathbf{\Sigma }}_w^{\ast }}{\lambda }+{\left({\mathit{X}}^{\ast }\right)}^{\mathsf{T}}\left({\mathit{X}}^{\ast }\right)\right)\mathit{r}+n-2{\mathit{r}}^{\mathsf{T}}{\left({\mathit{X}}^{\ast }\right)}^{\mathsf{T}}{\mathit{D}}_Y\mathbf{1}\end{array}$$

where, if we replace the definition of $\mathit{A}$, we get:

$$\begin{array}{c}\hfill min\frac{1}{2}{\mathit{r}}^{\mathsf{T}}\mathit{A}\mathit{r}+n-2{\mathit{r}}^{\mathsf{T}}{\left({\mathit{X}}^{\ast }\right)}^{\mathsf{T}}{\mathit{D}}_Y\mathbf{1}.\end{array}$$

(3)

From this, we can see that by essentially replacing $\mathit{r}$ with $\mathit{\varphi }+\mathit{\zeta }$ we have:

$$\begin{array}{c}\hfill min\frac{1}{2}{\mathit{\varphi }}^{\mathsf{T}}\mathit{A}\mathit{\varphi }+{\mathit{\zeta }}^{\mathsf{T}}\mathit{A}\mathit{\zeta }+n-2{\mathit{\varphi }}^{\mathsf{T}}{\left({\mathit{X}}^{\ast }\right)}^{\mathsf{T}}{\mathit{D}}_Y\mathbf{1}-2{\mathit{\zeta }}^{\mathsf{T}}{\left({\mathit{X}}^{\ast }\right)}^{\mathsf{T}}{\mathit{D}}_Y\mathbf{1}.\end{array}$$

where the term ${\mathit{\zeta }}^{\mathsf{T}}\mathit{A}\mathit{\zeta }=0$, since $\mathit{A}\mathit{\zeta }=0$ as $\mathit{\zeta }\in \mathcal{A}$.

Furthermore, it is important to note that if ${\mathit{\zeta }}^{\mathsf{T}}\mathit{A}\mathit{\zeta }=0$ then, because both $\frac{{\mathbf{\Sigma }}_w^{\ast }}{\lambda }$ and ${\left({\mathit{X}}^{\ast }\right)}^{\mathsf{T}}\left({\mathit{X}}^{\ast }\right)$ are nonnegative matrices, one can show that ${\mathit{\zeta }}^{\mathsf{T}}\frac{{\mathbf{\Sigma }}_w^{\ast }}{\lambda }\mathit{\zeta }=0$ and ${\mathit{\zeta }}^{\mathsf{T}}{\left({\mathit{X}}^{\ast }\right)}^{\mathsf{T}}\left({\mathit{X}}^{\ast }\right)\mathit{\zeta }=0$. From the latter, one can infer that all points ${\mathit{x}}_i^{\ast }$ are projected on the same point under $\mathit{\zeta }$, which leads to the fact that ${\mathit{\zeta }}^{\mathsf{T}}{\left({\mathit{X}}^{\ast }\right)}^{\mathsf{T}}=k$, which makes the last term a constant which can be ignored. Therefore, the optimization problem (3) is equivalent to:

$$\begin{array}{c}\hfill min\frac{1}{2}{\mathit{\varphi }}^{\mathsf{T}}\mathit{A}\mathit{\varphi }+n-2{\mathit{\varphi }}^{\mathsf{T}}{\left({\mathit{X}}^{\ast }\right)}^{\mathsf{T}}{\mathit{D}}_Y\mathbf{1}.\end{array}$$

(4)

Now, this means we can solve the problem in a space isomorphic to $\mathcal{A}$ and, essentially, we can choose to do it on the space spanned by the eigenvectors corresponding to the nonzero eigenvalues of $\mathit{A}$. These can be found using the matrix $(p+1)\times d$ matrix $\mathit{P}$ (where d is the number of nonzero eigenvalues of $\mathit{A}$), which has for columns the eigenvectors corresponding to the nonzero eigenvalues. This means the data can be projected to the new data ${\mathit{X}}^{†}={\mathit{X}}^{\ast }\mathit{P}$. Similarly, we can project $\mathit{\varphi }$ to get $\mathit{\eta }={\mathit{P}}^{\mathsf{T}}\mathit{\varphi }$. Therefore, we can show that the objective function (4) is equivalent to:

$$\begin{array}{c}\hfill min\frac{1}{2}{\mathit{\eta }}^{\mathsf{T}}{\mathit{P}}^{\mathsf{T}}\mathit{A}\mathit{P}\mathit{\eta }+n-2{\mathit{\eta }}^{\mathsf{T}}{\mathit{P}}^{\mathsf{T}}{\left({\mathit{X}}^{\ast }\right)}^{\mathsf{T}}{\mathit{D}}_Y\mathbf{1}.\end{array}$$

which simplifies to:

$$\begin{array}{c}\hfill min\frac{1}{2}{\mathit{\eta }}^{\mathsf{T}}{\mathit{A}}^{†}\mathit{\eta }+n-2{\mathit{\eta }}^{\mathsf{T}}{\left({\mathit{X}}^{†}\right)}^{\mathsf{T}}{\mathit{D}}_Y\mathbf{1}.\end{array}$$

(5)

where ${\mathit{A}}^{†}={\mathit{P}}^{\mathsf{T}}\mathit{A}\mathit{P}={\mathit{P}}^{\mathsf{T}}\frac{{\mathbf{\Sigma }}_w^{\ast }}{\lambda }\mathit{P}+{\mathit{P}}^{\mathsf{T}}{\left({\mathit{X}}^{\ast }\right)}^{\mathsf{T}}{\mathit{X}}^{\ast }\mathit{P}$ which is a $d\times d$ matrix and is also equal to $\frac{{\mathbf{\Sigma }}_w^{†}}{\lambda }+{\left({\mathit{X}}^{†}\right)}^{\mathsf{T}}{\mathit{X}}^{†}$. That is, ${\mathbf{\Sigma }}_w^{†}$ is the within variance when ${\mathit{X}}^{\ast }$ is replaced with ${\mathit{X}}^{†}$. Therefore, we have a projection problem which uses the projected data ${\mathit{X}}^{†}$ in the lower dimensional space (dimension d) only and it is equivalent to the original problem in (5). The solution of this problem is:

$$\begin{array}{c}\hfill \mathit{\eta }={\left({\mathit{A}}^{†}\right)}^{-1}{\left({\mathit{X}}^{†}\right)}^{\mathsf{T}}{\mathit{D}}_{\mathit{Y}}{\mathbf{1}}_n\end{array}$$

which uses the inverse of the ${\mathit{A}}^{†}$, which is nonsingular by construction.

5. Real Data Experiments

To demonstrate the performance of the new algorithm, we ran an analysis on eight datasets. All eight datasets are from the UCI Machine Learning repository. Since we do not discuss multicategory SVM approaches in this paper, we have chosen datasets that have only two classes, or, in the case of multiple classes, we merged together all the classes but the first to create two classes. The datasets we used are summarized in

Table 1. Dataset description. All source links start with ‘https://archive.ics.uci.edu/dataset’ and were valid on the final access on 25 January 2024.

Dataset	Observations	Features	Source, Citation
Iris	150	4	/53/iris, [11]
Haberman’s	306	3	/43/haberman+s+survival, [12]
Ionosphere	351	34	/52/ionosphere, [13]
Breast Cancer	699	9	/15/breast+cancer+wisconsin+original, [14]
Diabetes	≈253,000	21	/891/cdc+diabetes+health+indicators, [15]
Fertility	100	10	/244/fertility, [16]
Seeds	210	7	/dataset/236/seeds, [17]
Banknote	1372	5	/267/banknote+authentication, [18]

.

We split the data into 60% training, 20% testing, and 20% validation datasets and we reported the misclassification rates in the validation dataset, where we use the linear kernel to find an optimal hyperplane and calculate the quantities. The reported quantities in this paper are the average of 10 iterations.

Table 2. Overall misclassification errors (standard errors) for each algorithm in each dataset. The best algorithm for each dataset is highlighted in bold.

	Algorithm	SVM	LS SVM	MCV SVM	LSMCV SVM
Datasets
Iris		0 (0)	0 (0)	0 (0)	0 (0)
Haberman’s		0.27 (0.002)	0.25 (0.002)	0.27 (0.002)	0.25 (0.002)
Ionosphere		0.14 (0.002)	0.14 (0.002)	0.14 (0.002)	0.14 (0.002)
Breast Cancer		0.04 (0.000)	0.04 (0.000)	0.04 (0.000)	0.04 (0.000)
Diabetes		0.24 (0.001)	0.23 (0.001)	0.24 (0.001)	0.23 (0.001)
Fertility		0.17 (0.010)	0.12 (0.005)	0.17 (0.009)	0.12 (0.005)
Seeds		0.07 (0.001)	0.06 (0.001)	0.04 (0.001)	0.03 (0.001)
Banknote		0.01 (0.000)	0.03 (0.000)	0.01 (0.000)	0.03 (0.000)

summarizes the mean misclassification rates and the standard errors. We see that the four algorithms are relatively close, and the least squares approaches (either the classic LSSVM or our proposed methodology, which combines the least squares approach with the minimum class variance) to performed slightly better in most cases.

In addition to the misclassification rate, we calculated the average value of the precision, the recall, and the F1 score in each dataset. We present the results in Figure 1, Figure 2 and Figure 3. As we can see, the performance of the algorithms is very similar across the different metrics. In some cases, our method performs better than the rest (i.e., diabetes and seeds datasets) and in some cases, not as well (i.e., fertility). The differences are small, with our method being very close to the LSSVM performance. In the seeds dataset, our method is better than the LSSVM algorithm, but, in that case, it is very close to the MCVSVM algorithm. This is another indication that our method is able to simultaneously capture the advantages that both the LSSVM and MCVSVM algorithms offer.

Most importantly, we can see in

Table 3. Duration of each algorithm for each dataset (in seconds)—the faster algorithm for each dataset is highlighted in bold.

	Algorithm	SVM	LS SVM	MCV SVM	LSMCV
Datasets
Iris		0.229	0.00931	0.205	0.0207
Haberman’s		0.953	0.0305	0.783	0.0334
Ionosphere		1.29	0.0593	1.12	0.0784
Breast Cancer		8.79	0.241	6.34	0.217
Diabetes		10.4	0.342	8.003	0.307
Fertility		0.139	0.00786	0.125	0.016
Seeds		0.441	0.0164	0.372	0.0224
Banknote		59.3	1.67	49.1	1.42

that both the least square approaches are significantly faster. This difference is actually statistically significant. To demonstrate this, we ran two sample paired nonparametric tests, i.e., Wilcoxon signed-rank tests, for all six pairs of algorithms. The comparison between our algorithm, LSMCVSVM, and SVM gives a p-value of 0.0078, and the comparison between LSMCVSVM with MCVSVM gives a similar p-value, i.e., 0.0078. The comparison between LSSVM and LSMCVSVM gives a nonsignificant p-value (0.7422), which is expected as both algorithms use the least squares approach and have similar running times. The computational gains can be further signified if one extrapolates this difference in massive datasets where we may have a few million datapoints.

For the interested reader, the codes are available at [19].

6. Conclusions

In this work, we presented a new algorithm for classification which combines two existing algorithms. The first algorithm used is the LSSVM, which is one of the fastest algorithms in the SVM family of algorithms, as it has analytical solution. The second algorithm is the MCVSVM, which is an algorithm that generalizes better than the classic SVM algorithm, and it also allows for the variability in each class to be taken into account. The new algorithm, called LSMCVSVM, as demonstrated in our numerical section, has comparable performance with other classic SVM algorithms, like the SVM, LSSVM, and MCVSVM, but it runs in a fraction of time due to the fact that there is no need to solve a quadratic programming optimization problem as one can find an analytic solution. The computational gains of the new algorithms are similar to the computational gains of the LSSVM and the performance is very similar to the MCVSVM, demonstrating that the combination of the two algorithms creates a new algorithm which also combines the advantages of the two algorithms.

Another important aspect of the new method and an important contribution of this paper is the use of principal projections to address the singularity in the solution of $\mathit{r}$. Since the solution of LSMCVSVM requires the use of the inverse of a matrix, without this equivalence, we would not have been able to apply this algorithm to large p small n problems. This provides a way to do this, without the need to introduce other difficult and time-consuming methods for inverse matrix approximation to find the singularity in the matrix.

Other Approaches and Future Work

The SVM literature is full of variants of the classic SVM algorithm since its introduction by [1]. One can combine any of the existing algorithms and create new algorithms, which can be extremely valuable tools in the classification framework. For example, one of the many ideas that can be implemented is the combination of our algorithm with the two-cost alternative, which is an idea used to handle imbalanced classes, i.e., problems where one class has a lot more points than the other class. In this case, it makes sense that a misclassification from the small class should be more costly. Therefore, ref. [20] proposed the use of two different costs or penalties. By giving a bigger penalty to the smaller class, we try to minimize the effect of misclassifying one point may have. The theoretical development of this variation is similar to the development demonstrated in the previous section for LSMCVSVM; therefore, we present only a small introductory development and leave further development for future work. We start first by stating the optimization problem people need to optimize, which is the minimization of the following objective function:

$$\begin{array}{c}\hfill {\displaystyle min\frac{1}{2}{\mathit{\psi }}^{\mathsf{T}}{\mathbf{\Sigma }}_w\mathit{\psi }+{\lambda }_1\sum _{i:y_i=1}^n{\xi }_i^2+{\lambda }_{-1}\sum _{i:y_i=-1}^n{\xi }_i^2}\end{array}$$

under the equality constraints:

$$\begin{array}{c}\hfill y_i({\mathit{\psi }}^{\mathsf{T}}{\mathit{x}}_i-t)=1-{\xi }_i,\end{array}$$

If someone follows the proper procedure, then the solution will be:

$$\begin{array}{c}\hfill \mathit{r}={\left({\mathbf{\Lambda }}_D{\mathbf{\Sigma }}_w^{\ast }+{\left({\mathit{X}}^{\ast }\right)}^{\mathsf{T}}{\mathit{X}}^{\ast }\right)}^{-1}{\left({\mathit{X}}^{\ast }\right)}^{\mathsf{T}}{\mathit{D}}_{\mathit{Y}}{\mathbf{1}}_n\end{array}$$

where ${\Lambda }_D$ is the diagonal matrix that has, as the $i^{\mathrm{th}}$ entry on the main diagonal, the quantity $(1/{\lambda }_1)I(Y_i=1)+(1/{\lambda }_{-1})I(Y_i=-1)$, where $I(·)$ denotes the indicator function. There are different suggestions to select the two different costs, although the most frequently used in the literature (see for example [21]) is ${\lambda }_1/{\lambda }_{-1}=n_{-1}/n_1$ where $n_i$ is the number of observations in class $i=\{-1,1\}$. Here, we emphasize that the topic of imbalance is a very rich topic, in terms of literature, with hundred methods available on how to address imbalance in the SVM framework (see [22] for a comprehensive overview). Therefore, we prefer to address this topic separately, as this will give us a way to check the impact that our new algorithm may have in addressing imbalance.

In addition to the development of new classification methodologies by utilizing the SVM algorithm and its variants, in the literature, there are more ideas which use the SVM type of algorithms to implement new approaches. One such way of utilizing the new algorithm beyond the classification framework is its use in the sufficient dimension reduction (SDR) framework. Recently, SVMs have been introduced extensively in SDR (see for example [23,24,25]), and, therefore, much more developments may be studied in SDR by utilizing new algorithms. The least squares approach, which allows for the use of analytic solution to estimate the optimal hyperplane, might lead to the use of a real-time dimension reduction method, as was demonstrated by [23].

Finally, an anonymous reviewer pointed out to us that there was a similar work that was performed much earlier than our work. The work by [26] discusses a similar idea, as indicated by the title. Unfortunately, we were not able to find a version of that paper to read and compare it to our work. The only pointer to the paper we found online was in Chinese, which was impossible for us to read. It may be an interesting exercise for someone to compare our developments with their development and see if there are any differences.