Due to the advanced technology for data collection over the past decades, there has been a surge of data complexity in many research fields such as genomics, genetics, and finance, among others. Consequently, it is very common for the number of predictors in the dataset to be far larger than the number of observations (Donoho et al., 2000).
For example, in genomics it is crucial to build a classifier for the purpose of disease diagnosis, with thousands of candidate genes at hand but only tens of instances available for study. Such high dimensionality in data makes Le Zhou is with the School of Statistics, University of Minnesota, Minneapolis, MN 55455, USA (email: zhou0819@umn.edu).
Boxiang Wang is with the Department of Statistics and Actuarial Science, University of Iowa, Iowa City, IA 52242, USA (email: boxiangwang@uiowa.edu).
Yuwen Gu is with the Department of Statistics, University of Connecticut, Storrs, CT 06269, USA (email: yuwen.gu@uconn.edu).
Hui Zou is with the School of Statistics, University of Minnesota, Minneapolis, MN 55455, USA (email: zouxx019@umn.edu). Zou is supported in part by NSF grants DMS 1915842 and 2015120. traditional classification methods infeasible and poses new challenges from both theoretical and computational perspectives.
One method for performing high dimensional classification is the penalized large margin classifier. The standard support vector machine (SVM), initially proposed and investigated in Boser et al. (1992) and Vapnik (1995), has an objective equal to hinge loss plus an 2 penalty. It is also referred to as 2 -norm SVM. When the dimension greatly exceeds the sample size and there are many noisy features in the predictor set, it has been shown that it is more beneficial to use a sparse penalty such as the 1 norm penalty (a.k.a. the lasso) to replace the 2 norm penalty in order to perform classification and variable selection simultaneously in high dimensional setting (Zhu et al., 2003).
Consider the 1 norm SVM for example. It can be written as
where L(u) = (1 -u) + is the hinge loss. Just like in lasso regression, the 1 penalty induces sparsity in the solution and is thus capable of removing irrelevant features. More recently, Peng et al. (2016) investigated the rate of convergence of the 1 -norm SVM and an error bound of O( s log p n ) was established in their paper. The sparse penalized SVM can be computationally intensive especially when the number of predictors is huge in the dataset, owing to the non-differentiable loss function part. It is known that penalized problem in high dimensions with a smooth loss function can be efficiently computed by cyclical coordinate descent algorithm (Friedman et al., 2010). Nevertheless, the SVM is based on the non-differentiable hinge loss, which means that there is no convergence guarantee if one uses cyclical coordinate descent to solve the SVM. In principle, coordinate descent may not give the right solution due to the non-differentiability of the objective function (Luo and Tseng, 1992;Tseng, 2001). A similar problem under regression context is the quantile regression, in which the check loss is not differentiable (Fan et al., 2020). The typical method of solving quantile regression is the interior point algorithm. Since 1 -norm SVM can be transformed into linear programming, one may also consider interior point algorithm for solving it.
However, interior point algorithm may not scale well with high dimensional input and thus is not suitable for solving SVM in high dimensions.
Recently, Fernandes et al. (2021) studied an interesting smoothing technique for solving quantile regression with statistical guarantees. Later, Tan et al. (2021) further studied the smoothing quantile regression under high dimensional settings and showed that the statistical property of quantile regression is maintained after smoothing.
Motivated by their work, we develop a smooth version of SVM from statistical perspective, as opposed to trying to solve it exactly. Consider the first term in (I.1)
which is non-smooth. If we could replace it by some smooth loss such that the resulting estimator has nice theoretical properties, then we should focus on solving the smooth problem instead of the original problem. In fact, one may view (I.2) as the expectation of the hinge loss function with respect to the empirical measure assigning 1 n probability mass to each y i (x T i β + β 0 ), i = 1, . . . , n. The empirical measure is viewed as an estimator for the true distribution of the random variable y(x T β + β 0 ). Clearly, if we estimate the true distribution by using a smoothed kernel density estimator, then we can take the expectation of the hinge loss function with respect to the distribution determined by the smoothed kernel density estimator. This leads us to a new objective function
which we use to replace (I.2). Here h is the bandwidth of kernel density estimator and is used to index the new classifier. The resulting estimator is named as density convoluted support vector machine (DCSVM), since the kernel density estimator has a convolution interpretation. We study the following general form of penalized DCSVM in high dimensions,
The resulting estimator is called elastic-net DCSVM, which involves both 1 -DCSVM and 2 -DCSVM as special cases. By its convexity and smoothness, elastic-net DCSVM can be efficiently solved by using the generalized coordinate descent algorithm (Yang and Zou, 2013).
In this paper, we first study the theoretical properties of the elastic-net DCSVM. We show that the convergence rate of the elastic-net DCSVM is O p ( s log p n ) under the general random design setting. Furthermore, we develop novel efficient algorithm for computing elastic-net DCSVM. We use simulation studies and 8 benchmark datasets to demonstrate that elastic-net DCSVM delivers superior classification performance over its competitors, and the computational speed of DCSVM can be two orders of magnitude faster than that of SVM.
We first introduce some notation that is used throughout the paper. For an arbitrary index set A ⊂ {1, . . . , p}, any vector c = (c 1 , . . . , c p ) and any n × p matrix U, let c A = (c i , i ∈ A), and let U A be the submatrix with columns of U whose indices are in A. The complement of an index set A is denoted as A c = {1, . . . , p} \ A.
For any finite set B, let |B| be the number of elements in B. For a vector c ∈ R p and q ∈ [1, ∞), let c q = ( p j=1 |c j | q ) 1 q be its q norm, let c ∞ (or c max ) = max j |c j | be its ∞ norm, and let c min = min j |c j | be its minimum absolute value. For a matrix M, let λ min (M) and λ max (M) be its eigenvalue with smallest absolute value and largest absolute value, respectively. This is the common notation for eigenvalues of a matrix, and λ min , λ max should not be confused with the penalization parameter used in a penalty function. For any matrix G, let G = λ max (G T G) be its spectral norm. In particular, for a vector c, c = c 2 . For a, b ∈ R, let a ∧ b = min{a, b} and a ∨ b = max{a, b}. For a sequence {a n } and another nonnegative sequence {b n }, we write
For two sequences of random variables Z n and Z n , we write
, where x i = (x i1 , . . . , x ip ) T ∈ R p are predictors and y i ∈ {-1, 1} is the class label for the ith subject. We use X = (X 1 , . . . , X p ) to denote the design matrix, where X j = (x 1j , . . . , x nj ) T contains observations for the jth variable, and use y = (y 1 , . . . , y n ) T to represent the response vector. We focus on the general case where the observed data {(y i , x i )} n i=1 are i.i.d. samples from the distribution of a random vector (y, x). Let the jth component of the random vector x be denoted as x j . Meanwhile,
To perform the classification task, the support vector machine (SVM, Vapnik, 1995) seeks a separating hyperplane {x :
It is well known that the above problem can be equivalently formulated as a penalized empirical risk minimization problem:
where L(u) = (1 -u) + = max{1 -u, 0} is known as the SVM hinge loss and λ 0 > 0 is a tuning parameter that is one-to-one correspondent to the constant c in problem (II.1).
)≤t} is the empirical cdf based on i.i.d. realizations of U . The usage of the discontinuous empirical cdf here makes the objective in (II.2) to have the same degree of smoothness as the hinge loss L(•), i,e. continuous but nondifferentiable. This has motivated us to consider an alternative estimator for the cdf. If we use an estimator F (• ; β, β 0 ) that is smooth enough, the ∞ -∞ L(t)d F (t; β, β 0 ) shall lead us towards a new objective which is differentiable to certain degrees.
In particular, we consider the cdf from the kernel density estimator
where
and h > 0 is the bandwidth parameter to be tuned. Replacing F by F gives the new objective function,
where
) and " * " stands for convolution. As such, with the penalty term λ 0 β 2 2 , we obtain
We treat the classifier arisen from the above problem as a new classifier and coin it the density-convoluted SVM (DCSVM).
As discussed above, DCSVM originates from a statistical view of the SVM, while it shows merit from the computational perspective as it overcomes the non-differentiability of the original SVM problem. Smoothing a non-differentiable problem through convolution can be traced back to the idea of mollification (Friedrichs, 1944) and has also been studied in the optimization community, for example, Bertsekas (1973) and Rubinstein (1983).
The method was recently adopted to smooth the quantile regression by He et al. (2021), Fernandes et al. (2021) and Tan et al. (2021).
In this work, we focus on two most popular kernel functions, Gaussian kernel and Epanechnikov kernel in DCSVM, and we denote the corresponding convoluted loss function by L G h (v) and L E h (v), respectively. For the Gaussian kernel
where Φ(•) is the cumulative distribution function of the standard normal distribution.
For the Epanechnikov kernel K(u) = 3 4 (1 -u 2 )I(-1 ≤ u ≤ 1), where I(•) is the indicator function,
, the density-convoluted SVM loss functions with Gaussian kernel (left) and Epanechnikov kernels (right). Bottom row: plots of the first-order derivatives, L G h (v) and L E h (v).
The top row of Figure 1 depicts the DCSVM losses with Gaussian kernel and Epanechnikov kernel.
Intuitively, h should be small such that the DCSVM is very close to the SVM. According to density estimator, the optimal rate for h is
). So, we adopt h = Cn -1/5 in our implementation, where C is some numerical constant within the range (0.25, 3). The actual value of C in practice can be determined by cross-validation.
In high dimensions, we consider designing the estimator under a sparsity assumption that β * has many zero components. Let A = {j : β * j = 0, 1 ≤ j ≤ p} be the support set of β * , i.e., the set of indices of the important covariates. Let s = |A|. Throughout this paper, we allow p = pn and s = sn to diverge with n, and we assume sn ≥ 1 and pn goes to infinity as n goes to infinity. For convenience, we still use p and s to represent these quantities since no confusion is caused. In ultra-high dimensions, the dimension p is allowed to increase exponentially with the sample size n. We also assume that s is relatively of smaller order compared to n, which is necessary for the existence of a consistent estimator.
To perform the classification for high-dimensional data, we present sparse DCSVM with an additional 1-penalty term ( β0, β) := argmin
The 1-penalty is used to induce sparsity in the estimator. We also consider the following version of sparse DCSVM with only an 1-penalty term:
( β0, β) := argmin
Borrowing the commonly used terminologies for different penalties in high dimensional literature, we refer to the estimator in (II.4) as elastic-net DCSVM, and refer to the estimator in (II.5) as lasso DCSVM. Note that the lasso DCSVM is a special case of elastic-net DCSVM with λ0 = 0. In the statistics literature it is well-known than the elastic-net often yields much improved prediction than the lasso owing to the additional 2 regularization that can well handle the correlated covariates which occurs often in high-dimensional data. For example, in the case of the SVM Wang et al. (2006) showed that the elastic-net regularized SVM is more accurate than the 1-norm SVM. Therefore, we focus on the elastic-net penalized DCSVM in theory and in applications.
We now state the assumptions needed to establish our theoretical results. We first impose the following conditions on the random design.
Assumption 1. {(yi, xi)} n i=1 , (y, x) are independent and identically distributed on R × R p . x is a zero-mean sub-exponential random vector, i.e. E[x] = 0, and there exists a constant m0 > 0 such that sup
By definition of Orlicz norm and Markov's inequality, this further implies sup
For any index set A ⊂ {1, . . . , p}, consider the cone
be Hessian matrix of the population loss, or information matrix. We impose the following condition on the information.
Assumption 2. There exists a constant ρ > 0 such that
Assumption 1 is a general setting in the random design, which relaxes the classical condition that the components of x are bounded random variables (Peng et al., 2016). Assumption 2, which is a restricted eigenvalue type of condition, is needed to establish 2-type error bound for 1-penalized type of estimator. Similar conditions have been widely adopted in the literature (Bühlmann and Van De Geer, 2011;Fan et al., 2020).
Theorem 1. Assume assumptions 1-2 hold, and s log p = o(n). Choose the tuning parameters such that 8λ0 β * max < λ. Then there exists a large enough constant c0 > 0 such that with the choice λ = c0 log p n , the elastic-net DCSVM estimator ( β0, β) satisfies
Theorem 1 shows that the sparse DCSVM estimator can achieve the same rate of convergence as the 1-SVM (Peng et al., 2016). Meanwhile, the sparse DCSVM has better computational efficiency than penalized SVM due to the smoothness of its loss function, as shown in the next section.
In this section, we develop an efficient algorithm for computing the solution path of DCSVM.
At the outset, we present the first-order derivative of the density-convoluted SVM loss and show they are Lipschitz continuous in Lemma 1:
and L E h (v) be the DCSVM loss using Gaussian kernel and Epanechnikov kernel, respectively. For
where the Lipschitz constants are given as
and c E h = 3 4h .
The bottom row of Figure 1 depicts L G h (v) and L E h (v). Lemma 1 gives rise to the following quadratic majorization condition for the DCSVM:
where L h is exemplified by L G h and L E h and c h is the corresponding Lipschitz constant. Based on the Lipschitz condition, we develop a generalized coordinate descent (GCD) algorithm (Yang and Zou, 2013) to solve those sparse penalized DCSVMs. We first consider the adaptive lasso penalty. The algorithm can be easily adjusted to handle lasso and elastic net.
Without loss of generality, we assume each X j has zero mean and unit length. In a coordinate-wise manner, suppose the coordinate β 1 , β 2 , . . . , β j-1 have been updated and we now update β j . Denote by β0 and β by the current solution and let v i = y i ( β0 + x T i β). To update β j , instead of solving the coordinate-wise update function,
we solve its majorization function
Likewise, β 0 is updated to be β0 -1
In the appendix we provide theoretical analysis of the convergence of the generalized coordinate descent algorithm which is not in Yang and Zou (2013). In our implementation, we further apply the strong rule (Tibshirani et al., 2010), warm start, and active set strategy (Friedman et al., 2010) to further accelerate the algorithm.
In this section, we use several simulation examples to demonstrate the performance of DCSVM.
The response variables of all the simulated data are binary and the two classes are balanced, i.e., P (Y = 1) = P (Y = -1) = 0.5. In each example, define the p-dimensional mean vectors µ + = (0.7, 0.7, 0.7, 0.7, 0.7, 0, 0, . . . , 0) and µ -= -µ + , where p = 500 or 5000 in our experiments. Each observation from the positive class is drawn from N(µ + , Σ) and each observation from the negative class is drawn from N(µ -, Σ). We consider three different choices of Σ. In example 1, Σ = I p×p so the variables are independent. In both examples 2 and 3,
where Σ 5×5 have all diagonal elements of 1 and off-diagonal elements of ρ in example 2, and (Σ 5×5 ) i,j = ρ |i-j| in example 3. We use ρ = 0.2, 0.7, and 0.9.
We first compared elastic-net DCSVM with Gaussian kernel and Epanechnikov kernel with elastic-net SVM (Wang et al., 2006) and elasticnet logistic regression that is fitted using the R package gcdnet (Yang and Zou, 2013). For each example, the training size is 200 and we use five-fold cross-validation to select the best tuple of (h, λ 0 , λ) where h is chosen from 0.1, 0.25, 0.5, and 1, λ 0 is selected from 0.5 * (10 -4 , 10 -3 , 10 -2 , 10 -1 , 1, 5), and λ is searched along the solution path; for the SVM and logistic regression, we select λ 0 and λ in the same manner.
We record the prediction error and run time in Table I. The run time include all the time spent on tuning and training the models. We observe the DCSVM with Epanechnikov kernel has slightly better performance than DCSVM with Gaussian kernel, and both of them have better prediction accuracy than the other two methods. DCSVM with Epanechnikov kernel is the fastest while the elastic-net SVM is the slowest.
All the methods exhibited in Table I use elastic-net penalty. We now study the performance when using other sparse penalities. Due to the overall best performance, we stay with DCSVM with Epanechnikov kernel and we compare the prediction accuracy and variable selection when using lasso and elastic-net penalties. We present the results in Table II. In general, we find the elastic-net has the best performance in both prediction error and variable selection.
In this section, we demonstrate the performance of DCSVM using several benchmark data, which are available from UCI machine learning repository. We randomly split each data set into a training set and a test set with a 1:1 ratio. On the training set, we fit elastic-net DCSVM, elastic-net logistic regression, and elastic-net SVM, and tune each method using five-fold cross-validation. The prediction accuracy is computed based on the test set.
We present the result in Table III. We observe the elastic-net DCSVM has the best performance in general.
We first give some general formula regarding the loss function L h and its derivatives.
Proof of Theorem 1. By definition of the 1 -penalized CRR estimator and triangle inequality, we have
where we denote u := ββ * . On the other hand, by convexity of L h (•), we have
So by Hoeffding's inequality,
Meanwhile, we have E L h y(x T β * + β * 0 ) yx = 0 by the definition of β * and optimality condition. By the choice of tuning parameters we have
Notice that by assumption 1 and
. . , n}, ∀j ∈ {1, . . . , p}. By Theorem 1.4 in Götze et al. (2021), there exists an absolute constant η 0 > 0 such that
)n .
So following (A.0.5) we have
)n . Now, under E 1 ∩ E 2 , combining (A.0.2) and (A.0.3) we have
We give an upper bound for E[H(r)]. Let σ 1 , . . . , σn be i.i.d. Rademacher random variables (i.e. P(σ i = 1) = P(σ i = -1) = 1 2 ), which is independent from all the other random elements. By the symmetrization inequality (see for instance, Lemma 2.3.1 in Van Der Vaart and Wellner (1996)) and contraction inequality (see for instance, Theorem 4.12 in Ledoux and Talagrand (1991)), |L h (•)| ≤ 1 and Cauchy-Schwarz inequality, we have
By assumption 1 and definition of Orlicz norm, we know σ i y i x ij ψ 1 = x ij ψ 1 ≤ m 0 , ∀i ∈ {1, . . . , n}, ∀j ∈ {1, . . . , p}. Also, it is straightforward to see σ i y i ψ 1 = 1 log 2 . By Proposition 2.7.1 in Vershynin (2018), there exists a constant c 1 > 0 such that E[e tσ i y i
, ∀i ∈ {1, . . . , n}, ∀j ∈ {1, . . . , p}. By Jensen's inequality, we have for any
Consequently, for any
By the condition of Theorem 1, we know 1), so for large enough n,
for large enough n. Thus, combining (A.0.6) and (A.0.8) we get
This implies that H(r) = Op( rs log p n
). Define event G T := {H(r) ≤ T rs log p n } for any T > 0, then we have lim T →∞ lim sup n→∞ P (G c T ) = 0.
Next, for any (β 0 , β) ∈ (β * 0 , β * )+C(r), we derive a lower bound for E[G(β 0 , β)]. For large enough n, for any (β 0 , β) ∈ (β * 0 , β * )+C(r), by Taylor's theorem and assumption 2, there exists a ∈ [0, 1] such that
(A.0.9)
On the other hand, by our choice for tuning parameters, for any (β 0 , β) ∈ (β * 0 , β * ) + C(r) we have
Thus, combining (A.0.9), (A.0.10) and (A.0.11), under G T , we have for any
. By convexity of F (•) and norm functions and by (A.0.12), we further have
, which is a contradiction with the definition of ( β0 , β). So the claim is proved. By union bound, previous results and choice of tuning parameters, we have
Since log p n = o(1), as long as c 0 is large enough (for instance c 0 > 4
Combining this result and the previous claim, the proof of Theorem 1 is finished.
It is seen that
Thus inequality (IV.1) is obtained due to the mean value theorem.
We then prove inequality (IV.2). The inequality is trivial when
since L E h is twice differentiable between 1 -h and 1 + h, we see
where the second to the last inequality is due to
where the second to the last inequality is due to
denote the subvector of v with its kth component removed by
We also let ∂h be the sub-differential of a nonsmooth convex function h (see e.g., Bertsekas, 1999).
b) Iteration Complexity Analysis: Without loss of generality, we focus solely on the GCD algorithm for solving the weighted lasso penalized DCSVM
where w k ≥ 0 are the weights of the penalty. Indeed, this formulation covers all the sparsity patterns in Section II-C. Also, the intercept term β 0 can be absorbed into the formulation by setting x i1 = 1 for i = 1, . . . , n and w 1 = 0. For ease of exposition, let us rewrite (C.0.1) as the following unconstrained optimization problem
which implies that the gradient of g(•) is uniformly Lipschitz continuous with Lipschitz constant L = c h ρmax. When restricted to each coordinate, we have
which implies that the gradient of g(•) is coordinate-wise uniformly Lipschitz continuous with Lipschitz constants L k = c h X k 2 , k = 1, . . . , p.
In the GCD (cyclic coordinate descent) algorithm, let β r be the update of β after the rth cycle, r ≥ 0. For ease of notation, denote
is equivalent to
where the proximity operator prox does the soft-thresholding (Parikh and Boyd, 2013) and
Our analysis will be divided into three parts: the sufficient descent step, the cost-to-go estimate step, and the local error bound step. Similar techniques can be found in Luo and Tseng (1992), Luo andTseng (1993), Zhang et al. (2013) and Hong et al. (2013).
c) Sufficient Descent: Consider the proximal gradient method applied to solving the following problem
we have by (C.0.3)
d) Cost-to-go Estimate: Let X * := {β * |f (β * ) = min β f (β)} be the optimal solution set of problem (C.0.2). Let βr ∈ X * be the point in X * such that d X * (β r ) := min β∈X * ββ r = βrβ r . By optimality of
By the mean value theorem, there exists λ ∈ [0, 1] and ξ r = λβ r+1 + (1 -λ) βr such that g(β r+1 ) -g( βr ) = ∇g(ξ r ), β r+1 -βr .
It follows that
By the fact that ∇g(•) is Lipschitz continuous, it is implied that
It follows that
Here we handle the Gaussian and Epanechnikov kernels separately.
For the Gaussian kernel, that is, when L h (•) = L G h (•), according to (C.0.4) and (C.0.5), the GCD algorithm is descending along its iterations. We can thus restrict the domain of β to the sublevel set
We can see that g(β) = p(Xβ). It follows from Zhang et al. (2013) that for any ξ ≥ min β f (β), there exist κ, ε > 0 such that
for all β such that βprox h (β -∇g(β)) ≤ ε and f (β) ≤ ξ.
For the Epanechnikov kernel, that is, when L h (•) = L E h (•), one needs to add an additional ridge penalty µ β 2 for some small µ > 0 in order to achieve strong optimality. Thus, when the Epanechnikov kernel is used, we instead consider the following problem
and solve it using the GCD algorithm.
As a summary, we show in the following theorem that the GCD algorithm converges at least linearly.
Theorem 2. The GCD algorithm converges at least linearly to a solution in X * .
Proof. We first show that there exists some σ > 0 such that
For any r ≥ 1 and any 1 ≤ k ≤ p, by the optimality of
we have
It follows that
where L = max(1, L) and L = max(1, L -1 ). Therefore, when we take σ we obtain
which implies that ∆ r+1 ≤ c 3 1 + c 3 ∆ r , (C.0.9)
where c 3 = (c 2 κ 2 σ 2 + c 2 )c -1 1 . We can see from (C.0.9) that f (β r ) approaches f * with at least linear rate of convergence. From (C.0.5) again, this further implies that the sequence {β r } converges at least linearly.
Under the same settings introduced in our simulation section, we compared the performance of lasso DCSVM and elastic-net DCSVM, using Gaussian kernel. The result is shown in Table S.1. Again, we can see that the elastic-net DCSVM outperforms lasso DCSVM. We also conducted elastic-net DCSVM with Gaussian kernel on the same real datasets that we introduced in our real data section, and compared its performance with the performance of elastic-net SVM and elastic-net logistic regression. The result is displayed in Table S