High-Dimensional Penalized Bernstein Support Vector Machines

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).

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.

Notice, first, that the coordinate objective function of problem (5) can be written as follows

$$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}),$$

(6)

where $r_{i}=y_{i}(\tilde{\beta}_{0}+\bm{x}_{i}^{\top}\tilde{\bm{\beta}})$, and $\tilde{\beta}_{0}$ and $\tilde{\bm{\beta}}$ are the current estimates of $\beta_{0}$ and $\bm{\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 (6) 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\neq j$, fixed at the current iteration

$$\underset{\beta_{j})}{\min}Q_{\delta}(\beta_{j}|\tilde{\beta_{0}},\tilde{\beta}_{j}),$$

(7)

where

$$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}^{{}^{\prime}}(r_{i})y_{i}x_{ij}}{n}(\beta_{j}-\tilde{\beta}_{j})+\frac{L}{2}(\beta_{j}-\tilde{\beta}_{j})^{2}+P_{\lambda_{1},\lambda_{2}}(\beta_{j}).$$

(8)

The next proposition gives the explicit solution of (7) for the different forms of $P_{\lambda_{1}}(\beta_{j})$.

The proof is outlined in Appendix A.

Algorithm 1 (BSVM-GCD), given next, gives details of steps for solving the objective function in (8) 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 1 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.

Firstly, we consider the non-penalized BernSVM with loss function

$$\ell{(\beta_{0},\bm{\beta})}=\frac{1}{n}\sum_{i=1}^{n}B_{\delta}(y_{i}(\beta_{0}+\bm{x}_{i}^{\top}\bm{\beta}))=\frac{1}{n}\sum_{i=1}^{n}B_{\delta}(y_{i}\bm{X}_{i}^{\top}\bm{b}),$$

(11)

where, $\bm{b}=(\beta_{0},\bm{\beta}^{\top})^{\top}$ and $\bm{X}=(1,\bm{x}^{\top})^{\top}$. To solve

$$\frac{\partial{\ell}}{\partial\bm{b}}(\beta_{0},\bm{\beta})=\frac{1}{n}\sum_{i=1}^{n}y_{i}\bm{X}_{i}B^{{}^{\prime}}_{\delta}(y_{i}\bm{X}_{i}^{\top}\bm{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\bm{b}\partial\bm{b}^{\top}}(\beta_{0},\bm{\beta})=\frac{1}{n}\sum_{i=1}^{n}\bm{X}_{i}B^{{}^{\prime\prime}}_{\delta}(y_{i}\bm{X}_{i}^{\top}\bm{b})\bm{X}_{i}^{\top}.$$

The update of the coefficients is given then by

$$\bm{b}^{new}=\tilde{\bm{b}}-(\frac{\partial^{2}{\ell}}{\partial\bm{b}\partial\bm{b}^{\top}}(\beta_{0},\bm{\beta}))^{-1}\frac{\partial{\ell}}{\partial\bm{b}}(\beta_{0},\bm{\beta})\tilde{\bm{b}}.$$

Let $u_{i}=y_{i}.B^{{}^{\prime}}_{\delta}(y_{i}\bm{X}_{i}^{\top}\tilde{\bm{b}})$, $\phi_{ii}=B^{{}^{\prime\prime}}_{\delta}(y_{i}\bm{X}_{i}^{\top}\tilde{\bm{b}})$ and $\tilde{\bm{\Phi}}=diag(\phi_{ii})$ for $i=1\ldots,n$. Then, the update of the Newton-Raphson algorithm can be written as follows,

	$\displaystyle\bm{b}^{new}$	$\displaystyle=$	$\displaystyle\tilde{\bm{b}}-(\bm{X}^{\top}\bm{\Phi}\bm{X})^{-1}\bm{X}^{\top}\bm{u}$
		$\displaystyle=$	$\displaystyle(\bm{X}^{\top}\bm{\Phi}\bm{X})^{-1}\bm{X}^{\top}\bm{\Phi}[\bm{X}\tilde{\bm{b}}-\bm{\Phi}^{-1}\bm{u}]$
		$\displaystyle=$	$\displaystyle(\bm{X}^{\top}\bm{\Phi}\bm{X})^{-1}\bm{X}^{\top}\bm{\Phi}\bm{z},$

where $\bm{z}=\bm{X}\tilde{\bm{b}}-\bm{\Phi}^{-1}\bm{u}$. Thus, $\bm{b}^{new}$ is the update of a weighted least regression problem with response $\bm{z}$, with the loss function

$$\ell{(\beta_{0},\bm{\beta})}\approx\frac{1}{2n}\sum_{i=1}^{n}\phi_{ii}(z_{i}-\bm{X}_{i}^{\top}\bm{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 (5), can be solved by combining IRLS with a cyclic coordinate descent. Thus, for the Adaptive Elastic Net penalty, the solution is given by

$$\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}},$$

(12)

where $r_{i}=z_{i}-\tilde{\beta}_{0}-\bm{x}_{i}^{\top}\tilde{\bm{\beta}}$ and $(\tilde{\beta}_{0},\tilde{\bm{\beta}})$ is the estimate of $({\beta}_{0},\bm{\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).

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

$$(\tilde{\beta}_{0},\tilde{\bm{\beta}})^{new}=\arg\min_{(\beta_{0},\bm{\beta})}\frac{1}{n}\sum_{i=1}^{n}B_{\delta}(y_{i}(\beta_{0}+\bm{x}_{i}^{\top}\bm{\beta}))+\lambda_{1}||\hat{\bm{W}}\bm{\beta}||_{1},$$

(13)

where $\hat{\bm{W}}=diag\{\hat{w}_{j}\}$, with $\hat{w}_{j}=P^{{}^{\prime}}_{\lambda_{1}}(\tilde{\beta_{j}})/\lambda_{1},j=1,\ldots,p,$ are the current weights, and $P^{{}^{\prime}}_{\lambda_{1}}(.)$ is the first derivative of the SCAD or MCP penalty given respectively by

$$P^{{}^{\prime}}_{\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^{{}^{\prime}}_{\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 3 below.