An Accelerated Proximal Gradient Algorithm for Singly Linearly Constrained Quadratic Programs with Box Constraints

Abstract

Recently, the existed proximal gradient algorithms had been used to solve non-smooth convex optimization problems. As a special nonsmooth convex problem, the singly linearly constrained quadratic programs with box constraints appear in a wide range of applications. Hence, we propose an accelerated proximal gradient algorithm for singly linearly constrained quadratic programs with box constraints. At each iteration, the subproblem whose Hessian matrix is diagonal and positive definite is an easy model which can be solved efficiently via searching a root of a piecewise linear function. It is proved that the new algorithm can terminate at an 𝜀-optimal solution within 𝑂(1/√𝜀) iterations. Moreover, no line search is needed in this algorithm, and the global convergence can be proved under mild conditions. Numerical results are reported for solving quadratic programs arising from the training of support vector machines, which show that the new algorithm is efficient.

1. Introduction

In this paper, we mainly consider the following quadratic programming problem: min_𝑥∈Ω𝑓(𝑥)=12𝑥^𝑇𝐺𝑥−𝑏^𝑇𝑥, (1) where 𝐺∈𝑅^𝑛×𝑛 is symmetric and positive semidefinite, 𝑥,𝑏∈𝑅^𝑛, and the feasible region Ω is defined by Ω=𝑥∈𝑅^𝑛,𝑙≤𝑥≤𝑢,𝑐^𝑇𝑥=𝑑, (2) where 𝑙,𝑢,𝑐∈𝑅^𝑛 and 𝑑 is a scalar.

The singly linearly constrained quadratic programs with box constraints appear in a wide range of applications such as image processing, biological information and machine learning. Specifically, support vector machine (SVM) is one of the most classical models of (1). It is a promising technique for solving a variety of machine learning and function estimation problems. The SVM learning methodology has been shown to give good performance in a wide variety of problems such as face detection, text categorization, and handwritten character recognition. The number of variables in SVM is so huge that traditional optimization methods cannot be directly applied. Some decomposition method [1–4] with its subproblem being a special case of (1) is the main approach for large-scale SVM problems. The solution of the subproblem by generalized variable projection method (GVPM) and projected gradient method (PGM) is introduced in [5, 6], respectively. Moreover, Zanni, and so forth, also proposes the parallel decomposition algorithms based on these two methods in [7, 8]. For more general large-scale model, some of the parallelization pieces of literature in [9, 10] recently are proposed, but these methods cannot be applied specifically to SVM so far.

In this work, we will give an accelerated proximal gradient algorithm for (1). First, we consider the following nonsmooth convex optimization problem: min_{𝑥∈𝑅𝑛}𝐹(𝑥)=ℎ(𝑥)+𝑔(𝑥), (3) where 𝑔∶𝑅^𝑛→(−∞,∞] is a proper, lower semicontinuous (lsc), convex function and ℎ is convex smooth (i.e., continuously differentiable) on an open subset of 𝑅^𝑛 containing dom 𝑔={𝑥∣𝑔(𝑥)<∞}.∇ℎ is Lipschitz continuous on dom 𝑔. That is, ∇ℎ(𝑥)−∇ℎ𝑦≤𝐿𝑥−𝑦∀𝑥,𝑦∈dom𝑔 (4) for some 𝐿>0, where, and in what follows, ‖⋅‖ denotes the spectral norm. Obviously, problem (1) is a special case of (3) with ℎ(𝑥)=𝑓(𝑥) and 𝑔(𝑥) being the indicator function for the feasible region Ω defined by 𝑔(𝑥)=0if𝑥∈Ω,+∞else. (5)

Recently, great attention has been paid to the solution of (3). Nesterov and Nemirovski [11, 12] study the accelerated proximal gradient method for (1) with an attractive iteration complexity of 𝑂(1/√𝜀) for achieving 𝜀-optimality. Almost at the same time, Beck and Teboulle [13] give a fast iterative shrinkage-thresholding algorithm (FISTA) which achieves the same convergence rate. After that, Tseng [14] summarizes these algorithms and presents a unified treatment of these methods. All these algorithms have a good performance on large-scale problems, such as linear inverse problems, matrix game problems, and matrix completion. Motivated by the successful use of the accelerated proximal gradient method for (3), we extend Beck and Teboulle's algorithm to solve (1). In particular, the subproblem is solved by searching a root of a piecewise linear continuous function. Numerical results show that the new algorithm is efficient.

The paper is organized as follows. In Section 2, the proximal gradient algorithm and its accelerated version are presented for (1). The convergence of these algorithms is also discussed. In Section 3, the method to efficiently solve the subproblem is introduced. Numerical results on SVM classification problems generated by the random and real world data sets are shown in Section 4. Section 5 is devoted to some conclusions and further study.

2. Proximal Gradient Algorithms

In this section, we introduce the proximal gradient algorithm and its accelerated version which can be applied to solve (1). For any 𝑧∈𝑅^𝑛, consider the following quadratic approximation of 𝑓(𝑥) at 𝑧: 𝑄_𝐿()(𝑧)+(𝑧+𝐿𝑥,𝑧=𝑓𝑥−𝑧,∇𝑓)2‖𝑥−𝑧‖²=𝐿2𝑥^𝑇𝑥+(𝐺𝑧−𝐿𝑧−𝑏)^𝑇𝑥+const, (6)(7) where 𝐿>0 is a given parameter, const=(𝐿/2)𝑧^𝑇𝑧−∇𝑓(𝑧)^𝑇𝑧+𝑓(𝑧). Since (6) is a strongly convex function of 𝑥, 𝑄_𝐿(𝑥,𝑧) has a unique minimizer in Ω which we denote by 𝑃_𝐿(𝑧)=argmin𝑄_𝐿(𝑥,𝑧)∶𝑥∈Ω. (8)

First, we present a simple result which characterizes the optimality of 𝑃_𝐿(𝑧).

Lemma 1. For any 𝑧∈𝑅^𝑛, 𝑤=𝑃_𝐿(𝑧) if and only if there exists 𝛾(𝑧)∈𝑁_Ω(𝑤), such that ∇𝑓(𝑧)+𝐿(𝑤−𝑧)+𝛾(𝑧)=0, (9) where 𝑁_Ω(𝑤) is the normal cone of Ω at the point 𝑤 denoted by 𝑁_Ω(𝑤)=𝑤^∗∈𝑅^𝑛∣𝑤^∗,𝑣−𝑤≤0∀𝑣∈Ω. (10)

Proof. This result can be obtained directly from the optimality conditions of (8), which is a strongly convex problem.

Then, similarly to Lemma 2.3 in [13], we have the following key result.

Lemma 2. Let 𝑧∈𝑅^𝑛 and 𝐿>0 such that 𝑓𝑃_𝐿(𝑧)≤𝑄𝑝_𝐿(𝑧),𝑧. (11) Then for any 𝑥∈Ω, 𝑓(𝑥)−𝑓𝑃_𝐿(𝑧)≥𝐿2𝑃_𝐿(𝑧)−𝑧²+𝐿𝑧−𝑥,𝑃_𝐿(𝑧)−𝑧. (12)

Proof. Since 𝛾(𝑧)∈𝑁_Ω(𝑤) and 𝑓 is convex, it follows that 𝑓(𝑥)≥𝑓(𝑧)+𝑥−𝑧,∇𝑓(𝑧)+𝛾(𝑧),𝑥−𝑃_𝐿(𝑧). (13) Hence, from (9), (11), (13), and the definition of 𝑃_𝐿(𝑧) we have 𝑓(𝑥)−𝑓𝑃_𝐿(𝑧)≥𝑓(𝑥)−𝑄𝑝_𝐿(𝑧)𝐿,𝑧≥−2𝑃_𝐿(𝑧)−𝑧²+𝑥−𝑃_𝐿(𝑧)(𝑧)(𝑧𝐿,∇𝑓+𝛾)=−2𝑃_𝐿(𝑧)−𝑧²+𝐿𝑥−𝑃_𝐿(𝑧),𝑧−𝑃_𝐿(𝑧=𝐿)2𝑃_𝐿(𝑧)−𝑧²+𝐿𝑧−𝑥,𝑃_𝐿(𝑧)−𝑧. (14) This completes the proof.

Note that 𝑓(𝑥) is a convex quadratic function and (11) can always be satisfied for 𝑃_𝐿(𝑧) if 𝐿≥𝜆_max(𝐺), where 𝜆_max denotes the maximum eigenvalue of 𝐺. The proximal gradient algorithm for (1) is as follows.

Algorithm 3. We have the following steps.
Step 1. Choose 𝑥₀∈𝑅^𝑛.
Step 2. While (𝑥_𝑘−1 does not satisfy the terminal conditions) 𝑥_𝑘=𝑃_𝐿𝑥_𝑘−1𝑘=𝑘+1. (15)

The following theorem is about the iteration complexity of Algorithm 3. In what follows, 𝑋^∗ and 𝑓^∗ denote the set of optimal solutions and optimal value of (1), respectively.

Theorem 4. Let {𝑥_𝑘} and {𝑓(𝑥_𝑘)} be the sequences generated by Algorithm 3 with 𝐿≥𝜆_max(𝐺); then, for 𝑘≥1, the following results hold: (a){𝑓(𝑥_𝑘)} is nonincreasing and 𝑓𝑥_𝑘−𝑓𝑥_𝑘+1≥𝐿2𝑥_𝑘−𝑥_𝑘+1²; (16) (b)𝑓(𝑥_𝑘)−𝑓(𝑥^∗)≤𝐿‖𝑥₀−𝑥^∗‖²/2𝑘, for all 𝑥^∗∈𝑋^∗;(c) if {𝑥_𝑘^𝑗} is a convergent subsequence of {𝑥_𝑘} and lim_𝑗→∞𝑥_𝑘^𝑗=𝑥^∗∗ then 𝑥^∗∗∈𝑋^∗;(d) if 𝐺 is positive definite, then 𝑥_𝑘−𝑥^∗≤𝐿2Λ_min(𝐺)𝑘𝑥₀−𝑥^∗. (17)

Proof. (a) Let 𝑧=𝑥_𝑘, 𝑥_𝑘+1=𝑃_𝐿(𝑧); then from (11) and the optimality of 𝑃_𝐿(𝑧) we have 𝑓𝑥_𝑘+1≤𝑄𝑥_𝑘+1,𝑥_𝑘≤𝑄𝑥_𝑘,𝑥_𝑘=𝑓𝑥_𝑘. (18)
Furthermore, take 𝑥=𝑥_𝑘 in (12); we get 𝑓𝑥_𝑘−𝑓𝑥_𝑘+1≥𝐿2𝑥_𝑘−𝑥_𝑘+1². (19)
(b) We have 𝑘𝑓𝑥_𝑘−𝑓𝑥^∗=𝑘𝑓𝑥_𝑘−^𝑘−1_𝑛=0𝑓𝑥_𝑛+1+^𝑘−1_𝑛=0𝑓𝑥_𝑛+1−𝑘𝑓𝑥^∗≤^𝑘−1_𝑛=0𝑓𝑥_𝑛+1−𝑓𝑥^∗≤−^𝑘−1_𝑛=0𝐿2𝑥_𝑛+1−𝑥_𝑛²−𝐿𝑥_𝑛−𝑥^∗,𝑥_𝑛+1−𝑥_𝑛𝐿=−2^𝑘−1_𝑛=0𝑥^∗−𝑥_𝑛+1²−𝑥^∗−𝑥_𝑛²𝐿=−2𝑥^∗−𝑥_𝑘²−𝑥₀−𝑥^∗2≤𝐿2𝑥₀−𝑥^∗2, (20) where the first inequality uses the fact that {𝑓(𝑥_𝑘)} is nonincreasing and the second inequality is established by taking 𝑥=𝑥^∗, 𝑧=𝑥_𝑛 in (12). Dividing both sides by 𝑘 in the last inequality, we get 𝑓𝑥_𝑘−𝑓𝑥^∗≤𝐿𝑥₀−𝑥^∗2.2𝑘 (21)
(c) It follows from (b) that, for all 𝑥^∗∈𝑋^∗, 𝑓𝑥_𝑘^𝑗−𝑓𝑥^∗≤𝐿𝑥₀−𝑥^∗22𝑘^𝑗. (22)
Let 𝑗→∞ in the above inequality; we have 𝑓(𝑥^∗∗)≤𝑓(𝑥^∗). On the other hand, by the definition of 𝑥^∗, we have 𝑓(𝑥^∗∗)≥𝑓(𝑥^∗). Hence, 𝑓(𝑥^∗∗)=𝑓(𝑥^∗), which implies that 𝑥^∗∗∈𝑋^∗.
(d) It follows from the definition of 𝑓(𝑥) that, for all 𝑥^∗∈𝑋^∗, 𝑓𝑥_𝑘−𝑓𝑥^∗=∇𝑓𝑥^∗𝑇𝑥_𝑘−𝑥^∗+𝑥_𝑘−𝑥^∗𝑇𝐺𝑥_𝑘−𝑥^∗≥𝑥_𝑘−𝑥^∗𝑇𝐺𝑥_𝑘−𝑥^∗≥𝜆_min(𝐺)𝑥_𝑘−𝑥^∗2, (23) where the first inequality is established by the definition of 𝑥^∗ and the second inequality, with 𝜆_min(𝐺) denoting the minimum eigenvalue of 𝐺, is obtained from the positive definite property of 𝐺.
From (23) and (b), we get 𝑥_𝑘−𝑥^∗≤𝐿2𝜆_min(𝐺)𝑘𝑥₀−𝑥^∗. (24)

Theorem 4 shows that Algorithm 3 has a result complexity of 𝑂(1/𝑘). Now we introduce a new algorithm with an improved complexity result of 𝑂(1/𝑘²).

Algorithm 5. We have the following steps.
Step 1. Choose 𝑧₁=𝑥₀∈𝑅^𝑛, set 𝑘=1, 𝑡₁=1.
Step2. While (𝑥_𝑘−1 does not satisfy the terminal conditions) 𝑥_𝑘=𝑃_𝐿𝑧_𝑘,𝑡_𝑘+1=1+1+4𝑡²_𝑘2,𝑧_𝑘+1=𝑥_𝑘+𝑡_𝑘−1𝑡_𝑘+1𝑥_𝑘−𝑥_𝑘−1,𝑘=𝑘+1. (25)(26)(27)

Remark 6. In Algorithm 5, the operator 𝑃_𝐿(⋅) is employed on the point 𝑧_𝑘 which is a specific linear combination of the previous two points {𝑥_𝑘−1,𝑥_𝑘−2}. However, in Algorithm 3, the operator only uses the previous point 𝑥_𝑘−1. On the other hand, the computational effort of these two algorithms is almost the same except for the computation of (26) in Algorithm 5 which is negligible.
Now we give the promising improved complexity result for Algorithm 5.

Theorem 7. Let {𝑥_𝑘} and {𝑓(𝑥_𝑘)} be the sequences generated by Algorithm 5 with 𝐿≥𝜆_max(𝐺). Then, for 𝑘≥1, the following results hold: (a)𝑓(𝑥_𝑘)−𝑓(𝑥^∗)≤2𝐿‖𝑥₀−𝑥^∗‖²/(𝑘+1)², for all 𝑥^∗∈𝑋^∗;(b) if {𝑥_𝑘^𝑗} is a converging subsequence of {𝑥_𝑘} and lim_𝑗→∞𝑥_𝑘^𝑗=𝑥^∗∗, then 𝑥^∗∗∈𝑋^∗;(c) if G is positive definite, then 𝑥_𝑘−𝑥^∗≤𝐿2𝜆_min(𝐺)(𝑘+1)²𝑥₀−𝑥^∗. (28)

Proof. (a) Using Lemma 2, the proof is similar to the proof of Theorem 4.4 in [13]. (b) and (c) can be proved in the same way as that of (c) and (d) in Theorem 4, respectively.

3. How to Solve 𝑃_𝐿(⋅)

In this section, we extend the method studied by Dai and Flecher [6] to solve 𝑃_𝐿(⋅). In order to obtain 𝑃_𝐿(⋅) in (8), we need to solve the following problem: 𝐿min2𝑥^𝑇𝑥−𝑞^𝑇𝑥s.t.𝑐^𝑇𝑥=𝑑,𝑙≤𝑥≤𝑢, (29) where 𝑞=𝐿𝑧+𝑏−𝐺𝑧. For a given value of 𝜆, we consider the following box constrained QP problem: 𝐿min2𝑥^𝑇𝑥−𝑞^𝑇𝑥−𝜆𝑐^𝑇𝑥−𝑑s.t.𝑙≤𝑥≤𝑢 (30) and denote the minimizer of (30) as 𝑥(𝜆). Then, 𝑥(𝜆)=mid(𝑙,𝛼,𝑢), (31) where 𝛼=𝛼(𝜆) has components 𝛼_𝑖=(𝑞_𝑖+𝜆𝑐_𝑖)/𝐿 (𝑖=1,2,…𝑛) and mid(𝑙,𝛼,𝑢) is the componentwise operation that supplies the median of its three arguments. That is, [𝑥(𝜆)]_𝑖=𝑙_𝑖+𝛼_𝑖−𝑢_𝑖−max𝑙_𝑖,𝛼_𝑖,𝑢_𝑖−min𝑙_𝑖,𝛼_𝑖,𝑢_𝑖(𝑖=1,2,…𝑛). (32)

Define 𝑟(𝜆)=𝑐^𝑇𝑥(𝜆)−𝑑; (33) then from [6], we know that 𝑟(𝜆) is a piecewise linear continuous and monotonically increasing function of 𝜆. Furthermore, 𝑥(𝜆^∗) is the optimal of (29) if 𝜆^∗ is located such that 𝑟(𝜆^∗)=0. Hence, the main task of solving (29) is to find a 𝜆^∗ such that 𝑟(𝜆^∗)=0. For this purpose, we adopt the algorithm introduced in [6] which consists of a bracketing phase and a secant phase. Numerical experiments in the next section show that this algorithm achieves high efficiency.

4. Numerical Experiments

As we know, one of the important applications of SVM is classification. Given a training set 𝐷=𝑧_𝑖,𝑦_𝑖,𝑖=1,…,𝑛,𝑧_𝑖∈𝑅^𝑚,𝑦_𝑖∈{−1,1}, (34) we need to find a hyperplane to separate the two classes of points. This can be done by solving the following convex quadratic programming problem: 1min2𝑥^𝑇𝑄𝑥−𝑒^𝑇𝑥s.t.𝑦^𝑇𝑥=0,0≤𝑥≤𝐶𝑒, (35) where 𝑦=[𝑦₁,𝑦₂,…,𝑦_𝑛]^𝑇, 𝑒∈𝑅^𝑛 is the vector of all ones, and 𝐶 is a positive scalar. 𝑄 is a 𝑛×𝑛 symmetric and positive semidefinite matrix with entries 𝑄_𝑖𝑗=𝑦_𝑖𝑦_𝑗𝐾(𝑥_𝑖,𝑥_𝑗), 𝑖,𝑗=1,2,…,𝑛, where 𝐾∶𝑅^𝑚×𝑅^𝑚→𝑅 is some kernel function.

In this section, we illustrate the performance of Algorithms 3 and 5 with Matlab 7.10 on a Windows XP professional computer (Pentium R, 2.79 GHZ). We conduct a test on SVM classification problems with the random data sets and the real world data sets.

The generation of the random data sets is based on four parameters 𝑛, 𝑚, 𝛽, and 𝜂, where 𝑛 is the number of samples and 𝑚 is the dimension of 𝑧_𝑖(𝑖=1,…,𝑛). Each element of 𝑧_𝑖 is randomly generated in [𝛽,𝜂](𝑖=1,…,𝑛), and −1 or 1 randomly emerges in the 𝑖th entry of 𝑦(𝑖=1,…,𝑛).

We have generated three random data sets with 𝛽=−1, 𝜂=1, 𝑚=5, and 𝑛 = 200, 600, and 1000, respectively. The two real world data sets are the UCI adult data set and the heart data set. For the UCI adult, the versions with 𝑛 = 1605, 2265, and 3185, are considered. The Gauss radial basis function 𝐾𝑧_𝑖,𝑧_𝑗−𝑧=exp^𝑖−𝑧_𝑗²𝜎² (36) is used in our tests. The parameters 𝐶 in (35) and 𝜎 in (36) are set to (10,√40), (1,√13), and (1,1000) for random data sets, heart, and UCI adult data sets respectively.

For all test problems, the initial point is chosen as 𝑥₀=(0,…0)^𝑇 and the const 𝐿=∑^𝑛_𝑖=1𝑄_𝑖𝑖. The terminal rules for solving 𝑃_𝐿(⋅) are |𝑟(𝜆)|≤10⁻⁸.

For the test problem with random data set (𝑛=200) and heart data set, the running steps are 8000 and 1000, respectively. In order to make a visual comparison between Algorithms 3 and 5, we plot the values of 𝑓(𝑥_𝑘) with iterations in Figure 1. Furthermore, we compute 𝐶𝑅𝑗=^𝑀_𝑘=1𝑐𝑟𝑀×𝑗−1+𝑘𝑗=1,2,…𝐽, (37) where 𝑐𝑟(𝑀×(𝑗−1)+𝑘) denotes the number of calculating 𝑟(𝜆) in step 𝑀×(𝑗−1)+𝑘. In Figure 2, we plot the values of 𝐶𝑅(𝑗)(𝑗=1,2,…𝐽) with 𝑀=100, 𝐽=80 and 𝑀=20, 𝐽=50, respectively.

(a) Random data set

(b) Heart data set

Figure 1 

Function value curve with iterations.

(a) Random data set

(b) Heart data set

Figure 2 

Number of calculating 𝑟(𝜆) for every finite step.

For the test problems with the rest data sets, the terminal condition is based on the fulfilment of the KKT conditions within 0.001 (see [15]). The numerical results are shown in Tables 1 and 2, where the meaning of the indexes are as follows: (i)𝑛: the scale of the data sets,(ii)sec: the computing time (in seconds),(iii)𝑘: the number of iterations,(iv)raver: the average number of calculating 𝑟(𝜆) in each step.

Table 1 

Numerical results for Algorithm 3.

Table 2 

Numerical results for Algorithm 5.

From Figure 2, we can see that although the total number of calculating 𝑟(𝜆) every 100 (27) steps in Algorithm 3 is less than that the number in Algorithm 5. However, Algorithm 5 uses much less iterations and time to achieve the approximate solution, compared with Algorithm 3 from Figure 1 and Tables 1 and 2. For the subproblem, less than 4 iterations are used to achieve a very accurate solution on average, which shows that 𝑃_𝐿(⋅) can be solved efficiently. Moreover, from the column raver in Tables 1 and 2, we can also see that the average number of calculating 𝑟(𝜆) with real world data sets is slightly less than that with random data sets.

5. Conclusion and Future Work

We have extended the accelerated proximal point method to the solution of singly linearly constrained quadratic programming with box constraints. The new algorithm is proved to be globally convergent. Numerical results also show that the new algorithm performs well on medium-scale quadratic programs. On the other hand, Solving the subproblem by searching a root of a piecewise linear continuous function is a very cheap. Considering the good performance of the new Algorithm 5, we can apply it to the solution of the subproblems in decomposition methods for large-scale SVM problems. Moreover, a parallel version of this algorithm combined with the theory in [9, 10] is also a direction of our future research.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This work is supported by the National Natural Science Foundation of China (11101420, 71271204) and the Excellent Young Scientist Foundation of Shandong Province (2010BSE06047).