libsvm代码阅读：关于Solver类分析（一）

update：2014-2-27 LinJM @HQU  『 libsvm专栏地址：http://blog.csdn.net/column/details/libsvm.html 』

现在我们涉及到的Solver类是一个SVM优化求解的实现技术：SMO，即序贯最小优化算法。libsvm中最原始的Solver的代码有六百多行，再加上各种变形就上千行了，为了好理解，我们先来看看理论问题。

代码的开头如下：

// An SMO algorithm in Fan et al., JMLR 6(2005), p. 1889--1918
// Solves:
//
//	min 0.5(\alpha^T Q \alpha) + p^T \alpha
//
//		y^T \alpha = \delta
//		y_i = +1 or -1
//		0 <= alpha_i <= Cp for y_i = 1
//		0 <= alpha_i <= Cn for y_i = -1

即代码的实现主要参考文献：

Fan R E, Chen P H, Lin C J. Working set selection using second order information for training support vector machines[J]. The Journal of Machine Learning Research, 2005, 6: 1889-1918.
那么，我们就来读读这篇文章：

SVM的优化问题如下：



核心目标是：求出最优解alpha*。



分解算法只更新拉格朗日乘子alpha_i的一个固定大小的子集，其他保持不变。因此，每当更新一个新点加入到工作集，另一个点要被移除。在这个算法中，目标是每次在数据的一个小的子集上优化全局问题。

而SMO算法是将分解算法思想推向极致得出的，即每次迭代仅优化两个点的最小子集。



算法的流程如上所示。但是，如何find a two-element working set B，上面的算并没有讲。下面我们直接贴出本篇论文的Select_B 方法：







上面的算法的流程还是相对清晰的。下面，我们来看看SMO的伪代码（Algorithm 2）

Inputs:
y: array of {+1, -1}: class of the i-th instance
Q: Q[i][j] = y[i]*y[j]*K[i][j]; K: kernel matrix
len: number of instances

//parameters
eps = 1e-3 // stopping tolerance
tau = 1e-12

//main routine
initialize alpha array A to all zero
initialize gradient array G to all -1

while(1)
{
(i,j) = selectB()
if (j == -1)
break

//working set is (i,j)
a = Q[i][i]+Q[j][j]-2*y[i]y[j]*Q[i][j]
if (a <= 0)
a = tau
b = -y[i]*G[i]+y[j]*G[j]

//update alpha
oldAi = A[i], oldAj = A[j]
A[i] += y[i]*b/a
A[j] -= y[j]*b/a

//project alpha back to the feasible region
sum = y[i]*oldAi + y[j]*oldAj
if A[i] > C
A[i] = C
if A[i] < 0
A[i] = 0
A[j] = y[j]*(sum - y[i]*A[i])

if A[j] > C
A[j] = C
if A[j] < 0
A[j] = 0
A[i] = y[i]*(sum - y[j]*A[j])

//update gradient
deltaAi = A[i] - oldAi,   deltaAj = A[j] - oldAj
for t = 1 to len
G[t] += Q[t][i]*deltaAi + Q[t][j]*deltaAj
}

procedure selectB
//select i
i = -1
G_max = -inf
G_min = inf
for t = 1 to len
if(y[t]==+1 and A[t] < C) or (y[t]==-1 and A[t] >0)
{
if(-y[t]*G[t] >= G_max)
{
i = t
G_max = -y[t]*G[t]
}
}

//select j
j = -1
obj_min = inf
for t = 1 to len
{
if(y[t]==+1 and A[t] >0)or(y[t]==-1 and A[t] < C)
{
b = G_max + y[t]*G[t]
if (-y[t]*G[t] <= G_min)
G_min = -y[t]*G[t]
if (b > 0)
{
a = Q[i][i]+Q[t][t]-2*y[i]*y[t]*Q[i][t]
if (a <= 0)
a = tau
if (-(b*b)/a <= obj_min)
{
j = t
obj_min = -(b*b)/a
}
}
}
}

if (G_max-G_min < eps)
return (-1,-1)

return (i,j)
end procedure

本文地址：http://blog.csdn.net/linj_m/article/details/19698463

微博：林建民-机器视觉