Vector Machines
Hsuan-Tien Lin ([email protected]) Chih-Jen Lin ([email protected])
Department of Computer Science and Information Engineering, National Taiwan University, Taipei 106, Taiwan
Ruby C. Weng ([email protected]) Department of Statistics,
National Chengchi University, Taipei 116, Taiwan
Abstract. Platt’s probabilistic outputs for Support Vector Machines (Platt, 2000) has been popular for applications that require posterior class probabilities. In this note, we propose an improved algorithm that theoretically converges and avoids numerical difficulties. A simple and ready-to-use pseudo code is included.
Keywords: Support Vector Machine, Posterior Probability
1. Introduction
Given training examples xi ∈ Rn, i = 1, . . . , l, labeled by yi ∈ {+1, −1},
the binary Support Vector Machine (SVM) computes a decision func-tion f(x) such that sign(f(x)) can be used to predict the label of any test example x.
Instead of predicting the label, many applications require a posterior class probability Pr(y = 1|x). Platt (2000) proposes approximating the posterior by a sigmoid function
Pr(y = 1|x) ≈ PA,B(f ) ≡
1
1 + exp(Af + B), where f = f(x). (1) Let each fi be an estimate of f(xi). The best parameter setting z∗ =
(A∗, B∗) is determined by solving the following regularized maximum likelihood problem (with N+ of the yi’s positive, and N− negative):
min z=(A,B) F (z) = − l X i=1 tilog(pi) + (1 − ti) log(1 − pi) , (2)
for pi= PA,B(fi), and ti =
( N ++1 N++2 if yi= +1 1 N−+2 if yi= −1 , i = 1, . . . , l.
Platt (2000) gives a pseudo code for solving (2). In this note, we show how the pseudo code could be improved. We analyze (2) in Section 2,
and propose a more robust algorithm to solve it. Better implementation that avoids numerical difficulties is then discussed in Section 3. We compare our algorithm with Platt’s in Section 4. Finally, a ready-to-use pseudo code is in Appendix C.
2. Choice of Optimization Algorithm
We first introduce the simple optimization algorithm used in Platt’s pseudo code (Platt, 2000). Then, after proving that (2) is a convex optimization problem, we propose a more robust algorithm that enjoys similar simplicity, and theoretically converges.
2.1. Platt’s Approach: Levenberg-Marquardt Method
Platt (2000) uses a Levenberg-Marquardt (LM) algorithm from Press et al. (1992) to solve (2). The LM method was originally designed for solving nonlinear least-square problems. As an iterative procedure, at the k-th step, this method solves ( ˜Hk+ λkI)δk = −∇F (zk) to obtain
a direction δk, and moves the solution from zk to zk+1= zk+ δk if the
function value is sufficiently decreased. Here, ˜Hk is a special
approx-imation of the Hessian of the least-square problem, I is the identity matrix, and {zk}∞k=0 is the sequence of iteration vectors. When λk is
large, δk is close to the negative gradient direction. On the contrary, a
small λk leads δk to be more like a Newton’s direction.
In the pseudo code, Platt (2000) adapts the following rule for up-dating λk (Press et al., 1992):
If F (zk+ δk) < F (zk) then λk+1 ← 0.1 · λk ; Else λk+1 ← 10 · λk.
That is, if the new solution decreases the function value, λk is
reduced, and in the next iteration a more aggressive direction like Newton’s is tried. Otherwise, δk is unacceptable so we increase λk to
obtain a shorter vector which, more likely being a descent direction, may lower the function value.
Unfortunately, such an implementation may not converge to the minimum solution of (2). To the best of our knowledge, existing con-vergence proofs (e.g., Mor´e, 1978) all require more complicated or more robust rules for updating λk.
In fact, since (2) is not exactly a least-squares problem, the imple-mentation of Platt (2000) aims for general unconstrained optimization. It is known (e.g., Fletcher, 1987) that for unconstrained optimization we should avoid directly dealing with λk. Instead, the update of λk can
Thus, currently the optimization community uses trust-region methods for unconstrained optimization and the LM method is considered as its “progenitor” (Nocedal and Wright, 1999).
The LM-type implementation of Platt (2000) has one advantage: simplicity. However, the above discussion shows that it may not be the best choice for solving (2). Next, we propose an algorithm that is also simple, but enjoys better convergence properties.
2.2. Our Approach: Newton’s Method with Backtracking
As indicated by Platt (2000), any method for unconstrained optimiza-tion can be used for solving (2). Before we choose a suitable method, we analyze the optimization problem in more detail. First, the gradi-ent ∇F (z) and the Hessian matrix H(z) = ∇2F (z) are computed:
∇F (z) = Pl i=1fi(ti− pi) Pl i=1(ti− pi) , H(z) = Pl i=1fi2pi(1 − pi) Pli=1fipi(1 − pi) Pl i=1fipi(1 − pi) Pl i=1pi(1 − pi) .
Some analysis of this Hessian matrix is in the following theorem:
Theorem 1 The Hessian matrix H(z) is positive semi-definite. In ad-dition, H(z) is positive definite if and only if min
1≤i≤lfi 6= max1≤i≤lfi.
The proof is in Appendix A. Therefore, problem (2) is convex (and in general strictly convex). With such a nice property, we decide to use a simple Newton’s method with backtracking line search (Algorithm 6.2, Nocedal and Wright, 1999, and Section 10.5, Nash and Sofer, 1996). Though the trust-region type method mentioned in the end of Sec-tion 2.1 may be more robust, the implementaSec-tion is more complicated. For this two-variable optimization problem, simplicity is important, and hence trust-region methods are less favorable.
Our proposed algorithm is in Algorithm 1. As Hk = H(zk) may be
singular, a small positive diagonal matrix is added to the Hessian. With ∇F (zk)Tδk= −∇F (zk)T(Hk+ σI)−1∇F (zk) < 0,
the step size αk can be backtracked until the sufficient decrease
condi-tion (3) is satisfied.
If H(z) is positive definite for all z, the convergence of Algorithm 1 can be established from, for example, Theorem 10.2 by Nash and Sofer (1996). A simplified statement is shown in Theorem 2.
Input: Initial point z0, and parameter σ ≥ 0 such that H(z) + σI is
positive definite for all z 1: for k = 0, 1, 2, · · · do
2: Solve (Hk+ σI)δk= −∇F (zk)
3: Find αk as the first element of the sequence 1,12,14, · · · to satisfy
F (zk+ αkδk) ≤ F (zk) + 0.0001 · αk ∇F (zk)Tδk
(3)
4: Set zk+1 = zk+ αkδk
5: end for
Algorithm 1: Newton’s method with backtracking line search
Theorem 2 (Convergence of Algorithm 1 for general F (z))
If F (z) is twice continuously differentiable, H(z) is positive definite for all z, and F (z) attains an optimal solution at z∗, then limk→∞zk = z∗.
From Theorem 1, in some rare situations, H(z) is positive semi-definite but not positive semi-definite. Then, Theorem 2 cannot be directly applied. In Appendix B, we show that if σ > 0, Algorithm 1 still con-verges to an optimal solution. Therefore, we get the following theorem:
Theorem 3 (Convergence of Algorithm 1 for (2))
If Algorithm 1 is applied to (2) such that H(z) + σI is always positive definite, then limk→∞zk exists and is a global optimal solution.
3. Numerical Implementation
Next, we study the numerical difficulties that arise when solving (2) using Platt’s pseudo code. Then, we show our implementation that avoids the difficulties.
3.1. Platt’s Implementation
Platt (2000) uses the following pseudo code to calculate the objective value of (2) for a new pair of (A, B):
for i = 1 to len {
p = 1/(1+exp(deci[i]*A+B))
// At this step, make sure log(0) returns -200 err -= t*log(p)+(1-t)*log(1-p)
Here, len is l, the number of examples used, and err is the objective value. In addition, deci[i] is fi, and hence p stores the calculated pi.
However, t was lastly assigned to tlbefore this loop, and the calculation
does not use all ti, i = 1, . . . , l. Therefore, this pseudo code does not
correctly calculate the objective value of (2).
Furthermore, the above code assumes that log(0) returns −200, which reveals possible numerical difficulties:
1. log and exp could easily cause an overflow. If Afi + B is large,
exp(Afi + B) → ∞. In addition, when pi is near zero, log(pi) →
−∞. Although these problems do not always happen, consider-ing log(0) to be −200 is not a good solution.
2. 1−pi= 1−1+exp(Af1 i+B) is a “catastrophic cancellation” (Goldberg,
1991) when pi is close to one. That is, when subtracting two nearby
numbers that are already results of floating-point operations, the relative error can be so large that most digits are meaningless. For example, if fi = 1, and (A, B) = (−64, 0), in a simple C++
program with double precision, 1−pireturns zero but its equivalent
form exp(Afi+B)
1+exp(Afi+B) gives a more accurate result. This catastrophic
cancellation actually introduces most of the log(0) occurrences. Almost all algorithms that solve (2) need to face these issues. Next, we will discuss some techniques to resolve them.
3.2. Our Implementation
A problem of catastrophic cancellation can usually be resolved by reformulation: −tilog pi+ (1 − ti) log(1 − pi) (4) = (ti− 1)(Afi+ B) + log 1 + exp(Afi+ B) (5) = ti(Afi+ B) + log 1 + exp(−Afi− B) . (6)
With (5) or (6), 1−pidoes not appear. Moreover, log(0) never happens.1
Note, however, that even if (5) or (6) is used, the overflow problem may still occur. The problem is not serious if the IEEE floating-point standard is supported (Goldberg, 1991): an overflow leads to a special number INF, which can still be used in further operations. For example,
1
As pointed out by a reviewer, in many popular languages, log(1+...) can be replaced by log1p(...) to compute the result more accurately when the operand exp(Afi+ B) or exp(−Afi− B) is close to zero.
if a large αk in Line 3 of Algorithm 1 makes the exp operation of (5)
to overflow for some Afi + B, the new objective value would also be
evaluated as INF. Then, under the IEEE standard, INF is bigger than the current F (zk), and hence αk is reduced to a smaller value, with
which Afi+ B may not cause an overflow again.
Furthermore, regardless of whether the IEEE standard is supported, we can replace an overflow operation with an underflow one, a rule-of-thumb which has been frequently used in numerical computation. In general, an underflow is much less disastrous than an overflow. Therefore, we propose implementing (4) with the rule:
If Afi+ B ≥ 0 then use (6); Else use (5).
In addition, we can evaluate (1) by a similar trick:
If Af + B ≥ 0 then use 1+exp(−Af −B)exp(−Af −B) ; Else use (1).
The trick can be used in calculating ∇F (z) and H(z) as well: The term 1 − pi in H(z) can also cause a catastrophic cancellation. An easy
solution is to replace 1 − pi with the rule:
If Afi+ B ≥ 0 then use 1+exp(−Af1 i−B); Else use 1+exp(Afexp(Afi+B)i+B).
4. Experiment
We implemented Platt’s pseudo code (Platt, 2000), fixed the bug that was discussed in the beginning of Subsection 3.1, and compared it to our proposed algorithm. For fairness, both algorithms were realized in python, and were set with a stopping condition k∇F (zk)k∞< 10−5.
For the value of σ in Algorithm 1, we considered two approaches: − fixed: use a small fixed σ = 10−12.
− dynamic: apply Theorem 1 to check whether H(z) is positive def-inite, and set σ = 0 instead if the condition is true.
We compared the algorithms on two UCI data sets, sonar and shut-tle (Newman et al., 1998). Only classes 2 and 4 were taken from shutshut-tle to form a binary problem. The values fiwere generated with the scaled
data sets by LIBSVM using the RBF kernel (Chang and Lin, 2001). The soft-margin parameter log2C was varied in −5, −3, · · ·, 15, and the kernel parameter log2γ was varied in −15, −13, · · ·, 3. That is, 110 different problems (2) were tested for each data set.
Tables I and II list the average results for each data set. We first compared each algorithm based on the number of overflow errors en-countered, the number of iterations, and the final objective value F (z).
Table I. Average results of different algorithms for solving (2) on sonar algorithm # overflow errors # iterations final F (z) # backtracking steps per iteration
Platt’s 0 5.77 107.78 —
ours, fixed 0 5.56 107.78 0
ours, dynamic 0 5.56 107.78 0
Table II. Average results of different algorithms for solving (2) on shuttle algorithm # overflow
errors # iterations
final F (z)
# backtracking steps per iteration
Platt’s 589.30 8.00 158.62 —
ours, fixed 0 6.66 157.83 0.17
ours, dynamic 0 6.68 157.83 0.24
While Platt’s algorithm did reasonably well on sonar, it encountered numerous overflow errors on shuttle, needed more iterations, and some-times could not return a solution with decent F (z). On the other hand, our proposed algorithm worked well on both data sets.
The number of backtracking steps per iteration was also listed for the two approaches of setting σ. We can see that the fixed approach needed less backtracking steps per iteration on shuttle. The benefit came from the regularization on some nearly singular H(z). In addition, the fixed approach is simpler to implement in practice, and hence shall be preferred.
Finally, a simple and robust code is in Appendix C. It has been integrated into LIBSVM since version 2.6 (Chang and Lin, 2001). Source code in several popular languages can be downloaded at http://www. csie.ntu.edu.tw/~cjlin/libsvmtools.
Acknowledgment
We thank John Platt, S. Sathiya Keerthi, and the anonymous reviewers for helpful comments.
References
Chang, C.-C. and C.-J. Lin: 2001, ‘LIBSVM: a library for support vector machines’. Software available at http://www.csie.ntu.edu.tw/~cjlin/libsvm.
Fletcher, R.: 1987, Practical Methods of Optimization. John Wiley and Sons. Goldberg, D.: 1991, ‘What every computer scientist should know about
Mor´e, J. J.: 1978, ‘The Levenberg-Marquardt algorithm: Implementation and theory’. In: G. Watson (ed.): Numerical Analysis. New York, pp. 105–116. Nash, S. G. and A. Sofer: 1996, Linear and Nonlinear Programming. McGraw-Hill. Newman, D. J., S. Hettich, C. L. Blake, and C. J. Merz: 1998, ‘UCI Repository of machine learning databases’. Technical report, University of California, Irvine, Dept. of Information and Computer Sciences.
Nocedal, J. and S. J. Wright: 1999, Numerical Optimization. New York, NY: Springer-Verlag.
Platt, J.: 2000, ‘Probabilistic outputs for support vector machines and comparison to regularized likelihood methods’. In: A. Smola, P. Bartlett, B. Sch¨olkopf, and D. Schuurmans (eds.): Advances in Large Margin Classifiers. Cambridge, MA. Press, W. H., B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling: 1992, Numerical
Recipes: The Art of Scientific Computing. Cambridge (UK) and New York: Cambridge University Press, 2nd edition.
Appendix
A. Proof of Theorem 1
Since the definition of pi in (1) implies that 0 < pi < 1, we can
define vectors u and v with ui = fippi(1 − pi), and vi =ppi(1 − pi),
respectively. Then H(z) = u Tu uTv vTu vTv . By Cauchy inequality, detH(z)= l X i=1 u2i ! l X i=1 v2i ! − l X i=1 uivi !2 ≥ 0. (7)
Since the two diagonal terms and the determinant are all nonnegative, the matrix H(z) is positive semi-definite.
From (7), det
H(z)
= 0 if and only if u and v are parallel vectors. Since ui = fivi and vi > 0, this situation happens if and only if
all fi’s are equal. That is, the matrix H(z) is positive definite if and
only if min
1≤i≤lfi 6= max1≤i≤lfi.
B. Proof of Theorem 3
Case 1: H(z) is always positive definite. If one can prove that
S =n(A, B): F (A, B) ≤ F (A0, B0)
o
(8)
is bounded, then F (A, B) attains an optimal solution within S and The-orem 2 can be applied to show the convergence.
From Theorem 1, assume without loss of generality that f1 6= f2. Let ˆa = f1 1 f2 1 A B . Since f1 1 f2 1
is invertible, it suffices to show that ˆS = {ˆa: (A, B) ∈ S} is bounded. If not, there exists an infinite sequence {ˆak} ∞ k=1 in ˆS such that lim k→∞max |(ˆak)1|, |(ˆak)2| = ∞.
Then, without loss of generality, there exists an infinite subsequence K such that limk→∞,k∈K|(ak)1| = ∞. However, since F (Ak, Bk) is the
summation of positive terms,
F (Ak, Bk) ≥ −t1log
1
1 + e(ˆak)1 − (1 − t1) log
e(ˆak)1
1 + e(ˆak)1.
The right-hand-side above goes to ∞ as |(ˆak)1| → ∞. Therefore, there
exists some k such that F (Ak, Bk) > F (A0, B0), which somehow
con-tradicts ˆak∈ ˆS. Thus, ˆS is bounded and the proof is complete.
Case 2: When H(z) is only positive semi-definite for some z, from Theorem 1, all fi’s are equal. By considering fi = f for all i, we can
define a = Af +B and a single-variable function ¯F (a) = F (A, B). Then
¯ F0(a) = l X i=1 ti− l 1 + ea, ¯F 00 (a) = le a (1 + ea)2.
By simplifying (3), in Algorithm 1, (H(z) + σI)δ = −∇F (z) is
lea (1+ea)2 f2 f f 1 + σI (δ)1 (δ)2 = − l X i=1 ti−1+el a ! f 1 . (9)
If σ > 0, the solution δ satisfies (δ)1 = f · (δ)2. Then, the first (and the
second) equation of the linear system (9) is the same as ¯ F00(a) + σ f2+ 1 f · (δ)1+ (δ)2 = − ¯F0(a). (10)
Interestingly, if we apply Algorithm 1 to minimize ¯F (a) with f2σ+1
added to its Hessian ¯F00(a), equation (10) is exactly the linear system to be solved. Therefore, if a0= A0f + B0, then for all k,
ak+1 = ak+ αk f · (δk)1+ (δk)2 = Ak+ αk(δk)1 f +Bk+ αk(δk)2 . (11)
Since ¯F (a) is strictly convex from ¯F00(a) > 0, similar techniques in Case 1 can be used to prove that ¯F (a) attains an optimal solution. Therefore, from Theorem 2, the sequence {ak}∞k=0 globally converges.
Then, from (δk)1= f · (δk)2 and (11),
lim k→∞ak= (A0f + B0) + (f 2+ 1) ∞ X k=0 αk(δk)2
exists. Therefore, limk→∞Bk = B0+P∞k=0αk(δk)2 exists, and so does
limk→∞Ak. In addition, they form an optimal solution of
minimiz-ing F (A, B).
C. Pseudo Code of Algorithm 1
We recommend using double precision for the algorithm. Input parameters:
deci = array of SVM decision values
label = array of booleans: is the example labeled +1? prior1 = number of positive examples
prior0 = number of negative examples Outputs:
A, B = parameters of sigmoid
//Parameter setting
maxiter=100 //Maximum number of iterations minstep=1e-10 //Minimum step taken in line search sigma=1e-12 //Set to any value > 0
//Construct initial values: target support in array t,
// initial function value in fval
hiTarget=(prior1+1.0)/(prior1+2.0), loTarget=1/(prior0+2.0) len=prior1+prior0 // Total number of data
for i = 1 to len { if (label[i] > 0) t[i]=hiTarget else t[i]=loTarget }
A=0.0, B=log((prior0+1.0)/(prior1+1.0)), fval=0.0 for i = 1 to len {
fApB=deci[i]*A+B if (fApB >= 0)
fval += t[i]*fApB+log(1+exp(-fApB)) else
fval += (t[i]-1)*fApB+log(1+exp(fApB)) }
for it = 1 to maxiter {
//Update Gradient and Hessian (use H’ = H + sigma I) h11=h22=sigma, h21=g1=g2=0.0 for i = 1 to len { fApB=deci[i]*A+B if (fApB >= 0) p=exp(-fApB)/(1.0+exp(-fApB)), q=1.0/(1.0+exp(-fApB)) else p=1.0/(1.0+exp(fApB)), q=exp(fApB)/(1.0+exp(fApB)) d2=p*q h11 += deci[i]*deci[i]*d2, h22 += d2, h21 += deci[i]*d2 d1=t[i]-p g1 += deci[i]*d1, g2 += d1 }
if (abs(g1)<1e-5 && abs(g2)<1e-5) //Stopping criteria break
//Compute modified Newton directions det=h11*h22-h21*h21
dA=-(h22*g1-h21*g2)/det, dB=-(-h21*g1+h11*g2)/det gd=g1*dA+g2*dB
stepsize=1
while (stepsize >= minstep){ //Line search
newA=A+stepsize*dA, newB=B+stepsize*dB, newf=0.0 for i = 1 to len { fApB=deci[i]*newA+newB if (fApB >= 0) newf += t[i]*fApB+log(1+exp(-fApB)) else newf += (t[i]-1)*fApB+log(1+exp(fApB)) } if (newf<fval+0.0001*stepsize*gd){ A=newA, B=newB, fval=newf
break //Sufficient decrease satisfied }
else
stepsize /= 2.0 }
if (stepsize < minstep){ print ’Line search fails’
break }
}
if (it >= maxiter)
print ’Reaching maximum iterations’
