---
abstract: 'We consider classification tasks in the regime of scarce labeled training data in high dimensional feature space, where specific expert knowledge is also available. We propose a new hybrid optimization algorithm that solves the elastic-net support vector machine (SVM) through an alternating direction method of multipliers in the first phase, followed by an interior-point method for the classical SVM in the second phase. Both SVM formulations are adapted to knowledge incorporation. Our proposed algorithm addresses the challenges of automatic feature selection, high optimization accuracy, and algorithmic flexibility for taking advantage of prior knowledge. We demonstrate the effectiveness and efficiency of our algorithm and compare it with existing methods on a collection of synthetic and real-world data.'
author:
- 'Z. Qin'
- 'X. Tang'
- 'I. Akrotirianakis'
- 'A. Chakraborty'
bibliography:
- 'scr\_bib.bib'
title: 'HIPAD - A Hybrid Interior-Point Alternating Direction algorithm for knowledge-based SVM and feature selection'
---

\[section\] \[section\] \[section\] \[section\] \[section\]

Introduction {#sec:motivation}
============

Classification tasks on data sets with large feature dimensions are very common in real-world machine learning applications. Typical examples include microarray data for gene selection and text documents for natural language processing. Despite the large number of features present in the data sets, usually only small subsets of the features are relevant to the particular learning tasks, and local correlation among the features is often observed. Hence, feature selection is required for good model interpretability. Popular classification techniques, such as support vector machine (SVM) and logistic regression, are formulated as convex optimization problems. An extensive literature has been devoted to optimization algorithms that solve variants of these classification models with sparsity regularization [@pardalos-hansen-2008; @sra-nowozin-wright-2011]. Many of them are based on first-order (gradient-based) methods, mainly because the size of the optimization problem is very large. The advantage of first-order methods is that their computational and memory requirements at each iteration are low and as a result they can handle the large optimization problems occurring in classification problems. Their major disadvantage is their slow convergence, especially when a good approximation of the feature support has been identified. Second-order methods exhibit fast local convergence, but their computational and memory requirements are much more demanding, since they need to store and invert the Newton matrix at every iteration. It is therefore very important to be able to intelligently combine the advantages of both the first and the second order optimization methods in such a way that the resulting algorithm can solve large classification problems efficiently and accurately. As we will demonstrate in this paper such combination is possible by taking advantage of the problem structure and the change in its size during the solution process. In addition, we will also show that our algorithmic framework is flexible enough to incorporate prior knowledge to improve classification performance.

Related Work {#sec:related_work}
------------

The above requirements demand three features from a learning algorithm: 1. it should be able to automatically select features which are possibly in groups and highly correlated; 2. it has to solve the optimization problem in the training phase efficiently and with high accuracy; and 3. the learning model needs to be flexible enough so that domain knowledge can be easily incorporated. Existing methods are available in the literature that meet some of the above requirements [*individually*]{}. For enforcing sparsity in the solution, efficient optimization algorithms such as that proposed in [@koh2007interior] can solve large-scale sparse logistic regression. On the other hand, the $L_1$-regularization is unstable with the presence of highly correlated features - among a group of such features, essentially one of them is selected in a random manner. To handle local correlation among groups of features, the elastic-net regularization [@zou2005regularization] has been successfully applied to SVM [@wang2006doubly] and logistic regression [@ryali2010sparse]. However, incorporating domain knowledge into the logistic regression formulation is not straightforward. For SVM, including such knowledge in the optimization process has been demonstrated in [@fung2003knowledge]. Recently, an Alternating Direction Method of Multipliers (ADMM) has been proposed for the elastic-net SVM (ENSVM) [@yeefficient]. ADMM is quick to find an approximate solution to the ENSVM problem, but it is known to converge very slowly to high accuracy optimal solutions [@boyd2010distributed]. The interior-point methods (IPM) for SVM are known to be able to achieve high accuracy in their solutions with a polynomial iteration complexity, and the dual SVM formulation is independent of the feature space dimensionality. However, the classical $L_2$-norm SVM is not able to perform automatic feature selection. Although the elastic-net SVM can be formulated as a QP (in the primal form), its problem size grows substantially with the feature dimensionality. Due to the need to solve a Newton system in each iteration, the efficiency of IPM quickly deteriorates as the feature dimension becomes large.

Main Contributions
------------------

In this paper we propose a new hybrid algorithmic framework for SVM to address [*all*]{} of the above challenges and requirements [*simultaneously*]{}. Our framework combines the advantages of a first-order optimization algorithm (through the use of ADMM) and a second-order method (via IPM) to achieve both superior speed and accuracy. Through a novel algorithmic approach that is able to incorporate expert knowledge, our proposed framework is able to exploit domain knowledge to improve feature selection, and hence, prediction accuracy. Besides efficiency and generalization performance, we demonstrate through experiments on both synthetic and real data that our method is also more robust to inaccuracy in the supplied knowledge than existing approaches.

A Two-phase Hybrid Optimization Algorithm {#sec:two_phase_alg}
=========================================

As previously mentioned, for data sets with many features, the high dimensionality of the feature space still poses a computational challenge for IPM. Fortunately, many data sets of this kind are very sparse, and the resulting classifier $w$ is also expected to be sparse, i.e. only a small subset of the features are expected to carry significant weights in classification. Naturally, it is ideal for IPM to train a classifier on the most important features only.

Inspired by the Hybrid Iterative Shrinkage (HIS) [@shi2010fast] algorithm for training large-scale sparse logistic regression classifiers, we propose a two-phase algorithm to shrink the feature space appropriately so as to leverage the high accuracy of IPM while maintaining efficiency. Specifically, we propose to solve an elastic-net SVM (ENSVM) or doubly-regularized SVM (DrSVM) [@wang2006doubly] problem during the first phase of the algorithm. The elastic-net regularization performs feature selection with grouping effect and has been shown to be effective on data sets with many but sparse features and high local correlations [@zou2005regularization]. This is the case for text classification, microarray gene expression, and fMRI data sets. The support of the weight vector $w$ for ENSVM usually stabilizes well before the algorithm converges to the optimal solution. Taking advantage of that prospect, we can terminate the first phase of the hybrid algorithm early and proceed to solve a classical SVM problem with the reduced feature set in the second phase, using an IPM solver.

Solving the Elastic Net SVM using ADMM {#sec:admm}
--------------------------------------

SVM can be written in the regularized regression form as $$\label{eq:unconstr_svm}
    \min_{{\textbf{w}},b}\frac{1}{N}\sum_{i=1}^N(1 - (y_i({\textbf{x}_i}^T{\textbf{w}}+ b)))_+ + \frac{\lambda}{2}\|{\textbf{w}}\|_2^2,$$ where the first term is an averaged sum of the hinge losses and the second term is viewed as a ridge regularization on $w$. It is easy to see from this form that the classical SVM does not enforce sparsity in the solution, and $w$ is generally dense. The ENSVM adds an $L_1$ regularization on top of the ridge regularization term, giving $$\label{eq:en_svm}
    \min_{{\textbf{w}},b}\frac{1}{N}\sum_{i=1}^{N}(1-y_i({\textbf{x}_i}^T{\textbf{w}}+ b))_+ + \lambda_1\|{\textbf{w}}\|_1 + \frac{\lambda_2}{2}\|{\textbf{w}}\|_2^2.$$ Compared to the Lasso ($L_1$-regularized regression) [@tibshirani1996regression], the elastic-net has the advantage of selecting highly correlated features in groups (i.e. the grouping effect) while still enforcing sparsity in the solution. This is a particularly attractive feature for text document data, which is common in the hierarchical classification setting. Adopting the elastic-net regularization as in (\[eq:en\_svm\]) brings the same benefit to SVM for training classifiers.

To approximately solve problem (\[eq:en\_svm\]), we adopt the alternating direction method of multipliers (ADMM) for elastic-net SVM recently proposed in [@yeefficient]. ADMM has a long history dating back to the 1970s [@gabay1976dual]. Recently, it has been successfully applied to problems in machine learning [@boyd2010distributed]. ADMM is a special case of the inexact augmented Lagrangian (IAL) method [@rockafellar1973multiplier] for the structured unconstrained problem $$\label{eq:struct_unc}
    \min_x F(x) \equiv f(x) + g(Ax),$$ where both functions $f(\cdot)$ and $g(\cdot)$ are convex. We can decouple the two functions by introducing an auxiliary variable $y$ and convert problem (\[eq:struct\_unc\]) into an equivalent constrained optimization problem $$\begin{aligned}
\label{eq:struct_c}
  \min_{x,y} && f(x) + g(y), \quad s.t.  \>\> Ax = y.\end{aligned}$$ This technique is often called variable-splitting [@combettes2011proximal]. The IAL method approximately minimizes in each iteration the augmented Lagrangian of (\[eq:struct\_c\]) defined by $\mathcal{L}(x,y,\gamma) := f(x) + g(y) + \gamma^T(y-Ax) + \frac{\mu}{2}\|Ax-y\|_2^2$ , followed by an update to the Lagrange multiplier $\gamma \gets \gamma - \mu(Ax - y)$. The IAL method is guaranteed to converge to the optimal solution of (\[eq:struct\_unc\]), as long as the subproblem of approximately minimizing the augmented Lagrangian is solved with an increasing accuracy [@rockafellar1973multiplier]. ADMM can be viewed as a practical implementation of IAL, where the subproblem is solved approximately by minimizing $\mathcal{L}(x,y;\gamma)$ with respect to $x$ and $y$ [*alternatingly once*]{}. Eckstein and Bertsekas [@eckstein1992douglas] established the convergence of ADMM for the case of two-way splitting. Now applying variable-splitting and ADMM to problem (\[eq:en\_svm\]), [@yeefficient] introduced auxiliary variables $({\textbf{a}},{\textbf{c}})$ and linear constraints so that the non-smooth hinge loss and $L_1$-norm in the objective function are decoupled, making it easy to optimize over each of the variables. Specifically, problem (\[eq:en\_svm\]) is transformed into an equivalent constrained form $$\begin{aligned}
\label{eq:en_svm_constr}
    \min_{{\textbf{w}},b,{\textbf{a}},{\textbf{c}}} && \frac{1}{N}\sum_{i=1}^N (a_i)_+ + \lambda_1\|\textbf{c}\|_1 + \frac{\lambda_2}{2}\|{\textbf{w}}\|_2^2 \\
    \nonumber s.t. && \textbf{a} = {\textbf{e}}- Y(X{\textbf{w}}+ b{\textbf{e}}) \quad \mbox{and} \quad {\textbf{c}}= {\textbf{w}}$$ where ${\textbf{x}_i}^T$ is the $i$-th row of $X$, and $Y = \textrm{diag}({\textbf{y}})$. The augmented Lagrangian $$\begin{gathered}
    \mathcal{L}({\textbf{w}},b,{\textbf{a}},{\textbf{c}},\gamma_1,\gamma_2) := \frac{1}{N}\sum_{i=1}^N a_i + \lambda_1\|\textbf{c}\|_1 + \frac{\lambda_2}{2}\|{\textbf{w}}\|_2^2  + \gamma_1^T({\textbf{e}}- Y(X{\textbf{w}}+ b{\textbf{e}})-{\textbf{a}})\\
    +\gamma_2^T({\textbf{w}}-{\textbf{c}})  + \frac{\mu_1}{2}\|{\textbf{e}}- Y(X{\textbf{w}}+ b{\textbf{e}})-{\textbf{a}}\|_2^2 + \frac{\mu_2}{2}\|{\textbf{w}}-{\textbf{c}}\|_2^2\end{gathered}$$ is then minimized with respect to $({\textbf{w}},b), {\textbf{a}},$ and ${\textbf{c}}$ sequentially in each iteration, followed by an update to the Lagrange multipliers $\gamma_1$ and $\gamma_2$. The original problem is thus decomposed into three subproblems consisting of computing the proximal operator of the hinge loss function (with respect to ${\textbf{a}}$), solving a special linear system (with respect to $({\textbf{w}},b)$ ), and performing a soft-thresholding operation (with respect to ${\textbf{c}}$), which can all be done in an efficient manner. Due to lack of space in the paper, we have included the detailed solution steps in the Appendix (see Algorithm \[alg:admm-ensvm\] ADMM-ENSVM), where we define by $\mathcal{S}_\lambda(\cdot)$ the proximal operator associated with the hinge loss $$\begin{aligned}
  \mathcal{S}_\lambda(\omega) &=& \left\{
                                    \begin{array}{ll}
                                      \omega - \lambda, & \hbox{$\omega > \lambda$;} \\
                                      0, & \hbox{$0 \leq \omega \leq \lambda$;} \\
                                      \omega, & \hbox{$\omega < 0$.}
                                    \end{array}
                                  \right.\end{aligned}$$ and $\mathcal{T}_\lambda(\omega) = sgn(\omega)\max\{0,|\omega|-\lambda \}$ is the shrinkage operator.

SVM via Interior-Point Method {#sec:ipm}
-----------------------------

Interior Point Methods enjoy fast convergence rates for a wide class of QP problems. Their theoretical polynomial convergence ($O(n\log \frac{1}{\epsilon})$) was first established by Mizuno [@mizuno1994]. In addition, Andersen [*et al*]{} [@Andersen1996implement] showed that the number of iterations needed by IPMs to converge is $O(\log n)$, which demonstrates that their computational effort increases in a slower rate than the size of the problem.

Both the primal and the dual SVM are QP problems. The primal formulation of SVM [@vapnik2000nature] is defined as $$\begin{aligned}
  \textrm{(SVM-P)} \qquad   \min_{{\textbf{w}},b,\xi,{\textbf{s}}} && \frac{1}{2}{\textbf{w}}^T{\textbf{w}}+ c{\textbf{e}}^T\xi \\
  \qquad s.t. && y_i({\textbf{w}}^T{\textbf{x}_i}- b) + \xi_i - s_i = 1, i=1,\dots,N, \\
  && {\textbf{s}}\geq 0, \xi \geq 0.\end{aligned}$$ whereas the dual SVM has the form $$\begin{aligned}
  \textrm{(SVM-D)} \quad \min_{\alpha} && \frac{1}{2}\alpha^TQ\alpha - {\textbf{e}}^T\alpha \\
  s.t. && {\textbf{y}}^T\alpha = 0, \quad \mbox{and} \quad 0 \leq \alpha_i \leq c, \quad i = 1,\cdots,N,\end{aligned}$$ where $Q_{ij} = y_iy_j{\textbf{x}_i}^T{\textbf{x}_j}= \bar{X}\bar{X}^T$. By considering the KKT conditions of (SVM-P) and (SVM-D), the optimal solution is given by ${\textbf{w}}= \bar{X}^T\alpha = \sum_{i\in SV}\alpha_i y_i {\textbf{x}_i}$, where $SV$ is the set of sample indices corresponding to the support vectors. The optimal bias term $b$ can be computed from the complementary slackness condition $\alpha_i(y_i({\textbf{x}_i}^T{\textbf{w}}+ b)-1 + \xi_i) = 0$.

Whether to solve (SVM-P) or (SVM-D) for a given data set depends on its dimensions as well as its sparsity. Even if $X$ is a sparse matrix, $Q$ in (SVM-D) is still likely to be dense, whereas the Hessian matrix in (SVM-P) is the identity. The primal problem (SVM-P), however, has a larger variable dimension and more constraints. It is often argued that one should solve (SVM-P) when the number of features is smaller than the number of samples, whereas (SVM-D) should be solved when the number of features is less than that of the samples. Since in the second phase of Algorithm \[alg:hipad\] we expect to have identified a small number of promising features, we have decided to solve (SVM-D) by using IPM. Solving (SVM-D) is realized through the OOQP [@gertz2003object] software package that implements a primal-dual IPM for convex QP problems.

The Two-phase Algorithm
-----------------------

Let us keep in mind that the primary objective for the first phase is to appropriately reduce the feature space dimensionality without impacting the final prediction accuracy. As we mentioned above, the suitability of ADMM for the first phase depends on whether the support of the feature vector converges quickly or not. On an illustrative dataset from [@yeefficient], which has 50 samples with 300 features each, ADMM converged in 558 iterations. The output classifier ${\textbf{w}}$ contained only 13 non-zero features, and the feature support converged in approximately 50 iteration (see Figure \[fig:fea\_supp\_plots\] in the Appendix for illustrative plots showing the early convergence of ADMM). Although the remaining more than 500 iterations are needed by ADMM in order to satisfy the optimality criteria, they do not offer any additional information regarding the feature selection process. Hence, it is important to monitor the change in the signs and indices of the support and terminate the first phase promptly. In our implementation, we adopt the criterion used in [@shi2010fast] and monitor the relative change in the iterates as a surrogate of the change in the support, i.e. $$\label{eq:hipad_transition}
    \frac{\|{\textbf{w}}^{k+1}-{\textbf{w}}^k\|}{\max(\|{\textbf{w}}^k\|,1)} < \epsilon_{tol}.$$ We have observed in our experiments that when the change over the iterates is small, the evolution of the support indices stabilizes too.

Upon starting the second phase, it is desirable for IPM to warm-start from the corresponding sub-vector of the solution returned by ADMM. It should also be noted that we apply IPM during the second phase to solve the classical $L_2$-regularized SVM (\[eq:unconstr\_svm\]), instead of the ENSVM (\[eq:en\_svm\]) in the first phase. There are two main reasons for this decision. First, although ENSVM can be reformulated as a QP, the size of the problem is larger than the classical SVM due to the additional linear constraints introduced by the $L_1$-norm. Second, since we have already identified (approximately) the feature support in the first phase of the algorithm, enforcing sparsity in the reduced feature space becomes less critical. The two-phase algorithm is summarized in Algorithm \[alg:hipad\].

\[alg:hipad\] Given ${\textbf{w}}^0, b^0, {\textbf{a}}^0, {\textbf{c}}^0, {\textbf{u}}^0,$ and ${\textbf{v}}^0$. **PHASE 1: ADMM for ENSVM** $({\textbf{w}}^{\textrm{ADMM}},b^{\textrm{ADMM}}) \gets \textrm{ADMM-ENSVM}({\textbf{w}}^0, b^0, {\textbf{a}}^0, {\textbf{c}}^0, {\textbf{u}}^0,{\textbf{v}}^0)$ **PHASE 2: IPM for SVM** $\widetilde{{\textbf{w}}} \gets $ non-zero components of ${\textbf{w}}^{\textrm{ADMM}}$ $(\widetilde{X}, \widetilde{Y}) \gets $ sub-matrices of $(X, Y)$ corresponding to the support of ${\textbf{w}}^{\textrm{ADMM}}$ $({\textbf{w}},b) \gets \textrm{SVM-IPM}(\widetilde{X}, \widetilde{Y}, \widetilde{{\textbf{w}}}$), through (SVM-D). $({\textbf{w}},b)$

Domain Knowledge Incorporation {#sec:enk-svm}
==============================

Very often, we have prior domain knowledge for specific classification tasks. Domain knowledge is most helpful when the training data does not form a comprehensive representation of the underlying unknown population, resulting in poor generalization performance of SVM on the unseen data from the same population. This often arises in situations where labeled training samples are scarce, while there is an abundance of unlabeled data.

For high dimensional data, ENSVM performs feature selection along with training to produce a simpler model and to achieve better prediction accuracy. However, the quality of the feature selection depends entirely on the training data. In pathological cases, it is very likely that the feature support identified by ENSVM does not form a good representation of the population. Hence, when domain knowledge about certain features is available, we should take it into consideration during the training phase and include the relevant features in the feature support should them be deemed important for classification.

In this section, we explore and propose a new approach to achieve this objective. We consider domain knowledge in the form of class-membership information associated with features. We can incorporate such information (or enforce such rules) in SVM by adding equivalent linear constraints to the SVM QP problem (KSVM) [@fung2003knowledge; @lauer2008incorporating]. To be specific, we can model the above information with the linear implication $$\label{eq:linear_implication}
    B{\textbf{x}}\leq \textbf{d} \quad \Rightarrow \quad {\textbf{w}}^T{\textbf{x}}+ b \geq 1,$$ where $B \in \mathbb{R}^{k_1\times m}$ and ${\textbf{d}}\in \mathbb{R}^{k_1}$. It is shown in [@fung2003knowledge] that by utilizing the non-homogeneous Farkas theorem of the alternative, (\[eq:linear\_implication\]) can be transformed into the following equivalent system of linear inequalities with a solution ${\textbf{u}}$ $$\label{eq:equiv_lin_constr}
    B^T{\textbf{u}}+ {\textbf{w}}= \textbf{0}, \quad \textbf{d}^T{\textbf{u}}- b + 1 \leq 0, \quad \textbf{u} \geq \textbf{0}.$$ Similarly, for the linear implication for the negative class membership we have: $$\label{eq:lin_imp_neg}
    D{\textbf{x}}\leq {\textbf{g}}\Rightarrow  {\textbf{w}}^T{\textbf{x}}+ b \leq -1, \quad D \in \mathbb{R}^{k_2\times m}, g\in\mathbb{R}^{k_2},$$ which can be represented by the set of linear constraints in ${\textbf{v}}$ $$\label{eq:equiv_lin_constr_neg}
    D^T{\textbf{v}}- {\textbf{w}}= {\textbf{0}}, \quad {\textbf{g}}^T{\textbf{v}}+ b + 1 \leq 0, \quad {\textbf{v}}\geq {\textbf{0}}.$$ Hence, to incorporate the domain knowledge represented by (\[eq:linear\_implication\]) and (\[eq:lin\_imp\_neg\]) into SVM, Fung [*et al*]{} [@fung2003knowledge] simply add the linear constraints (\[eq:equiv\_lin\_constr\]) and (\[eq:equiv\_lin\_constr\_neg\]) to (SVM-P). Their formulation, however, increases both the variable dimension and the number of linear constraints by at least $2m$, where $m$ is the number of features in the classification problem we want to solve. This is clearly undesirable when $m$ is large, which is the scenario that we consider in this paper.

In order to avoid the above increase in the size of the optimization probelm, we choose to penalize the quadratics $\|B^T{\textbf{u}}+{\textbf{w}}\|_2^2$ and $\|D^T{\textbf{v}}-{\textbf{w}}\|_2^2$ instead of their $L_1$ counterparts considered in [@fung2003knowledge]. By doing so the resulting problem is still a convex QP but with a much smaller size. Hence, we consider the following model for domain knowledge incorporation. $$\begin{aligned}
   \textrm{(KSVM-P)} \qquad
    \min_{{\textbf{w}},b,\xi,{\textbf{u}},{\textbf{v}},\eta_u,\eta_v} && \frac{1}{2}{\textbf{w}}^T{\textbf{w}}+ c{\textbf{e}}^T\xi + \frac{\rho_1}{2}\|B^T{\textbf{u}}+{\textbf{w}}\|_2^2 \\
    && + \rho_2\eta_u + \frac{\rho_3}{2}\|D^T{\textbf{v}}-{\textbf{w}}\|_2^2 + \rho_4\eta_v \\
  \nonumber s.t. && y_i({\textbf{w}}^T{\textbf{x}_i}+ b) \geq 1 - \xi_i, \quad i = 1,\cdots,N, \\
   \nonumber && \textbf{d}^T{\textbf{u}}- b + 1 \leq \eta_u, \\
   \nonumber && \textbf{g}^T{\textbf{v}}+ b + 1 \leq \eta_v,\\
   \nonumber && (\xi,{\textbf{u}},{\textbf{v}},\eta_u,\eta_v) \geq {\textbf{0}}.\end{aligned}$$

We are now ready to propose a novel combination of ENSVM and KSVM, and we will explain in the next section how the combined problem can be solved in our HIPAD framework. The main motivation behind this combination is to exploit domain knowledge to improve the feature selection, and hence, the generalization performance of HIPAD. To the best of our knowledge, this is the first method of this kind.

ADMM Phase
----------

Our strategy for solving the elastic-net SVM with domain knowledge incorporation is still to apply the ADMM method. First, we combine problems (\[eq:en\_svm\]) and (KSVM-P) and write the resulting optimization problem in an equivalent unconstrained form (by penalizing the violation of the inequality constraints through hinge losses in the objective function) $$    \textrm{(ENK-SVM)} \quad \min_{{\textbf{w}},b,{\textbf{u}}\geq {\textbf{0}},{\textbf{v}}\geq {\textbf{0}}} F({\textbf{w}},b,{\textbf{u}},{\textbf{v}}),$$ where $F({\textbf{w}},b,{\textbf{u}},{\textbf{v}}) \equiv \frac{\lambda_2}{2}\|{\textbf{w}}\|_2^2 + \lambda_1\|{\textbf{w}}\|_1 + \frac{1}{N}\sum_{i=1}^{N}(1-y_i(x_i^Tw + b))_+
    + \frac{\rho_1}{2}\|B^T{\textbf{u}}+{\textbf{w}}\|_2^2 + \rho_2({\textbf{d}}^T{\textbf{u}}- b + 1)_+ + \frac{\rho_3}{2}\|D^T{\textbf{v}}-{\textbf{w}}\|_2^2
    + \rho_4({\textbf{g}}^T{\textbf{v}}+ b + 1)_+$. We then apply variable-splitting to decouple the $L_1$-norms and hinge losses and obtain the following equivalent constrained optimization problem: $$\begin{aligned}
\label{eq:ensvm_ksvm_c}
  \min_{{\textbf{w}},b,{\textbf{u}},{\textbf{v}},{\textbf{a}},{\textbf{c}},p,q} &&\>\>\>  F({\textbf{w}},b,{\textbf{u}},{\textbf{v}},{\textbf{a}},{\textbf{c}},p,q) \\
  \nonumber s.t. && \>\>\> {\textbf{a}}= {\textbf{e}}- (\bar{X}{\textbf{w}}+ {\textbf{y}}b), \>\>\>  {\textbf{c}}= {\textbf{w}}, \\
  \nonumber && \>\>\> q = {\textbf{d}}^T{\textbf{u}}- b + 1,   \>\>\> p = {\textbf{g}}^T{\textbf{v}}+ b + 1, \\
  \nonumber && \>\>\> {\textbf{u}}\geq {\textbf{0}}, \>\>\> {\textbf{v}}\geq {\textbf{0}}.\end{aligned}$$ with $F({\textbf{w}},b,{\textbf{u}},{\textbf{v}},{\textbf{a}},{\textbf{c}},p,q)
  \equiv \frac{\lambda_2}{2}\|{\textbf{w}}\|_2^2 + \lambda_1\|{\textbf{c}}\|_1 + \frac{1}{N}{\textbf{e}}^T({\textbf{a}})_+
    + \frac{\rho_1}{2}\|B^T{\textbf{u}}+{\textbf{w}}\|_2^2 + \rho_2(q)_+
+ \frac{\rho_3}{2}\|D^T{\textbf{v}}-{\textbf{w}}\|_2^2 + \rho_4(p)_+$. As usual, we form the augmented Lagrangian $\mathcal{L}$ of problem (\[eq:ensvm\_ksvm\_c\]), $$\begin{gathered}
    \mathcal{L} := F({\textbf{w}},b,{\textbf{u}},{\textbf{v}},{\textbf{a}},{\textbf{c}},p,q)
    + \gamma_1^T({\textbf{e}}-({\bar{X}}{\textbf{w}}+{\textbf{y}}b)-{\textbf{a}})
    + \frac{\mu_1}{2}\|{\textbf{e}}-({\bar{X}}+{\textbf{y}}b)-{\textbf{a}}\|_2^2\\
    + \gamma_2^T({\textbf{w}}-{\textbf{c}}) + \frac{\mu_2}{2}\|{\textbf{w}}-{\textbf{c}}\|_2^2
    + \gamma_3({\textbf{d}}^T{\textbf{u}}-b+1-q) + \frac{\mu_3}{2}\|{\textbf{d}}^T{\textbf{u}}-b+1-q\|_2^2 \\
    + \gamma_4({\textbf{g}}^T{\textbf{v}}+b+1-p) + \frac{\mu_4}{2}\|{\textbf{g}}^T{\textbf{v}}+b+1-p\|_2^2\end{gathered}$$ and minimize $\mathcal{L}$ with respect to ${\textbf{w}},b,{\textbf{c}},{\textbf{a}},p,q,{\textbf{u}},{\textbf{v}}$ individually and in order. For the sake of readability, we do not penalize the non-negative constraints for ${\textbf{u}}$ and ${\textbf{v}}$ in the augmented Lagrangian.

Given $({\textbf{a}}^k, {\textbf{c}}^k, p^k, q^k)$, solving for $({\textbf{w}},b)$ again involves solving a linear system $$\label{eq:ksvm_admm_linsys_wb}
        \left(
          \begin{array}{cc}
            \kappa_1 I + \mu_1 X^TX & \mu_1 X^T{\textbf{e}}\\
            \mu_1{\textbf{e}}^TX & \mu_1 N + \kappa_2 \\
          \end{array}
        \right)
         \left(
          \begin{array}{c}
            {\textbf{w}}^{k+1} \\
            b^{k+1} \\
          \end{array}
        \right) =
        \left(
          \begin{array}{c}
             \textbf{r}_{{\textbf{w}}}\\
             \textbf{r}_b\\
          \end{array}
        \right),$$ where $\kappa_1 = \lambda_2 + \mu_2 + \rho_1 + \rho_3, \kappa_2 = \mu_3 + \mu_4, \textbf{r}_{{\textbf{w}}} = X^TY\gamma_1^k + \mu_1 X^TY({\textbf{e}}-{\textbf{a}}^k) - \gamma_2^k + \mu_2{\textbf{c}}^k + \rho_3D^T{\textbf{v}}^k - \rho_1B^T{\textbf{u}}^k$ and $\textbf{r}_{b} = {\textbf{e}}^TY\gamma_1^k + \mu_1{\textbf{e}}^TY({\textbf{e}}-{\textbf{a}}^k) + \gamma_3^k + \mu_3({\textbf{d}}^T{\textbf{u}}^k + 1 - q^k) - \gamma_4^k - \mu_4({\textbf{g}}^T{\textbf{v}}^k+1-p^k)$. Similar to solving the linear system in Algorithm \[alg:admm-ensvm\] ADMM-ENSVM, we can compute the solution to the above linear system through a few PCG iterations, taking advantage of the fact that the left-hand-side matrix is of low-rank.

To minimize the augmented Lagrangian with respect to ${\textbf{u}}$, we need to solve a convex quadratic problem with non-negative constraints $$\label{eq:ksvm_admm_usub}
    \min_{{\textbf{u}}\geq {\textbf{0}}} \quad \frac{\rho_1}{2}\|B^T{\textbf{u}}+ {\textbf{w}}^{k+1}\|_2^2 + \gamma_3^k{\textbf{d}}^T{\textbf{u}}+ \frac{\mu_3}{2}\|{\textbf{d}}^T{\textbf{u}}- b^{k+1} + 1 - q^{k}\|_2^2.$$ Solving problem efficiently is crucial for the efficiency of the overal algorithm. We describe a novel way to do so. Introducing a slack variable ${\textbf{s}}$ and transferring the non-negative constraint on ${\textbf{u}}$ to ${\textbf{s}}$, we decompose the problem into two parts which are easy to solve. Specifically, we reformulate (\[eq:ksvm\_admm\_usub\]) as $$\begin{aligned}
    \min_{{\textbf{u}},{\textbf{s}}\geq {\textbf{0}}} && \>\>\> \frac{\rho_1}{2}\|B^T{\textbf{u}}+ {\textbf{w}}^{k+1}\|_2^2 + \gamma_3^k{\textbf{d}}^T{\textbf{u}}+ \frac{\mu_3}{2}\|{\textbf{d}}^T{\textbf{u}}- b^{k+1} + 1 - q^{k}\|_2^2 \\
    \nonumber s.t. && \>\>\>  {\textbf{u}}- {\textbf{s}}= {\textbf{0}}.\end{aligned}$$ Penalizing the linear constraint ${\textbf{u}}- {\textbf{s}}= {\textbf{0}}$ in the new augmented Lagrangian, the new subproblem with respect to $({\textbf{u}},{\textbf{s}})$ is $$\begin{gathered}
\label{eq:ksvm_admm_ussub}
    \min_{{\textbf{u}},{\textbf{s}}\geq {\textbf{0}}} \quad \frac{\rho_1}{2}\|B^T{\textbf{u}}+ {\textbf{w}}^{k+1}\|_2^2 + \gamma_3^k{\textbf{d}}^T{\textbf{u}}\\
    + \frac{\mu_3}{2}\|{\textbf{d}}^T{\textbf{u}}- b^{k+1} + 1 - q^k\|_2^2 + \gamma_5^T({\textbf{s}}-{\textbf{u}}) + \frac{\mu_5}{2}\|{\textbf{u}}-{\textbf{s}}\|_2^2.\end{gathered}$$ Given an ${\textbf{s}}^k \geq {\textbf{0}}$, we can compute ${\textbf{u}}^{k+1}$ by solving a $k_1 \times k_1$ linear system $$\label{eq:ksvm_admm_linsys_u}
    (\rho_1 BB^T + \mu_3 {\textbf{d}}{\textbf{d}}^T + \mu_5 I){\textbf{u}}^{k+1} = \textbf{r}_{\textbf{u}},$$ where $\textbf{r}_{\textbf{u}}= -\rho_1B{\textbf{w}}^{k+1} + \mu_3 {\textbf{d}}b^{k+1} + \mu_3{\textbf{d}}(q^k-1) - {\textbf{d}}\gamma_3^k + \gamma_5 + \mu_5{\textbf{s}}^k$. We assume that $B$ has full row-rank. This is a reasonable assumption since otherwise there is at least one redundant domain knowledge constraint and we can simply remove it. The number of domain knowledge constraints ($k_1$ and $k_2$) are usually small, so the system (\[eq:ksvm\_admm\_linsys\_u\]) can be solved exactly and efficiently by Cholesky factorization.

Solving for ${\textbf{s}}^{k+1}$ corresponding to ${\textbf{u}}^{k+1}$ is easy, observing that problem (\[eq:ksvm\_admm\_ussub\]) is separable in the elements of ${\textbf{s}}$. For each element $s_i$, the optimal solution to the one-dimensional quadratic problem with a non-negative constraint on $s_i$ is given by $\max(0,u_i - \frac{(\gamma_5)_i}{\mu_5})$. Writing in the vector form, ${\textbf{s}}^{k+1} = \max({\textbf{0}}, {\textbf{u}}^{k+1} - \frac{\gamma_5^k}{\mu_5})$. Similarly, we solve for ${\textbf{v}}^{k+1}$ by introducing a non-negative slack variable ${\textbf{t}}$ and solve the linear system $$\label{eq:ksvm_admm_linsys_v}
  (\rho_3DD^T + \mu_4 {\textbf{g}}{\textbf{g}}^T + \mu_6 I){\textbf{v}}^{k+1} = {\textbf{r}}_{\textbf{v}},$$ where ${\textbf{r}}_{\textbf{v}}= \rho_3D{\textbf{w}}^{k+1} - \mu_4{\textbf{g}}b^{k+1} - {\textbf{g}}\gamma_4^k - \mu_4{\textbf{g}}(1-p^k) + \gamma_6 + \mu_6{\textbf{t}}^k$, and ${\textbf{t}}^{k+1} = \max(0,{\textbf{v}}^{k+1} - \frac{\gamma_6^k}{\mu_6})$.

Now given $({\textbf{w}}^{k+1},b^{k+1},{\textbf{u}}^{k+1},{\textbf{v}}^{k+1})$, the solutions for ${\textbf{a}}$ and ${\textbf{c}}$ are exactly the same as in Lines \[line:a\_sol\] and \[line:c\_sol\] of Algorithm \[alg:hipad\], i.e. $$\begin{aligned}
  {\textbf{a}}^{k+1} &=& \mathcal{S}_{\frac{1}{N\mu_1}}\left( {\textbf{e}}+ \frac{\gamma_1^k}{\mu_1} - Y(X{\textbf{w}}^{k+1} + b^{k+1}{\textbf{e}}) \right), \\
  {\textbf{c}}^{k+1} &=& \mathcal{T}_{\frac{\lambda_1}{\mu_2}}\left( \frac{\gamma_2^k}{\mu_2}+{\textbf{w}}^{k+1} \right).\end{aligned}$$ The subproblem with respect to $q$ is $$\begin{aligned}
\label{eq:ksvm_admm_qsub}
    \nonumber \min_q && \quad \rho_2(q)_+ - \gamma_3^k q  + \frac{\mu_3}{2}\|{\textbf{d}}^T{\textbf{u}}^k - b^{k+1} + 1 - q\|_2^2 \equiv \\
     && \rho_2(q)_+
     + \frac{\mu_3}{2}\|q - ({\textbf{d}}^T{\textbf{u}}^k - b^{k+1} + 1 + \frac{\gamma_3^k}{\mu_3})\|_2^2.\end{aligned}$$ The solution is given by a (one-dimensional) proximal operator associated with the hinge loss $$\label{eq:ksvm_admm_qsol}
    q^{k+1} = \mathcal{S}_{\frac{\rho_2}{\mu_3}}\left( {\textbf{d}}^T{\textbf{u}}^k - b^{k+1} + 1 + \frac{\gamma_3^k}{\mu_3} \right).$$ Similarly, the subproblem with respect to $p$ is $$\min_p \quad\quad \rho_4(p)_+ - \gamma_4^k p + \frac{\mu_4}{2}\|{\textbf{g}}^T{\textbf{v}}^k + b^{k+1} + 1 - p\|_2^2,$$ and the solution is given by $$\label{eq:ksvm_admm_psol}
    p^{k+1} = \mathcal{S}_{\frac{\rho_4}{\mu_4}}\left( {\textbf{g}}^T{\textbf{v}}^k + b^{k+1} + 1 + \frac{\gamma_4^k}{\mu_4} \right).$$ Due to lack of space in the paper, we summarize the detailed solution steps in the Appendix (see Algorithm \[alg:admm\_enk\] ADMM-ENK)

Although there appears to be ten additional parameters (six $\rho$’s and four $\mu$’s) in the ADMM method for ENK-SVM, we can usually set the $\rho$’s to the same value and do the same for the $\mu$’s. Hence, in practice, there is only one additional parameter to tune, and our computational experience in Section \[sec:exp\_ksvm\] is that the algorithm is fairly insensitive to the $\mu$’s and $\rho$’s.

IPM Phase
---------

The second phase for solving the knowledge-based SVM problem defined by (KSVM-P) follows the same steps as that described in section \[sec:ipm\]. Note that in the knowledge-based case we have decided to solve the primal problem. This decision was based on extensive numerical experiments with both the primal and dual formulation which revealed that the primal formulation is more efficient.

We found in our experiments that by introducing slack variables and transforming the above problem into a linearly equality-constrained QP, Phase 2 of HIPAD usually requires less time to solve.

HIPAD with domain knowledge incorporation
-----------------------------------------

We formally state the new two-phase algorithm for the elastic-net KSVM in Algorithm \[alg:hipad\_ksvm\].

\[alg:hipad\_ksvm\] Given ${\textbf{w}}^0, b^0, {\textbf{a}}^0, {\textbf{c}}^0, {\textbf{u}}^0,{\textbf{v}}^0,p^0,q^0,{\textbf{s}}^0\geq{\textbf{0}},{\textbf{t}}^0\geq{\textbf{0}}$. **PHASE 1: ADMM for ENK-SVM** $({\textbf{w}},b,{\textbf{u}},{\textbf{v}}) \gets $ADMM-ENK$({\textbf{w}}^0, b^0, {\textbf{a}}^0, {\textbf{c}}^0, {\textbf{u}}^0,{\textbf{v}}^0,p^0,q^0,{\textbf{s}}^0,{\textbf{t}}^0)$ **PHASE 2: IPM for KSVM** $\widetilde{{\textbf{w}}} \gets $ non-zero components of ${\textbf{w}}$ $(\widetilde{X}, \widetilde{Y}) \gets $ sub-matrices of $(X, Y)$ corresponding to the support of ${\textbf{w}}$ $\eta_u^0 \gets {\textbf{d}}^T{\textbf{u}}-b+1$ $\eta_v^0 \gets {\textbf{g}}^T{\textbf{v}}+b+1$ $({\textbf{w}},b) \gets \textrm{SVM-IPM}$$(\widetilde{X}, \widetilde{Y}, \widetilde{{\textbf{w}}},b,{\textbf{u}},{\textbf{v}},\eta_u^0,\eta_v^0$) $({\textbf{w}},b)$

Numerical results
=================

We present our numerical experience with the two main algorithms proposed in this paper: HIPAD and its knowledge-based version HIPAD-ENK. We compare their performance with their non-hybrid counterparts, i.e., ADMM-ENSVM and ADMM-ENK, which use ADMM to solve the original SVM problem. The transition condition at the end of Phase 1 is specified in (\[eq:hipad\_transition\]), with $\epsilon_{tol} = 10^{-2}$. The stopping criteria for ADMM are as follows: $\frac{|F^{k+1}-F^k|}{\max\{1,|F^k|\}} \leq \epsilon_1, \|{\textbf{a}}-({\textbf{e}}-\bar{X}{\textbf{w}}-{\textbf{y}}b)\|_2\leq \epsilon_1, \|{\textbf{c}}-{\textbf{w}}\|_2\leq \epsilon_1$ and $ \frac{\|{\textbf{w}}^{k+1}-{\textbf{w}}\|_2}{\|{\textbf{w}}^k\|_2} \leq \epsilon_2$, with $\epsilon_1 = 10^{-5}$, and $\epsilon_2 = 10^{-3}$.

HIPAD vs ADMM
-------------

To demonstrate the practical effectiveness of HIPAD, we tested the algorithm on nine real data sets which are publicly available. **rcv1** [@lewis2004rcv1] is a collection of manually categorized news wires from Reuters. The original multiple labels have been consolidated into two classes for binary classification. **real-sim** contains UseNet articles from four discussion groups, for simulated auto racing, simulated aviation, real autos, and real aviation. Both **rcv1** and **real-sim** have large feature dimensions but are highly sparse. The rest of the seven data sets are all dense data. **rcv1** and **real-sim** are subsets of the original data sets, where we randomly sampled 500 training instances and 1,000 test instances. **gisette** is a handwriting digit recognition problem from NIPS 2003 Feature Selection Challenge, and we also sub-sampled 500 instances for training. (For testing, we used the original test set of 1,000 instances.) **duke**, **leukemia**, and **colon-cancer** are data sets of gene expression profiles for breast cancer, leukemia, and colon cancer respectively. **fMRIa**, **fMRIb**, and **fMRIc** are functional MRI (fMRI) data of brain activities when the subjects are presented with pictures and text paragraphs. The data was compiled and made available by Tom Mitchell’s neuroinformatics research group [^1]. Except the three fMRI data sets, all the other data sets and their references are available at the LIBSVM website [^2].

The parameters of HIPAD, ADMM-ENSVM, and LIBSVM were selected through cross validation on the training data. We summarize the experiment results in Table \[tab:real\_data\_results\]. Clearly, HIPAD produced the best overall predication performance. In order to test the significance of the difference, we used the test statistic in [@iman1980approximations] based on Friedman’s $\chi^2_F$, and the results are significant at $\alpha=0.1$. In terms of CPU-time, HIPAD consistently outperformed ADMM-ENSVM by several times on dense data. The feature support sizes selected by HIPAD were also very competitive or even better than the ones selected by ADMM-ENSVM. In most cases, HIPAD was able to shrink the feature space to below 10$\%$ of the original size.

Simulation for Knowledge Incorporation {#sec:exp_ksvm}
--------------------------------------

We generated synthetic data to simulate the example presented at the beginning of Section \[sec:enk-svm\] in the high dimensional feature space. Specifically, four groups of multi-variate Gaussians $K_1,\cdots,K_4$ were sampled from $\mathcal{N}(\mu_1^+, \Sigma_1), \cdots, \mathcal{N}(\mu_4^+, \Sigma_4)$ and $\mathcal{N}(\mu_1^-, \Sigma_1), \cdots, \mathcal{N}(\mu_4^-, \Sigma_4)$ for four disjoint blocks of feature values (${\textbf{x}}_{K_1}, \cdots, {\textbf{x}}_{K_4}$). For positive class samples, $\mu_1^+ = 2, \mu_2^+ = 0.5, \mu_3^+ = -0.2, \mu_4^+ = -1$; for negative class samples, $\mu_1^- = -2, \mu_2^- = -0.5, \mu_3^- = 0.2, \mu_4^- = 1$. All the covariance matrices have 1 on the diagonal and 0.8 everywhere else. The training samples contain blocks $K_2$ and $K_3$, while all four blocks are present in the test samples. A random fraction (5%-10%) of the remaining entries in all the samples are generated from the standard Gaussian distribution.

The training samples are apparently hard to separate because the values of blocks $K_2$ and $K_3$ for the two classes are close to each other. However, blocks $K_1$ and $K_4$ in the test samples are well-separated. Hence, if we are given information about these two blocks as general knowledge for the entire population, we could expect the resulting SVM classifier to perform much better on the test data. Since we know the mean values of the distributions from which the entries in $K_1$ and $K_4$ are generated, we can supply the following information about the relationship between the block sample means and class membership to the KSVM: $\frac{1}{L_1}\sum_{i\in K_1}x_i \geq 4 \Rightarrow {\textbf{x}}\in A^+$, and $\frac{1}{L_4}\sum_{i\in K_4}x_i \geq 3 \Rightarrow {\textbf{x}}\in A^-$ where $L_j$ is the length of the $K_j$, $j=1,\cdots,4$, $A^+$ and $A^-$ represent the positive and negative classes, and the lowercase $x_i$ denotes the $i$-th entry of the sample ${\textbf{x}}$. Translating into the notation of (KSVM-P), we have $$\small
  B = \left(
        \begin{array}{cccccc}
          {\textbf{0}}& \underbrace{-\frac{{\textbf{e}}}{L_1}^T } & {\textbf{0}}& {\textbf{0}}& {\textbf{0}}&{\textbf{0}}\\
          & K_1 & & & & \\
        \end{array}
      \right)
   , d = -4, \label{eq:enk_data_knowledge1} \mbox{ and }
 D = \left(
        \begin{array}{cccccc}
          {\textbf{0}}& {\textbf{0}}& {\textbf{0}}& {\textbf{0}}& \underbrace{-\frac{{\textbf{e}}}{L_4}^T }  & {\textbf{0}}\\
          & & & & K_4 & \\
        \end{array}
      \right)
   , g = -3.$$

The information given here is not precise, in that we are confident that a sample should belong to the positive (or negative) class only when the corresponding block sample mean well exceeds the distribution mean. This is consistent with real-world situations, where the domain or expert knowledge tends to be conservative and often does not come in exact form.

We simulated two sets of synthetic data for ENK-SVM as described above, with $(N_{train} = 200, N_{test} = 400, m_{train} = 10,000)$ for **ksvm-s-10k** and $(N_{train} = 500, N_{test} = 1,000, m_{test} = 50,000)$ for **ksvm-s-50k**. The number of features in each of the four blocks ($K_1,K_2,K_3,K_4$) is 50 for both data sets. Clearly, HIPAD-ENK is very effective in terms of speed, feature selection, and prediction accuracy on these two data sets. Even though the features in blocks $K_1$ and $K_4$ are not discriminating in the training data, HIPAD-ENK was still able to identify all the 200 features in the four blocks correctly and exactly. This is precisely what we want to achieve as we explained at the beginning of Section \[sec:enk-svm\]. The expert knowledge not only helps rectify the separating hyperplane so that it generalizes better on the entire population, but also makes the training algorithm realize the significance of the features in blocks $K_1$ and $K_4$.

Conclusion
==========

We have proposed a two-phase hybrid optimization framework for solving the ENSVM, in which the first phase is solved by ADMM, followed by IPM in the second phase. In addition, we have proposed a knowledge-based extension of the ENSVM which can be solved by the same hybrid framework. Through a set of experiments, we demonstrated that our method has significant advantage over the existing method in terms of computation time and the resulting prediction accuracy. The algorithmic framework introduced in this paper is general enough and potentially applicable to other regularization-based classification or regression problems.

Appendix
========

\[alg:admm-ensvm\] Given ${\textbf{w}}^0, b^0, {\textbf{a}}^0, {\textbf{c}}^0, \gamma_1^0,$ and $\gamma_2^0$. $({\textbf{w}}^{k+1}, b^{k+1}) \gets$ PCG solution of: $\left(
          \begin{array}{cc}
            (\lambda_2 + \mu_2)I + \mu_1 X^TX & \mu_1 X^T{\textbf{e}}\\
            \mu_1{\textbf{e}}^TX & \mu_1 N \\
          \end{array}
        \right)
        \left(
          \begin{array}{c}
            {\textbf{w}}^{k+1} \\
            b^{k+1} \\
          \end{array}
        \right) =
        \left(
          \begin{array}{c}
            X^TY\gamma_1^k - \mu_1 X^TY({\textbf{a}}^k-{\textbf{e}}) - \gamma_2^k + \mu_2{\textbf{c}}^k \\
            {\textbf{e}}^TY\gamma_1^k - \mu_1{\textbf{e}}^TY({\textbf{a}}^k - {\textbf{e}}) \\
          \end{array}
        \right)$ \[line:admm\_linsys\]

${\textbf{a}}^{k+1} \gets \mathcal{S}_{\frac{1}{N\mu_1}}\left( {\textbf{e}}+ \frac{\gamma_1^k}{\mu_1} - Y(X{\textbf{w}}^{k+1} + b^{k+1}{\textbf{e}}) \right)$ \[line:a\_sol\]

${\textbf{c}}^{k+1} \gets \mathcal{T}_{\frac{\lambda_1}{\mu_2}}\left( \frac{\gamma_2^k}{\mu_2}+{\textbf{w}}^{k+1} \right)$ \[line:c\_sol\] $\gamma_1^{k+1} \gets \gamma_1^k + \mu_1({\textbf{e}}- Y(X{\textbf{w}}^{k+1} + b^{k+1}{\textbf{e}}) - {\textbf{a}}^{k+1})$ $\gamma_2^{k+1} \gets \gamma_2^k + \mu_2({\textbf{w}}^{k+1} - {\textbf{c}}^{k+1})$ $({\textbf{w}}^K,b^K)$

![Illustration of the early convergence (in approximately 50 iterations) of the feature support for ADMM.[]{data-label="fig:fea_supp_plots"}](fea_supp_plots.png){width="99.00000%"}

\[alg:admm\_enk\] Given ${\textbf{w}}^0, b^0, {\textbf{a}}^0, {\textbf{c}}^0, {\textbf{u}}^0,{\textbf{v}}^0,p^0,q^0,{\textbf{s}}^0\geq{\textbf{0}},{\textbf{t}}^0\geq{\textbf{0}},\gamma_i^0$, and the parameters $\lambda_1,\lambda_2,\rho_i, i = 1,\cdots,6$. $({\textbf{w}}^{k+1}, b^{k+1}) \gets $ PCG solution of the structured linear system \[line:admm\_enk\_linsys\_wb\] ${\textbf{u}}^{k+1} \gets $ the solution of the linear system ; $\>\>\>$ ${\textbf{s}}^{k+1} \gets \max({\textbf{0}}, {\textbf{u}}^{k+1} - \frac{\gamma_5^k}{\mu_5})$ ${\textbf{v}}^{k+1} \gets $ the solution of the linear system ; $\>\>\>$ ${\textbf{t}}^{k+1} \gets \max(0,{\textbf{v}}^{k+1} - \frac{\gamma_6^k}{\mu_6})$ ${\textbf{a}}^{k+1} \gets \mathcal{S}_{\frac{1}{N\mu_1}}\left( {\textbf{e}}+ \frac{\gamma_1^k}{\mu_1} - Y(X{\textbf{w}}^{k+1} + b^{k+1}{\textbf{e}}) \right)$; $\>\>\>$ ${\textbf{c}}^{k+1} \gets \mathcal{T}_{\frac{\lambda_1}{\mu_2}}\left( \frac{\gamma_2^k}{\mu_2}+{\textbf{w}}^{k+1} \right)$ $q^{k+1} \gets \mathcal{S}_{\frac{\rho_2}{\mu_3}}\left( {\textbf{d}}^T{\textbf{u}}^k - b^{k+1} + 1 + \frac{\gamma_3^k}{\mu_3} \right)$; $\>\>\>$ $p^{k+1} \gets \mathcal{S}_{\frac{\rho_4}{\mu_4}}\left( {\textbf{g}}^T{\textbf{v}}^k + b^{k+1} + 1 + \frac{\gamma_4^k}{\mu_4} \right)$ $\gamma_1^{k+1} \gets \gamma_1^k + \mu_1({\textbf{e}}- Y(X{\textbf{w}}^{k+1} + b^{k+1}{\textbf{e}}) - {\textbf{a}}^{k+1})$; $\>\>\>$ $\gamma_2^{k+1} \gets \gamma_2^k + \mu_2({\textbf{w}}^{k+1} - {\textbf{c}}^{k+1})$ $\gamma_3^{k+1} \gets \gamma_3^k + \mu_3({\textbf{d}}^T{\textbf{u}}^{k+1}-b^{k+1}+1-q^{k+1})$; $\>\>\>$ $\gamma_4^{k+1} \gets \gamma_4^k + \mu_4({\textbf{g}}^T{\textbf{v}}^{k+1}+b^{k+1}+1-p^{k+1})$ $\gamma_5^{k+1} \gets \gamma_5^k + \mu_5({\textbf{s}}^{k+1} - {\textbf{u}}^{k+1})$; $\>\>\>$ $\gamma_6^{k+1} \gets \gamma_6^k + \mu_6({\textbf{t}}^{k+1} - {\textbf{v}}^{k+1})$ $({\textbf{w}}^K,b^K)$

[^1]: http://www.cs.cmu.edu/ tom/fmri.html

[^2]: http://www.csie.ntu.edu.tw/ cjlin/libsvmtools/datasets
