0

0

Keywords : SVM, Classification, Bernstein polynomial, Variables Selection, Non asymptotic Error Bound.

1 Introduction↩︎

The SVM, [1], is a powerful classifier used for binary classification to perform the prediction accuracy. Its motivation comes from geometric considerations and it seeks a hyperplane that separates two classes of data points by the largest margins (i. e. minimal distances of observations to a hyperplane). The classification rule depends only on a subset of observations, called support vectors, which lie along the lines indicating the width of the margin. It consists of classifying a test observation based on which side of the maximal margin hyperplane it lies. However, SVM can be less efficient in high dimensional setting with a large number of variables with only few of them are relevant for classification. Nowadays, problems with sparse scale data are common in applications such as finance, document classification, image analysis and gene expression analysis. Many sparse regularized SVM approaches are proposed to control the sparsity of the solution and to achieve both variable selection and classification. To list a few, SVM with \(\ell_1\) penalty ([2]; [3]), SVM with elastic net [4], SVM with the SCAD penalty and SVM with the combination of SCAD and \(\ell_2\) penalty [5].

The objective function of these regularized SVM approaches is not differentiable at 1, so standard optimization techniques cannot be directly applied. To overcome this problem, recent alternatives are developed based on the main idea of approximating the hinge loss function with a modified smooth function which is differentiable. Then, the use of maximization-minimization (MM) principle together with a coordinate descent algorithm can be employed efficiently to update each component of the vector parameter of the SVM model. This idea was recently implemented in the generalized coordinate descent algorithm (gCDA) by [6] and the SVM cluster correlation-network by [7]. Both authors considered the Huber loss function as a smooth approximation of the hinge loss function. They used the fact that the Huber loss is differentiable with a Lipschitz first derivative to build an upper-bound quadratic surrogate function for the Huber coordinate-wise objective function. Then, they solved the corresponding problem by a generalized coordinate descent algorithm using the MM principle to ensure the descent property. Other smooth approximations of the hinge loss can be found in the \(\ell_2\) loss linear SVM of [8] and [9]. Another algorithm that is known to be computationally efficient and might be adapted to solve penalized SVM is the iterative re-weighted least squares (IRLS-type) algorithm, which solves the regularized logistic regression [10]. However, all the objective functions of the aforementioned methods are not twice differentiable, and thus, efficient IRLS-type algorithms can not be designed to solve such approximate sparse SVM problems.

On the other side, some works on theoretical guaranties of sparse SVM (or \(\ell_1\) SVM) are recently proposed, and are essentially based on the assumption that the Hessian of the SVM theoretical (expected) hinge loss is well-defined and continuous in order to use Taylor expansion or the continuity definition. For example, [11] studied the theoretical SVM loss function to overcome the differentiability of the hinge loss. Then, some appropriate assumptions about the distribution of the data \((\boldsymbol{x}_i,y_i)\), where \(\boldsymbol{x}_i\in\mathbb{R}^p\) and \(y_i\in\{-1,1\}, i=1,...,n\), were considered to ensure the existence and the continuity of the Hessian with respect to the vector of coefficients of the SVM regression model. [11] established asymptotic properties of the coefficients in low-dimension (\(p<n\)). [12] used the same assumptions to establish the variable selection consistency in high-dimensions when the true vector of coefficients is sparse. [13] exploited the same assumptions on the theoretical hinge loss in order to develop an upper bound of the error estimation of sparse SVM in high-dimension.

Motivated by the computation problem of the SVM hinge loss and the discontinuity of the second derivative of the hessian of Huber hinge loss, we propose the BernSVM loss function, which is based on the Bernstein polynomial approximation. Its second derivative is well-defined, continuous and bounded, which simplifies the implementation of a generalized coordinate descent algorithm with strict descent property and the implementation of an IRLS-type algorithm to solve penalized SVM. By construction, the proposed BernSVM loss function has suitable properties, which allows its Hessian to satisfy a restricted strong convexity [14]. This is essential for non asymptotic theoretical guaranties developed in this work. More precisely and contrary to the works cited earlier, we consider the empirical loss function to establish an upper bound of the \(\ell_2\) norm of the estimation error of the BernSVM weighted Lasso estimator. Similar to [15], we achieve an estimation error rate of order \(\sqrt{s\log(p)/n}\).

This paper is organized as follows. In Section 2, we give more details about the construction of the BernSVM loss function, using the Bernstein polynomials. Some properties of this function are oulined. Then we describe, in the same section, the two algorithms to solve the solution path of the penalized BernSVM. The first algorithm is a generalized coordinate descent and the second one is an IRLS type algorithm. In Section 3, we undertake the theoretical properties of the penalized BernSVM with weighetd lasso penalty. In particular, we derive an upper bound of the \(\ell_2\) norm of the error estimation. Then, we extend this result to penalized BernSVM with a non-convex penalty (SCAD or MPC), using local linear approximation (LLA) algorithm. The empirical performance of BernSVM based on simulation studies and application to real datasets is detailed in Section 4. Finally, a discussion and some conclusions are given in Section 5.

2 The penalized BernSVM↩︎

In this section we consider the penalized BernSVM in high dimension for binary classification with convex and non-convex penalties. We first present the BernSVM loss function, which is a spline polynomial of degree fourth, as a smooth approximation to the hinge loss function. Then, we present two efficient algorithms for computing the solution path of a solution to the penalized BernSVM problem. The first one is based on coordinate descent scheme and the second one is an IRLS type algorithm.

Assume we observe a training data of \(n\) pairs, \(\{(y_1,{\mathbf{x}}_1),\ldots,(y_n,{\mathbf{x}}_n)\}\), where \({\mathbf{x}}_i\in I\!\!R^p\) and \(y_i\in\{-1,1\}, i=1,...,n,\) denotes class labels.

2.1 The fourth degree approximation spline↩︎

We consider the problem of smoothing the SVM loss function \(v\) : \(t\mapsto (1-t)_{+}\). The idea is to fix some \(\delta>0\), and to construct the simplest polynomial spline \[\label{loss} B_{\delta}(t)= \left\{\begin{array}{ll} v(t)\;, & \mathrm{i}\mathrm{f}\;|t-1|\;>\delta,\\ g_{\delta}(t)\;, & \mathrm{i}\mathrm{f}\;|t-1|\;\leq\delta, \end{array}\right.\tag{1}\] with \(g_{\delta}\) is a decreasing polynomial, convex, and such that \(B_{\delta}\) is twice continuously differentiable. We have in particular \(B_{\delta}(t) \rightarrow v(t)\), as \(\delta\rightarrow 0\), for all \(t \in \mathbb{R}\). These conditions can be rewritten in terms of the polynomial \(g_{\delta}\) alone as

• \(g_{\delta}(1-\delta)=\delta\), \(g_{\delta}^{'}(1-\delta)=-1\), and \(g_{\delta}^{''}(1-\delta)=0\);

• \(g_{\delta}(1+\delta)=0\), \(g_{\delta}^{'}(1+\delta)=0\), and \(g_{\delta}^{''}(1+\delta)=0\);

• \(g_{\delta}^{''}\geq 0.\)

Consider the affine transformation \(t\mapsto(t-1+\delta)/2\delta\), which maps the interval \([1-\delta,\;1+\delta]\) on \([0\), 1\(]\), and let \[q_{\delta}(x)=g_{\delta}(2\delta x+1-\delta)\;,\;x\in\;[0,\;1].\] In terms of \(q_{\delta}\), the three conditions above become

• \(q_{\delta}(0)=\delta\), \(q_{\delta}^{'}(0)=-2\delta\), and \(q_{\delta}^{''}(0)=0\);

• \(q_{\delta}(1)=0\), \(q_{\delta}^{'}(1)=0\), and \(q_{\delta}^{''}(1)=0\);

• \(q_{\delta}^{''}\geq 0.\)

These conditions can be dealt with in a simple manner in terms of a Bernstein basis.

2.1.1 The BernSVM loss function↩︎

The members of the Bernstein basis of degree \(m\), \(m\geq 0\), are the polynomials \[b_{k,m}(x)=\;\left(\begin{array}{l} m\\ k \end{array}\right)x^{k}(1-x)^{m-k},\;x\in\;[0,\;1],\] for \(0\leq k\leq m\). The coefficients of a polynomial \(P\) in the Bernstein basis of degree \(m\), will be denoted by \(\{c(k,\;m;P)\;:\;k=0,\;.\;.\;.\;,\;m\}\), and so we have \[P(x)=\sum_{k=0}^{m}c(k,\;m;P)b_{k,m}(x)\;,\;x\in\;[0,\;1].\] Let \(\triangle\) be the forward difference operator. Here, when applied to a function \(h\) of two arguments: \((k,\;m) \mapsto h(k,\;m)\), it is understood that \(\triangle\) operates on the first argument: \[\triangle h(k,\;m)=h(k+1,\;m)-h(k,\;m)\quad and\] \[\triangle^{2}h(k,\;m)=h(k+2,\;m)-2h(k+1,\;m)+h(k,\;m), \quad for all \quad (k, m).\] A basic fact related to the Bernstein basis is the following for \(0\leq k\leq m,\) we have \[b_{k,m}'=m(b_{k-1,m-1}-b_{k,m-1})=-m\triangle b_{k-1,m-1} \quad for\quad m\geq 1 \quad qnd \quad 0\leq k\leq m,\] with the convention \(b_{j,m}=0\) for \(j\not\in\{0,\;.\;.\;.\;,\;m\}.\)

A useful consequence shows how derivatives of polynomials represented in the Bernstein basis act on the coefficients. Indeed, we have

• \(c(k,\;m-1;P')=m\triangle c(k,\;m;P)\), \(0\leq k\leq m-1\), \(m\geq 1\), so that \[P^{'}(x)=m\sum_{k=0}^{m-1}\triangle c(k,\;m;P)b_{k,m-1}(x),\;x\in\;[0,\;1],\]

• \(c(k,\;m-2;P) =m(m-1)\triangle^{2}c(k,\;m;P)\), \(0\leq k\leq m-2, m\geq 2\), so that \[P^{''}(x)=m(m-1)\sum_{k=0}^{m-2}\triangle^{2}c(k,\;m;P)b_{k,m-2}(x),\;x\in\;[0,\;1].\]

Then, it is easy to see that it is hopeless to find a cubic spline \(B_{\delta}\) satisfying the constraints. If such a spline exists, its second derivative will be a polynomial spline of degree at most one, which must equal zero at the endpoints \(1-\delta\) and \(1+\delta\). By continuity of the second derivative, it therefore equals \(0\) everywhere. The first derivative is therefore a degree zero polynomial on the entire axis and must equal \(-1\) by the first condition. This contradicts the second condition.

The proof is differed to Appendix A.

The following proposition summarizes the properties of the first and the second derivative of the BernSVM loss function \(B_{\delta}(.)\) defined in equation (1).

2.1.2 Graphical comparaison of the three loss functions↩︎

In this section, we examine some properties of the BernSVM loss function through a graphical illustration. We also compare its behavior with the standard SVM loss function \(v(t)\) and the hinge loss function considered in HHSVM [6], which is given by \[\phi_c(t) = \left\{ \begin{array}{ll} 0 & ,t \geq 1, \\ \frac{(1-t)^2}{2\delta} & ,1 -\delta <t\leq 1, \\ 1-t-\frac{\delta}{2} & ,t \leq 1-\delta. \\ \end{array} \right.\]

Figure 1 shows this graphical illustration. When \(\delta\) is small (i.e. close to zero), the three loss functions are almost identical (Panel (a)). When \(\delta \geq 1\) (Panel (b)), both HHSVM and BernSVM loss functions have the same shape that the SVM loss. However, BernSVM approximates better the SVM loss function. Moreover, BernSVM loss function is twice differentiable everywhere.

2.2 Algorithms for solving penalized BernSVM↩︎

In this section, we propose the penalized BernSVM as an alternative to the penalized SVM and the penalized HHSVM for binary classification in high dimension settings. We are assuming \(n\) pairs of training data \((\boldsymbol{x}_i,y_i)\), for \(i=1,...,n\), with \(\boldsymbol{x}_i\in \mathbb{R}^p\) predictors and \(y_i \in\{-1,1\}\). We assume that the predictors are standardized: \[\frac{1}{n}\sum_{i=1}^{n}x_{ij} = 0, \quad \frac{1}{n}\sum_{i=1}^{n}x_{ij}^2 = 1.\]

We propose to solve the penalized BernSVM problem given by \[\label{problem} \underset{(\beta_0,\boldsymbol{\beta})}{\min} \frac{1}{n} \sum_{i=1}^{n}B_{\delta}(y_i(\beta_0 + \boldsymbol{x}_i^{\top}\boldsymbol{\beta})) + P_{\lambda_1,\lambda_2}(\boldsymbol{\beta}),\tag{4}\] where the objective function \(B_{\delta}(.)\) is defined in (1 ) and \[P_{\lambda_1,\lambda_2}(\boldsymbol{\beta}) = \sum_{j=1}^{p} P_{\lambda_1}(|\beta_j|) + \frac{\lambda_2}{2} ||\boldsymbol{\beta}||_2^2\] is a penalty function added for obtaining sparse coefficients’ estimates. The penalty \(P_{\lambda_1}(|\beta_j|)\) takes different forms:

• If \(P_{\lambda_1}(|\beta_j|) = \lambda_1 |\beta_j|\), \(P_{\lambda_1,\lambda_2}(\boldsymbol{\beta})\) is the elastic net penalty (EN);

• If \(P_{\lambda_1}(|\beta_j|) = \lambda_1w_j|\beta_j|\), where the weight \(w_j > 0\) for all \(j\), then \(P_{\lambda_1,\lambda_2}(\boldsymbol{\beta})\) is the adaptive Lasso + L2 norm penalty (AEN);

• If one set \[P_{\lambda_1}(|\beta_j|) = \lambda_1 |\beta_j| \mathbb{1}_{|\beta_j|\leq \lambda_1} - \frac{|\beta_j|^2 - 2a\lambda_1|\beta_j| + \lambda_1^2}{2(a-1)}\mathbb{1}_{\lambda_1 < |\beta_j|\leq a\lambda_1} + \frac{(a+1)\lambda_1^2}{2}\mathbb{1}_{|\beta_j| > a \lambda_1},\] with \(a > 2\), then \(P_{\lambda_1,\lambda_2}(\boldsymbol{\beta})\) is the SCAD + L2 norm penalty (SCADEN);

• If one set \[P_{\lambda_1}(|\beta_j|)=\lambda_1 (|\beta_j| - \frac{|\beta_j|^2}{2\lambda_1 a}) \mathbb{1}_{|\beta_j|< \lambda_1 a} + \frac{\lambda_1^2 a}{2}\mathbb{1}_{|\beta_j| \geq \lambda_1 a},\] where \(a > 1,\) then in this case, \(P_{\lambda_1,\lambda_2}(\boldsymbol{\beta})\) is the MCP + L2 norm penalty (MCPEN).

Before going through details of the proposed algorithms to solve BernSVM optimization problem, we would like to emphasize that, for small \(\delta\), the minimizer of the penalized BernSVM in (4 ) is a good approximation to the minimizer of the penalized SVM. This is outlined in the following proposition.

In the remaining of section 2, we propose two competitor algorithms for solving (4 ) with convex penalties. The first one is a generalized coordinate descent algorithm based on MM principle, while the second one is an IRLS-type algorithm. Moreover, we briefly describe a third algorithm based on local linear approximation for solving (4 ) with non-convex penalties.

2.2.1 The coordinate descent algorithm for penalized BernSVM↩︎

In this section, we present a generalized coordinate descent (GCD) algorithm to solve penalized BernSVM, called BSVM-GCD, similar to the GCD approach proposed by [6]. Because the second derivative of the BernSVM objective function is bounded, we approximate the coordinate wise objective function by a surrogate function and we use the coordinate descent to solve the optimization problem in (4 ). The proposed approach is described next.

Notice, first, that the coordinate objective function of problem (4 ) can be written as follows \[\label{3} F(\beta_j|\beta_0,\tilde{\beta}_j) := \frac{1}{n} \sum_{i=1}^{n}B_{\delta}\{r_i + y_{i}x_{ij}(\beta_j - \tilde{\beta}_j)\} + P_{\lambda_1,\lambda_2}(\beta_j),\tag{5}\] where \(r_i = y_i (\tilde{\beta}_0 + \boldsymbol{x}_i^{\top}\tilde{\boldsymbol{\beta}})\), and \(\tilde{\beta}_0\) and \(\tilde{\boldsymbol{\beta}}\) are the current estimates of \(\beta_0\) and \(\boldsymbol{\beta}\), respectively. As the loss function, \(B_{\delta}(.)\), is convex, twice differentiable and has a second derivative bounded by \(L = 3/4\delta\), one can rely on the mean value theorem and approximate the objective function in (5 ) by a (surrogate) quadratic function, \(Q_{\delta}(.)\), and then use the MM principle to solve the new problem for each \(\beta_j\) holding \(\tilde{\beta}_0\) and \(\beta_k=\tilde{\beta}_j\), for \(k\ne j\), fixed at the current iteration \[\label{QFunct} \underset{\beta_j)}{\min} Q_{\delta}(\beta_j|\tilde{\beta_0},\tilde{\beta}_j),\tag{6}\] where \[\label{obj} Q_{\delta}(\beta_j|\beta_0,\tilde{\beta}_j) = \frac{\sum_{i=1}^{n}B_{\delta}(r_i)}{n} + \frac{\sum_{i=1}^{n}B_{\delta}^{'}(r_i)y_ix_{ij}}{n}(\beta_j - \tilde{\beta}_j) + \frac{L}{2}(\beta_j - \tilde{\beta}_j)^2 + P_{\lambda_1,\lambda_2}(\beta_j).\tag{7}\] The next proposition gives the explicit solution of (6 ) for the different forms of \(P_{\lambda_1}(\beta_j)\).

The proof is outlined in Appendix A.

Algorithm 2 (BSVM-GCD), given next, gives details of steps for solving the objective function in (7 ) using both coordinate descent and MM principle.

The upper bound of the second derivative of the BernSVM loss function can be reached in \(t = 1\). This means that \(F(t|\beta_0,\tilde{\beta}_j) \leq Q(t|\beta_0,\tilde{\beta}_j)\). To reach a strict inequality, one can relax the upper bound \(L\) and use \(\tilde{L} = (1 + \epsilon)L\), where \(\epsilon = 1e^{-6}\). Hence, Algorithm 2 is implemented using \(\tilde{L}\) in the place of \(L\). This leads to the strict descent property of the BSVM-GCD algorithm, which is given in the next proposition.

This proposition shows that the objective function \(F(.)\) decreases after each majorization update, which means that for \(\tilde{\beta}^{new}_j \neq \tilde{\beta}_j\) we have \[F(\beta^{new}_j|\beta_0,\tilde{\beta}_j) < F(\tilde{\beta}_j|\beta_0,\tilde{\beta}_j).\] This result shows that the BSVM-GCD algorithm enjoys strict descent property.

2.2.2 The IRLS algorithm for penalized BernSVM↩︎

In this section, we present an IRLS-type algorithm combined with coordinate descent, termed BSVM-IRLS, to solve the penalized BernSVM in (4 ). In their approach, [10] used IRLS to solve the regularized logistic regression because the second derivative exist (glmnet R package). On the other hand, [6] used a GCD algorithm to solve the regularized HHSVM because the Huber loss does not have the second derivative everywhere (gcdnet R packge). Our BernSVM loss function, in contrast, is very smooth and thus one can rely on this advantage and solve (4 ) using IRLS trick combined with coordinate descent. The BSVM-IRLS approach is described below.

Firstly, we consider the non-penalized BernSVM with loss function \[\label{2} \ell{(\beta_0,\boldsymbol{\beta})} = \frac{1}{n} \sum_{i=1}^{n}B_{\delta}(y_i(\beta_0 + \boldsymbol{x}_i^{\top}\boldsymbol{\beta})) = \frac{1}{n} \sum_{i=1}^{n}B_{\delta}(y_i\boldsymbol{X}_i^{\top}\boldsymbol{b}),\tag{8}\] where, \(\boldsymbol{b} = (\beta_0,\boldsymbol{\beta}^{\top})^{\top}\) and \(\boldsymbol{X} = (1,\boldsymbol{x}^{\top})^{\top}\). To solve \[\frac{\partial{\ell}}{\partial \boldsymbol{b}}(\beta_0,\boldsymbol{\beta}) = \frac{1}{n} \sum_{i=1}^{n}y_i\boldsymbol{X}_iB^{'}_{\delta}(y_i\boldsymbol{X}_i^{\top}\boldsymbol{b}) = 0,\] one can use the Newton-Raphson algorithm because the loss function has a continuous second derivative defined by \[\frac{\partial^{2}{\ell}}{\partial \boldsymbol{b}\partial \boldsymbol{b}^{\top}}(\beta_0,\boldsymbol{\beta}) = \frac{1}{n} \sum_{i=1}^{n}\boldsymbol{X}_iB^{''}_{\delta}(y_i\boldsymbol{X}_i^{\top}\boldsymbol{b})\boldsymbol{X}_i^{\top}.\] The update of the coefficients is given then by \[\boldsymbol{b}^{new} = \tilde{\boldsymbol{b}} - (\frac{\partial^{2}{\ell}}{\partial \boldsymbol{b}\partial \boldsymbol{b}^{\top}}(\beta_0,\boldsymbol{\beta}))^{-1} \frac{\partial{\ell}}{\partial \boldsymbol{b}}(\beta_0,\boldsymbol{\beta})\tilde{\boldsymbol{b}}.\] Let \(u_i = y_i.B^{'}_{\delta}(y_i\boldsymbol{X}_i^{\top}\tilde{\boldsymbol{b}})\), \(\phi_{ii} = B^{''}_{\delta}(y_i\boldsymbol{X}_i^{\top}\tilde{\boldsymbol{b}})\) and \(\tilde{\boldsymbol{\Phi}} = diag(\phi_{ii})\) for \(i=1\ldots, n\). Then, the update of the Newton-Raphson algorithm can be written as follows, \[\begin{align} \boldsymbol{b}^{new}& = &\tilde{\boldsymbol{b}} - (\boldsymbol{X}^{\top}\boldsymbol{\Phi}\boldsymbol{X})^{-1}\boldsymbol{X}^{\top}\boldsymbol{u} \\ &= &(\boldsymbol{X}^{\top}\boldsymbol{\Phi}\boldsymbol{X})^{-1}\boldsymbol{X}^{\top}\boldsymbol{\Phi}[\boldsymbol{X}\tilde{\boldsymbol{b}} - \boldsymbol{\Phi}^{-1}\boldsymbol{u}]\\ &= &(\boldsymbol{X}^{\top}\boldsymbol{\Phi}\boldsymbol{X})^{-1}\boldsymbol{X}^{\top}\boldsymbol{\Phi}\boldsymbol{z}, \end{align}\] where \(\boldsymbol{z} = \boldsymbol{X}\tilde{\boldsymbol{b}} - \boldsymbol{\Phi}^{-1}\boldsymbol{u}\). Thus, \(\boldsymbol{b}^{new}\) is the update of a weighted least regression problem with response \(\boldsymbol{z}\), with the loss function \[\ell{(\beta_0,\boldsymbol{\beta})} \approx \frac{1}{2n} \sum_{i=1}^{n} \phi_{ii}(z_i - \boldsymbol{X}_i^{\top}\boldsymbol{b})^2.\] This problem can be solved by an IRLS algorithm in the case \(p < n\). In high dimensional settings, the penalized BernSVM, defined in (4 ), can be solved by combining IRLS with a cyclic coordinate descent. Thus, for the Adaptive Elastic Net penalty, the solution is given by \[\label{EN} \hat{\beta}^{AEN}_j = \frac{S(\frac{\sum_{i=1}^{n}\phi_{ii}x_{ij} r_i}{n} + \frac{\sum_{i=1}^{n}\phi_{ii}x_{ij}^2 }{n}\tilde{\beta_j},w_j\lambda_1)}{\lambda_2 + \frac{\sum_{i=1}^{n}\phi_{ii}x_{ij}^2 }{n}} \quad \text{and} \quad \hat{\beta}^{AEN}_0 = \frac{\sum_{i=1}^{n}\phi_{ii}r_i}{\sum_{i=1}^{n}\phi_{ii}},\tag{9}\] where \(r_i = z_i - \tilde{\beta}_0 - \boldsymbol{x}_i^{\top}\tilde{\boldsymbol{\beta}}\) and \((\tilde{\beta}_0,\tilde{\boldsymbol{\beta}})\) is the estimate of \(({\beta}_0,\boldsymbol{\beta})\) obtained on the current iteration. Of note, to deal with the zero weights (i.e, \(\phi_{ii}=0\) for some \(i\)), one can replace them by a constant \(\xi = 0.001\) or use a modified Neweton-Raphson (NR) algorithm and replace all the weights by an upper bound of the second derivative. In our work, we adopted the modified NR algorithm and set \(\phi_{ii} = L\) for \(i=1,\ldots,n\), where \(L = \frac{3}{4\delta}\) is the upper bound of the second derivative of the proposed loss function \(B_{\delta}(.)\) given by (1 ).

2.2.3 Local linear approximation algorithm for BernSVM with non-convex penalties↩︎

We consider in this section the penalized BernSVM with non-convex penalty SCAD or MCP and we propose to solve the resulting optimization problem using the local linear approximation (LLA) of [16].

Without loss of generality, we assume in this section that \(\lambda_2 = 0\). The LLA approach is an appropriate approximation of SCAD or MCP penalty. It is based on the first order Taylor expansion of the MCP or SCAD penalty functions around \(|\tilde{\beta}_j|\), and adapts both penalties to be solved as an iterative adaptive Lasso penalty; the weights are updated/estimated at each iteration. In our context, this means that the penalized BernSVM with SCAD or MCP penalty can be solved iteratively as follows \[\label{weighted1} (\tilde{\beta}_0,\tilde{\boldsymbol{\beta}})^{new} = \arg\min_{(\beta_0,\boldsymbol{\beta})} \frac{1}{n} \sum_{i=1}^{n}B_{\delta}(y_i(\beta_0+\boldsymbol{x}_i^{\top}\boldsymbol{\beta})) + \lambda_1 ||\hat{\boldsymbol{W}}\boldsymbol{\beta}||_1,\tag{10}\] where \(\hat{\boldsymbol{W}}= diag\{\hat{w}_j\}\), with \(\hat{w}_j = P^{'}_{\lambda_1}(\tilde{\beta_j})/\lambda_1, j=1,\ldots,p,\) are the current weights, and \(P^{'}_{\lambda_1}(.)\) is the first derivative of the SCAD or MCP penalty given respectively by \[P^{'}_{\lambda_1}(t) = \lambda_1 \mathbb{1}_{\{t \leq \lambda_1\}} + \frac{(a\lambda_1 - t)_+}{a - 1}\mathbb{1}_{\{t > \lambda_1\}},\] and \[P^{'}_{\lambda_1}(t) = (\lambda_1 - \frac{t}{a})_+.\] To sum up, penalized BernSVM with LLA penalty is an iterative algorithm, whcih solves at each iteration BSVM-IRLS (or BSVM-GCD). The steps of the LLA algorithm for penalized BernSVM are described in Algorithm 4 below.

3 Theoretical guaranties for penalized BernSVM estimator↩︎

In the next section, we develop an upper bound of the \(\ell_2\) norm of the estimation error \(||\hat{\boldsymbol{\beta}} - \boldsymbol{\beta}^{\star}||_2\) for penalized BernSVM with weighted lasso (\(\lambda_2 = 0\)) and SCAD or MCP penalties. Here, \(\hat{\boldsymbol{\beta}}\) is the minimizer of the corresponding penalized BernSVM and \(\boldsymbol{\beta}^{\star}\) is the minimizer of the population version of (4 ) without the penalty term; i.e., \(\boldsymbol{\beta}^{\star}\) is the minimizer of \(\mathbb{E}[B_{\delta}(y\boldsymbol{x}^\top\boldsymbol{\beta})]\). For seek of simplicity, the intercept is omitted for the theoretical development.

3.1 An upper bound of the estimation error for weighted lasso penalty↩︎

In this section, we derive an upper bound of the estimation error \(||\hat{\boldsymbol{\beta}} - \boldsymbol{\beta}^{\star}||_2\) for penalized BernSVM with weighted lasso penalty using the properties of the BernSVM loss function.

The weighted penalized BernSVM optimization problem is defined as follows \[\label{weighted} \arg\min_{\boldsymbol{\beta}\in \mathbb{R}^p}L(\boldsymbol{\beta}): = \frac{1}{n} \sum_{i=1}^{n}B_{\delta}(y_i\boldsymbol{x}_i^{\top}\boldsymbol{\beta}) + \lambda_1 ||\boldsymbol{W}\boldsymbol{\beta}||_1,\tag{11}\] where \(\boldsymbol{W} = diag(\boldsymbol{w})\) is a diagonal matrix of known weights, with \(w_j \geq 0\), for \(j=1\ldots p\). The (unweighted) Lasso is a special case of (11 ), with \(w_j = 1, j\in \{1,\ldots,p\}\).
Let \(\hat{\boldsymbol{\beta}}\) be a minimizer of (11 ), and assume that the population-level minimizer \(\boldsymbol{\beta}^{\star}\) defined above is sparse and unique. Define \(S\subset \{1,2,...,p\}\) the set of non-zero elements of \(\boldsymbol{\beta}^{\star}\), with \(|S| = s\), and \(S^{c}\) is the complement of \(S\).
In order to obtain an upper bound of the estimation error, the BernSVM loss function needs to satisfy the strongly convex assumption. This means that the minimum of the eigenvalues of its corresponding Hessian matrix, defined as \(\boldsymbol{H} = \boldsymbol{X}^{\top}\boldsymbol{D}\boldsymbol{X}/n\), is lower bounded by a constant \(\kappa > 0\); the matrix \(\boldsymbol{D} = diag\{B_{\delta}^{"}(y_i\boldsymbol{x}_i^{\top}\boldsymbol{\beta}^{\star})\}\), \(i=1,\ldots,n\). In high-dimensional settings, with \(p > n\), \(\boldsymbol{H}\) has rank lower or equal to \(n\). Therefore, we can not reach the strong convexity assumption of the Hessian matrix. To overcome this problem, one can look for a restricted strong convexity assumption to be verified on a restricted set \(\mathcal{C}\subset \mathbb{R}^p\) [14], which is defined next.

We also state the following assumptions:

• Let \[\mathcal{C}(S) = \{\boldsymbol{h}\in \mathbb{R}^{p}:||\boldsymbol{W}_{\boldsymbol{S}^{c}}\boldsymbol{h}_{\boldsymbol{S}^{c}}||_1 \leq \gamma ||\boldsymbol{h}_S||_1 \},\] where \(\gamma > 0\) and \(\boldsymbol{h}_{\mathcal{A}}=(h_j:j\in \mathcal{A})\), with \(\mathcal{A}\subset \{1,\ldots,p\}\), is a sub-vector of \(\boldsymbol{h}\). We assume that the design matrix, for all \(\boldsymbol{h} \in \mathcal{C}\), satisfies the condition \[\frac{||\boldsymbol{X}\boldsymbol{h}||_2^2}{n} \geq \kappa_1 ||\boldsymbol{h}||_2^2 - \kappa_2 \frac{\log(p)}{n}||\boldsymbol{h}||^2_1,\] with probability \(1- \exp(-c_0n)\), for some \(c_0>0\), where \(\boldsymbol{x}_i, i= 1,\ldots, n,\) are i.i.d Gaussian or Sub-Gaussian.

• There exist a ball, \(\boldsymbol{B}(\boldsymbol{x}_0,r_0)\), centered at a point \(\boldsymbol{x}_0\), with radius \(r_0\), such as \(B^{''}_{\delta}\circ f(\boldsymbol{x}_0) > 0\), with \(f(\boldsymbol{x}_0) = y \boldsymbol{x}^{\top}_0 \boldsymbol{\beta}^{\star}\), and for any \(\boldsymbol{x} \in \boldsymbol{B}(\boldsymbol{x}_0,r_0)\), we have \(B^{''}_{\delta}\circ f(\boldsymbol{x}) > \kappa_3\), for some constant \(\kappa_3> 0\).

• The columns of the design matrix \(\boldsymbol{X}\) are bounded: \[\forall j\in\{1,...,p\} \quad ||\boldsymbol{x}_{j}||_2 \leq M.\]

• The densities of \(\boldsymbol{X}\) given \(\boldsymbol{y}=1\) and \(\boldsymbol{y}=-1\) are continuous and have at least the first moment.

[17] show that Assumption A1 is satisfied with high probability for Gaussian random matrices. [18] extend this result for Sub-Gaussian random matrices. The existence of \(\boldsymbol{x}_0\) in Assumption A2 is given by the behavior of the second derivative \(B^{''}_{\delta}(.)\) (see Proposition 1). In fact we have \(B^{''}_{\delta}(1) = \frac{3}{4\delta}> 0\), then if we let \(f(\boldsymbol{x}_0) = 1\), we can choose \(\boldsymbol{x}_0 \in f^{-1}(\{1\}) \subset \mathbb{R}^{p}\). Hence by the continuity of \(B^{''}_{\delta}\circ f(.)\) we can choose \(\kappa_3 = \frac{3}{4\delta}-\epsilon>0\) for some small \(\epsilon>0\). Note that this assumption is similar to Assumptions A2 and A4 in [11], to ensure the positive-definiteness of the Hessian around \(\boldsymbol{\beta}^{\star}\). Assumption A3 is needed to ensure that the first derivative of the loss function is a bounded random variable around \(\boldsymbol{\beta}^{\star}\). Moreover, Assumption A4 is needed to ensure the existence of the first derivative of the population version of the loss function. Assumption A4 is similar to A1 in [11]. However, in Assumption A2 of [11], one needs at least two moments to ensure the continuity of the Hessian around \(\boldsymbol{\beta}^{\star}\).

The next lemma establishes that the hessian matrix \(\boldsymbol{H}\) satisfies the RSC condition with high probability.

The proof of Lemma 1 is given in Appendix B.

The following theorem establishes an upper bound of the penalized BernSVM with Adaptive Lasso.

The proof of Theorem 2 is given in Appendix B. In Theorem 2, the upper bound depends on \(\lambda_1\). Therefore, a choice of \(\lambda_1\) is necessary to capture the near-oracle property of the estimator [15]. Note that in Theorem 2, we have used the fact that the first derivative of the loss function is bounded. In the next lemma, we use a probabilistic upper bound of the first derivative of the BernSVM loss function with aim to give a best choice of \(\lambda_1\) and to guarantee a rate of \(\mathbf{O}_{P}(\sqrt{s\log(p)/n})\) for the error \(||\hat{\boldsymbol{\beta}} - \boldsymbol{\beta}^{\star}||_2\).

The proof of this lemma is given in Appendix C. The following theorem outlines the rate of convergence of the error \(||\hat{\boldsymbol{\beta}} - \boldsymbol{\beta}^{\star}||_2\).

The proof of this theorem is similar to the proof of Theorem 2. It is detailed in Appendix C.

In the next section we will derive an upper bound for BernSVM with weighted Lasso, where the weights are estimated.

3.2 An upper bound of the estimation error for weighted Lasso with estimated weights↩︎

In this section, let \(w_j\) be unknown positive weights and \(\hat{w}_j\) their corresponding estimates, for \(j=1,\ldots,p\), which are assumed to be non-negative. Let \(\hat{\boldsymbol{\beta}}\) be a minimizer of the weighted Lasso with weights \(\hat{w}_j\), defined as follows \[\hat{\boldsymbol{\beta}} = \arg\min_{\boldsymbol{\beta}\in \mathbb{R}^p} \frac{1}{n} \sum_{i=1}^{n}B_{\delta}(y_i\boldsymbol{x}_i^{\top}\boldsymbol{\beta}) + \lambda_1 \sum_{i=1}^{p}\hat{w}_j|\beta_j|.\] To derive an upper bound of \(||\hat{\boldsymbol{\beta}} - \boldsymbol{\beta}^{\star}||_2\), we suppose that the following event \(\Omega_0\) is satisfied [19] \[\Omega_0 = \{\hat{w}_j \leq w_j, \quad \forall j \in \boldsymbol{S}\} \cap \{ w_j \leq \hat{w}_j,\quad \forall j \in \boldsymbol{S}^{c}\}.\] The next theorem gives an upper bound of the BernSVM with weighted Lasso, where the weights are estimated.

The proof of this theorem is outlined in Appendix D.

In the next section, we extend Theorem 3 for non-convex penalties.

3.3 An extension of the upper bound for non-convex penalties↩︎

We study here the penalized BernSVM estimator with the non-convex penalties SCAD and MCP. We establish an upper bound of its error estimation under some regularity conditions.
Recall that Algorithm 3 describes steps to compute the solution of BernSVM with the non-convex penalties, as defined in Equation (10 ). The weights \(\hat{w}_j, j= 1,\ldots, p,\) are computed iteratively using the first derivative of SCAD or MCP penalty. To tackle such a problem, we search first for an upper bound of \(||\hat{\boldsymbol{\beta}}- \boldsymbol{\beta}^{\star} ||_2\), where \(\hat{\boldsymbol{\beta}}\) is an estimator of Equation (10 ). The next theorem states this result.

The proof of Theorem 5 is given in Appendix E. We are now able to derive an upper bound of \(||\hat{\boldsymbol{\beta}}^{(l)}- \boldsymbol{\beta}^{\star} ||_2\), where \(\hat{\boldsymbol{\beta}}^{(l)}\) is the estimator of the \(l-\)th iteration of Algorithm 3.

The proof of Theorem 6 is straightforward and can be outlined as follows. Since \(P^{'}_{\lambda_1}(t)\) is concave in \(t\), we have \(||P^{'}_{\lambda_1}(|\boldsymbol{\beta}^{\star}_S|)||_2 \leq \sqrt{s}\lambda_1\), which means that \(||P^{'}_{\lambda_1}(|\boldsymbol{\beta}^{\star}_S|)||_2 = \boldsymbol{O}(\sqrt{s\log(p)/n}).\) Then, by induction using Equation ( ?? ), for each iteration, we have \(||\hat{\boldsymbol{\beta}}^{(l)}- \boldsymbol{\beta}^{\star} ||_2 = \boldsymbol{O}_{P}(\sqrt{s\log(p)/n}).\)

4 Simulation and empirical studies↩︎

We study the empirical behavior of BernSVM with its competitors in terms of computational time, and we evaluate the methods estimation accuracy through different measures of performance via simulations and application to real genetic data sets.

4.1 Simulation study↩︎

In this section, we first compare the computation time of three algorithms: BSVM-GCD, BSVM-IRLS and HHSVM. In the second step, we examine the finite sample performance of eight regularized SVM methods with different penalties in terms of five measures of performance.

4.1.1 Simulation study designs↩︎

The data sets are generated following three scenarios described bellow.

1. Scenario 1: In this scenario the columns of \(\boldsymbol{X} \in \mathbb{R}^{n\times p}\), with \(n=100\) and \(p=5000\), are generated from a multivariate normal with mean 0 and correlation matrix \(\boldsymbol{\Sigma}\) with compound symmetry structure, where the correlation between all covariates is fixed to \(\rho = 0.5\) in this case. The coefficients are set to be as follows \[\beta^{\star}_j = (-1)^j\times \exp(-(2\times j - 1)/20)\] for \(j = 1,2,...,50\) and \(\beta^{\star}_j = 0\) for \(j = 51,...,p\). Firstly, we generated \[z_i = \boldsymbol{x}_i^{\top}\boldsymbol{\beta}^{\star} + \epsilon_i, \quad i=1,\ldots,n,\] where \(\epsilon_i \sim N(0,\sigma^2)\). The variance \(\sigma^2\) is chosen such as a signal to noise ratio (\(SNR\)) is equal to \(3\), where \(SNR =\boldsymbol{\beta}^{{\star}^{\top}}\boldsymbol{\Sigma}\boldsymbol{\beta}^{\star}/\sigma^2\). Secondly, the response variable \(y_i\) is obtained by passing \(z_i\) through an inverse-logistic to obtain the probability \(pr_i = P(y_i=1|\boldsymbol{x}_i)= 1/(1+exp(-z_i))\). Then \(y_i\) is generated from a binomial distribution with probability equal to \(pr_i\). By choosing \(\delta=2, 1, 0.5 , 0.1, 0.01\), we compare the computation time of the three algorithms mentioned above using lasso penalty (\(\lambda_2=0\)). The whole process is repeated 100 times and the average run time is reported in Table 1.

2. Scenario 2: In this scenario, we aimed to study the impact of the correlation magnitude on the run time. As in scenario 1, we generated a data set \((\boldsymbol{x}_i,y_i), i=1,\ldots,n\) with \(n=100\) and \(\boldsymbol{x}_i \in \mathbb{R}^{ p}\) with \(p=1000\). We consider different values for the correlation coefficient \(\rho\). The coefficients are generated as follows \(\beta^{\star}_j \sim U(0.9,1.1)\) for \(1\leq j \leq 25\), \(\beta^{\star}_j = -1\) for \(51\leq j \leq 75\) and \(\beta^{\star}_j = 0\) otherwise. Here, \(s = 50\) is the number of significant coefficients. By choosing in this scenario \(\delta=0.01, 0.5 , 2\) and \(\rho = 0.2,0.5,0.75,0.95\), we compare the computation time of three algorithms using lasso penalty (\(\lambda_2=0\)). The whole process is repeated 100 times and The average run time is reported in Table 2.

3. Scenario 3: This scenario was suggested by [20]. Here, it is considered to compare the performance of sparse classification methods (BernSVM, sparse logistic, HHSVM and SparseSVM) with lasso and elastic net penalties. We generate our data from a logistic model \[z_i = \log(\frac{p_i}{1-p_i}) = \beta_0 + \boldsymbol{x}^{\top}_{S,i}\boldsymbol{\beta}_{S}, \quad i=1....,n,\] where \(\boldsymbol{X} \in \mathbb{R}^{n\times p}\) with \(n = 50\) and \(p = 800\). \(S\) is the set of active variables. The columns of \(\boldsymbol{X}\) are generated from a multivariate normal with mean \(0\) and correlation matrix \(\boldsymbol{\Sigma}\). We take \(\Sigma_{ij} = \rho\) for \(i \neq j\) and \(\Sigma_{ii} = 1, \quad 1\leq i,j \leq p\), which means that all the variables are equally correlated with \(\rho = 0.2,0.5,0.8\), so the correlation matrix is different in each scenario. The number of the active predictors is taken equal to \(s = [p\xi]\) with \(\xi = 0.05, 0.3\). The coefficients of the active variables are generated as \((-1)^vu\) where \(v\sim Bernoulli(0.3)\) and \(u\) is uniformly distributed on \((0,0.5)\). The response variable \(y_i\), \(i=1,\ldots,n\), where \(y_i \in \{-1,1\}\), is generated by passing \(z_i\) through an inverse-logistic. The value of intercept \(\beta_0\) is fixed such as the conditional probability \(P(y_i=1|\boldsymbol{x}_i) = 0.3\).

   In this scenario, for each Monte Carlo replication, we simulate two data sets: 1) a training data of \(n = 50\) observations, which is used by the different methods to perform a 10-fold cross-validation in order to choose the best model, for each method; 2) a test data of \(n_{test} = 200\) observations, which is used to evaluate the methods performance. The latter is based on five measures, which are given as follows:

   • The misclassification rate, \(MR\) is defined by \[MR = \frac{\sum_{i=1}^{n_{test}} [({y_{test}}_{i}- \hat{y}_i)/2]^2}{n_{test}};\]

   • The sensitivity \(SE\) is defined as \[SE = 1 - \frac{\sum_{i\in A^{1}} [({y_{test}}_{i}- \hat{y}_i)/2]^2}{n_{test}},\] where \(A^1=\{i: {y_{test}}_{i} =1\}\);

   • The specificity \(SP\) is defined as \[SP = 1 - \frac{\sum_{i\in A^{-1}} [({y_{test}}_i- \hat{y}_i)/2]^2}{n_{test}},\] where \(A^{-1}=\{i: {y_{test}}_{i} = -1\}\);

   • The precision \(PR\) is defined by \[PR = \frac{ \neq \{j: {\beta}^{\star}_j \neq 0, \hat{\beta}_j \neq 0\}}{\neq \{\hat{\beta}_j \neq 0\}},\]

   • The recall \(RC\) is defined by \[RC = \frac{\neq \{j: {\beta}^{\star}_j \neq 0, \hat{\beta}_j \neq 0\}}{\neq \{{\beta}^{\star}_j \neq 0\}},\] where \(\hat{y}_i\) is the \(i-\)th predicted for the response variable, \(\beta^{\star}_j\) is \(j-\)th coordinate of the original coefficients and \(\hat{\beta}_j\) is the \(j-\)th coordinate of the estimated coefficients. We note that the performance of a given method in terms of SE (or SP or PR or RC) is indicated by its largest value.

The eight sparse classification methods to be compared in scenario 3 are as follows:

• BernSVM-Lasso: the BSVM-IRLS with Lasso;

• BernSVM-EN: the BSVM-IRLS with Elastic Net;

• Logistic-Lasso: the binomial model with Lasso computed using glmnet R package [10];

• Logistic-EN: the binomial model with Elastic Net computed using glmnet R package [10];

• HHSVM\(_L\) : HHSVM with Lasso computed using gcdnet R package [6];

• HHSVM\(_{EN}\): HHSVM with Elastic Net computed using gcdnet R package [6];

• SparseSVM-Lasso: the sparse SVM with Lasso using the sparseSVM R package [21];

• SparseSVM-EN: the sparse SVM with Elastic Net using the sparseSVM R package [21].

Of note, we have used two algorithms to solve the penalized BernSVM optimization problem: the BSVM-GCD algorithm based on an MM combined with coordinate descent and the BSVM-IRLS algorithm based on an IRLS scheme combined with coordinate descent.

4.1.2 Simulation results↩︎

Summing up all the scenarios considered here, we have the following observations:

\(\bullet\) Table 1 summarizes the average run time of the HHSVM and the penalized BernSVM with lasso penalty in scenario 1. Figure \(1\) provides averaged computation time curve of 100 runs of the three algorithms as a function of \(\delta\) values in the same scenario. Our proposed BSVM-IRLS approach obtains the best computation time followed by HHSVM for any \(\delta\) value. While BSVM-GCD produces the highest computation time for each delta value. These results are well summarized in Figure 1 by the gaits of the algorithms’ computation time curves, where the curve for HHSVM is intermediate between the other two.

\(\bullet\) Table 2 summarizes the average run time of the HHSVM and Bernstein penalized SVM with lasso penalty in scenario 2 with different values of \(\delta\) and \(\rho\). The result is similar to scenario 1 where BSVM-IRLS provides the best run time followed by HHSVM for any value of \(\delta\) and \(\rho\). It is interesting to notice here that the three algorithms produce their best time for the same \(\delta=0.5\) for each value of \(\rho\).

\(\bullet\) Table \(3\) report the average of the five performance measures of our methods and the competitor’s on Scenario 3, for \(\xi = 0.05,0.3\) and \(\rho = 0.2,0.5,0.8\). The results obtained for scenario 3 can be summarized as follows:

• For a fixed number of active variables, increasing rho leads to a decrease of about \(5-14\%\) in MR, a growth of about \(21-40\%\) in \(SE\), and a slight increase in \(RC\) for methods using the EN penalty(except SparseSVM-EN). While all methods are relatively stable in terms of \(SP\) and \(PR\) to any increase in \(\rho\) and/or in \(\xi\).

• On the other hand, for any fixed value of \(\rho\) and for any increase in \(\xi\), \(MR\) decreases by about \(8-16\%\), \(SE\) decreases by about \(8-40\%\) and \(RC\) increases slightly for only the methods based on the EN penalty and increases greatly of about \(11-28\%\) for BernSVM-EN and HHSVM\(_{EN}\).

• We observe that the methods are relatively comparable in terms of \(MR\) with a slight advantage to those using the EN penalty for any pair \((\rho,\xi)\). While for any fixed \(\rho\), methods based on standard penalized logistic and SVM provide good performance than the others in terms of \(SE\) for any value of \(\xi\). The difference with methods based on the smooth approximation of the hinge loss function decreases for increasing the value of \(\rho\) and a fixed value of \(\xi\). In addition, the previous methods are still relatively better in terms of \(SP\). While the methods with the lasso penalty (except SparseSVM) and logistic regression with EN penalty are better in terms of \(RC\) for any pair \((\rho,\xi)\). Finally, the methods are relatively comparable in terms of \(PR\).

Therefore, it can be concluded that our approach produced satisfactory and relatively similar results to HHSVM.

4.2 Empirical study↩︎

To illustrate the effectiveness of our BernSVM approach, we also compare classification accuracy, sensitivity and specificity of the latter methods on three large-scale samples of real data sets. The first data set is the DNA methylation measurements studied in [7], and two standard data sets, usually used to evaluate classifiers performance in the literature, namely, the prostate cancer [22] and the leukemia [23] data sets.
The DNA methylation measurements are collected around the BLK gene located in chromosome 8 to detect differentially methylated regions (DMRs, refer to genomic regions with significantly different methylation levels between two groups of samples, e.g. : cases-controls). The dataset contains measurements of DNA methylation levels of 40 different samples of cell types: B cells (8 samples), T cells (19 samples), and monocytes (13 samples). Methylation levels are measured in, \(p=5896\), CpG sites (predictors). Since methylation levels are known to be different between B-cell types compared to the T- and Monocyte-cell types around the BLK gene [7], we code the cell types as \(y=\{-1,1\}\) binary response, with \(y=1\) corresponds to B-cells and \(y=-1\) corresponds to T- and Monocyte-cell types. The second data-set is available in a publicly R package, spls. The prostate data consists of \(n = 102\) subjects, 52 patients with prostate tumor and 50 patients used as control subjects. The data contains expression of \(p = 6033\) genes across the subjects. We code the response variable \(y=\{-1,1\}\) binary response, with \(y=1\) corresponds to subjects with tumor prostate and \(y = -1\) corresponds to normal subjects. The third one contains \(n=72\) observations and \(p = 6817\) predictors and comes from a study of gene expression in two types of acute leukemias: acute lymphoblastic leukemia (ALL) and acute myeloid leukemia (AML). Gene expression levels were measured using Affymetrix high density oligonucleotide arrays. Subjects with AML are coded as \(y=-1\) class and subjects with ALL as \(y=1\) class. In total, we have 25 subjects in class 1 and 47 subjects in class -1. All the datasets were subjects to pre-screening procedure using a two sample t-tests [24] and only \(1000\) predictors are retained.

For the Elastic Net regularization, we take \(\lambda_2 = \alpha = 3/4\). The whole process was repeated 100 times and for each data-set we train the methods using 40\(\%\) of the observations to perform a \(10\) cross-validation to choose regularization parameters of each method and 60\(\%\) of the observations to test the methods and to compute the performance accuracy measures. The value of \(\delta\) is fixed to 2 for BernSVM and HHSVM.

We see from Table 4 that on the prostate data, the methods with the EN penalty provide a better result in terms of MR and SP (except SparseSVM). While in terms of SE, all methods provided a comparable result with little advantage to logistic with lasso penalty. For leukemia data, methods with EN penalty (except SparseSVM) dominate for all three performance measures. For methylation data, the methods with EN penalty produce the best result in terms of MR. While SparseSVM-EN followed by Logistic-EN dominate in terms of SE. More over all methods are relatively comparable in terms of SP, except SparseSVM-Lasso.

5 Discussion↩︎

In this work we aimed to better investigate the binary classification in high dimensional setting using the SVM classifier. We proposed a new smooth loss function, called BernSVM, to approximate the SVM hinge loss function. The BernSVM loss function has nice properties, which allows for efficient implementation of the penalized SVM and helps to derive appropriate theoretical results of the model parameter estimators. We have proposed two algorithms to solve the penalized BernSVM. The first one termed BSVM-GCD and combines the coordinate descent algorithm and the MM principle. The second one, called BSVM-IRLS, and uses an IRLS-type algorithm to solve the underlying optimization problem. We have also derived non asymptotic results of the penalized BernSVM estimator with weighted Lasso, with known weights, and we have extended this result for unknown weights. Furthermore, we have showed that the estimation error of penalized BernSVM with weighted Lasso achieves a rate of order \(\sqrt{s\log(p)/n}\). We compared our approach with its competitors through a simulation study and the results showed that our method outperforms its competitors in terms of the computational time while maintaining good performance in terms of prediction and variable selection. The proposed work have shown also accurate results when analyzing three high-dimensional real data sets. The penalized BernSVM is implemented in an R package BSVM, which is publicly available for use from Github (lien).
Finally, given the attractive properties of the BernSVM approach to approximate penalized SVM efficiently, one can adapt the proposed method for group penalties, such as group-Lasso/SCAD/MCP. Moreover, recently an interesting penalty has been developed [20], which fits linear regression models that splits the set of covariates into groups using a new diversity penalty. The authors provided also interesting properties of their proposed estimator. Adaptation of the diversity penalty to penalized SVM using the proposed BernSVM approach would be an interesting avenue to investigate. This is left for future work.

6 Appendix↩︎

6.1 Appendix \(\boldsymbol{A}\)↩︎

6.2 Appendix \(\boldsymbol{B}\)↩︎

6.3 Appendix \(\boldsymbol{C}\)↩︎

6.4 Appendix \(\boldsymbol{D}\)↩︎

6.5 Appendix \(\boldsymbol{E}\)↩︎

References↩︎