Shark源码分析（十二）：线性SVM

本文主要是介绍Shark源码分析（十二）：线性SVM，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

Shark源码分析（十二）：线性SVM

关于svm算法，这个在我关于机器学习的博客中已经描述的比较详实了，这里就不再赘述。svm主要有三种类型，这里我所介绍的是线性svm算法的代码。相较于使用核函数的svm算法，代码的整体框架应该是一样的，只是在对偶问题的求解上所使用的方法可能是不一样的。

LinearClassifier类

这个类所表示的是算法的决策平面，是一个多分类的线性分类模型。定义在<include/shark/Models/LinearClassifier.h>中。

template<class VectorType = RealVector>
class LinearClassifier : public ArgMaxConverter<LinearModel<VectorType> >
{
public:LinearClassifier(){}std::string name() const{ return "LinearClassifier"; }
};

相当简单的一个类，并没有什么好说明的地方。

ArgMaxConverter类

该类是LinearClassifier的基类，其作用是将一个输出的向量通过arg_max操作转变为一个类标记，就是输出分量最大的那一维。该类定义在<include/shark/Models/Converter.h>。

template<class Model>
class ArgMaxConverter : public AbstractModel<typename Model::InputType, unsigned int>
{
private:typedef typename Model::BatchOutputType ModelBatchOutputType;
public:typedef typename Model::InputType InputType;typedef unsigned int OutputType;typedef typename Batch<InputType>::type BatchInputType;typedef Batch<unsigned int>::type BatchOutputType;ArgMaxConverter(){ }ArgMaxConverter(Model const& decisionFunction): m_decisionFunction(decisionFunction){ }std::string name() const{ return "ArgMaxConverter<"+m_decisionFunction.name()+">"; }RealVector parameterVector() const{return m_decisionFunction.parameterVector();}void setParameterVector(RealVector const& newParameters){m_decisionFunction.setParameterVector(newParameters);}std::size_t numberOfParameters() const{return m_decisionFunction.numberOfParameters();}Model const& decisionFunction()const{return m_decisionFunction;}Model& decisionFunction(){return m_decisionFunction;}// 计算输入数据的类标签void eval(BatchInputType const& input, BatchOutputType& output)const{ModelBatchOutputType modelResult;m_decisionFunction.eval(input,modelResult);std::size_t batchSize = shark::size(modelResult);output.resize(batchSize);if(modelResult.size2()== 1) //对于二分类的情况{for(std::size_t i = 0; i != batchSize; ++i){// 如果输出大于0表示正类，否则为负类output(i) = modelResult(i,0) > 0.0;}}else{for(std::size_t i = 0; i != batchSize; ++i){output(i) = static_cast<unsigned int>(arg_max(row(modelResult,i)));}}}void eval(BatchInputType const& input, BatchOutputType& output, State& state)const{eval(input,output);}void eval(InputType const & pattern, OutputType& output)const{typename Model::OutputType modelResult;m_decisionFunction.eval(pattern,modelResult);if(modelResult.size()== 1){output = modelResult(0) > 0.0;}else{output = static_cast<unsigned int>(arg_max(modelResult));}}void read(InArchive& archive){archive >> m_decisionFunction;}void write(OutArchive& archive) const{archive << m_decisionFunction;}private:Model m_decisionFunction;
};

在LinearClassifier类的代码中，该模板类的模板参数是LinearModel，这个模板类我们之前已经介绍过了。

AbstractLinearSvmTrainer类

这个类是所有线性svm训练方法的基类。该类定义在<include/shark/Algorithms/Trainers/AbstractSvmTrainer.h>中。

template <class InputType>
class AbstractLinearSvmTrainer
: public AbstractTrainer<LinearClassifier<InputType>, unsigned int>
, public QpConfig
, public IParameterizable
{
public:typedef AbstractTrainer<LinearClassifier<InputType>, unsigned int> base_type;typedef LinearClassifier<InputType> ModelType;AbstractLinearSvmTrainer(double C, bool unconstrained = false): m_C(C), m_unconstrained(unconstrained){ RANGE_CHECK( C > 0 ); }double C() const{ return m_C; }void setC(double C) {RANGE_CHECK( C > 0 );m_C = C;}bool isUnconstrained() const{ return m_unconstrained; }RealVector parameterVector() const{RealVector ret(1);ret(0) = (m_unconstrained ? std::log(m_C) : m_C);return ret;}void setParameterVector(RealVector const& newParameters){SHARK_ASSERT(newParameters.size() == 1);setC(m_unconstrained ? std::exp(newParameters(0)) : newParameters(0));}size_t numberOfParameters() const{ return 1; }// 对于以下的两个类成员，在QpConfig的构造函数中没有为它们赋值。稍后可以看到，它们自己的构造函数中是有默认值的using QpConfig::m_stoppingcondition; // 算法训练的停止条件using QpConfig::m_solutionproperties; // 当前解的一些性质using QpConfig::m_verbosity; // 冗长程度（字面翻译，在后面的代码中体现的不是这个意思），默认值是0，protected:double m_C; // 目标函数中正则化项的系数bool m_unconstrained; // 是否使用log C 代替了C，如果是的话则摆脱了C > 0的限制，并不知道这个有什么用
};

QpStoppingCondition类、QpStopType类和QpSolutionProperties类

这三个类都定义在<include/shark/Algorithms/QP/QuadraticProgram.h>中。

struct QpStoppingCondition
{QpStoppingCondition(double accuracy = 0.001, unsigned long long iterations = 0xffffffff, double value = 1e100, double seconds = 1e100){minAccuracy = accuracy;maxIterations = iterations;targetValue = value;maxSeconds = seconds;}//违反KKT条件的阈值下限double minAccuracy;//最大迭代次数unsigned long long maxIterations;//目标函数值的阈值double targetValue;//算法运行的最长时间double maxSeconds;
};

enum QpStopType
{QpNone = 0,QpAccuracyReached = 1,QpMaxIterationsReached = 4,QpTimeout = 8,
};

struct QpSolutionProperties
{QpSolutionProperties(){type = QpNone;accuracy = 1e100;iterations = 0;value = 1e100;seconds = 0.0;}QpStopType type;double accuracy;     //当前解违反KKT条件的程度unsigned long long iterations; //当前循环的次数double value;    // 当前目标函数的值double seconds; // 当前程序的运行时间
};

LinearCSvmTrainer类

该类就是用于训练线性SVM的，定义在<include/shark/Algorithms/Trainers/CSvmTrainer.h>。

template <class InputType>
class LinearCSvmTrainer : public AbstractLinearSvmTrainer<InputType>
{
public:typedef AbstractLinearSvmTrainer<InputType> base_type;LinearCSvmTrainer(double C, bool unconstrained = false) : AbstractLinearSvmTrainer<InputType>(C, unconstrained){}std::string name() const{ return "LinearCSvmTrainer"; }void train(LinearClassifier<InputType>& model, LabeledData<InputType, unsigned int> const& dataset){std::size_t dim = inputDimension(dataset);QpBoxLinear<InputType> solver(dataset, dim);RealMatrix w(1, dim, 0.0);row(w, 0) = solver.solve(base_type::C(),0.0,QpConfig::stoppingCondition(),&QpConfig::solutionProperties(),QpConfig::verbosity() > 0);model.decisionFunction().setStructure(w);}
};

从代码中可以看出，主要还是调用QpBoxLinear类的solve方法来求解。

QpBoxLinear类

该类是利用矩阵分解的方法来求解目标函数是hinge损失函数的线性svm，定义在<include/shark/Algorithms/QP/QpBoxLinear.h>。光看代码你可能会不清楚一些操作的具体含义，需要看一下”A Dual Coordinate Descent Method for Large-scale Linear SVM”这篇论文。

template <class InputT>
class QpBoxLinear
{
public:typedef LabeledData<InputT, unsigned int> DatasetType;typedef typename LabeledData<InputT, unsigned int>::const_element_reference ElementType;QpBoxLinear(const DatasetType& dataset, std::size_t dim): m_data(dataset), m_xSquared(m_data.size()), m_dim(dim){SHARK_ASSERT(dim > 0);for (std::size_t i=0; i<m_data.size(); i++){ElementType x_i = m_data[i];m_xSquared(i) = inner_prod(x_i.input, x_i.input);}}// 参数reg相当于论文中的D_{ii}RealVector solve(double bound,double reg,QpStoppingCondition& stop,QpSolutionProperties* prop = NULL,bool verbose = false){SHARK_ASSERT(bound > 0.0);SHARK_ASSERT(reg >= 0.0);Timer timer;timer.start();std::size_t ell = m_data.size();RealVector alpha(ell, 0.0); // 表示拉格朗日乘子RealVector w(m_dim, 0.0); // 权值向量RealVector pref(ell, 1.0);          // measure of success of individual stepsdouble prefsum = ell;               // normalization constantstd::vector<std::size_t> schedule(ell); // 更新每一个拉格朗日乘子的顺序// prepare countersstd::size_t epoch = 0;std::size_t steps = 0;// prepare performance monitoring for self-adaptationdouble max_violation = 0.0;const double gain_learning_rate = 1.0 / ell;double average_gain = 0.0;bool canstop = true;// outer optimization loopwhile (true){// 计算更新的下标顺序，至于这种算法的原理我就不是很清楚了，论文里也没有说明double psum = prefsum;prefsum = 0.0;std::size_t pos = 0;for (std::size_t i=0; i<ell; i++){double p = pref[i];double num = (psum < 1e-6) ? ell - pos : std::min((double)(ell - pos), (ell - pos) * p / psum);std::size_t n = (std::size_t)std::floor(num);double prob = num - n;if (Rng::uni() < prob) n++;for (std::size_t j=0; j<n; j++){schedule[pos] = i;pos++;}psum -= p;prefsum += p;}SHARK_ASSERT(pos == ell);for (std::size_t i=0; i<ell; i++) std::swap(schedule[i], schedule[Rng::discrete(0, ell - 1)]);// inner loop// 算法的符号与论文中是相反的，包括g，pg和new_a的计算max_violation = 0.0;for (std::size_t j=0; j<ell; j++){// active variablestd::size_t i = schedule[j];ElementType e_i = m_data[i];double y_i = (e_i.label > 0) ? +1.0 : -1.0;// compute gradient and projected gradientdouble a = alpha(i);double wyx = y_i * inner_prod(w, e_i.input);double g = 1.0 - wyx - reg * a;double pg = (a == 0.0 && g < 0.0) ? 0.0 : (a == bound && g > 0.0 ? 0.0 : g);// update maximal KKT violation over the epochmax_violation = std::max(max_violation, std::abs(pg));double gain = 0.0;// 更新参数的过程if (pg != 0.0){// SMO-style coordinate descent stepdouble q = m_xSquared(i) + reg;double mu = g / q; // 该参数同时也表示了a的两次更新之间的差值double new_a = a + mu;// numerically stable updateif (new_a <= 0.0){mu = -a;new_a = 0.0;}else if (new_a >= bound){mu = bound - a;new_a = bound;}// 更新参数alpha(i) = new_a;w += (mu * y_i) * e_i.input;gain = mu * (g - 0.5 * q * mu);steps++;}// update gain-based preferences{if (epoch == 0) average_gain += gain / (double)ell;else{double change = CHANGE_RATE * (gain / average_gain - 1.0);double newpref = std::min(PREF_MAX, std::max(PREF_MIN, pref(i) * std::exp(change)));prefsum += newpref - pref(i);pref[i] = newpref;average_gain = (1.0 - gain_learning_rate) * average_gain + gain_learning_rate * gain;}}}epoch++;if (stop.maxIterations > 0 && ell * epoch >= stop.maxIterations) //这里的最大循环次数指的是内部循环的次数{if (prop != NULL) prop->type = QpMaxIterationsReached;break;}if (timer.stop() >= stop.maxSeconds){if (prop != NULL) prop->type = QpTimeout;break;}if (max_violation < stop.minAccuracy){if (verbose) std::cout << "#" << std::flush;if (canstop){if (prop != NULL) prop->type = QpAccuracyReached;break;}else{// prepare full sweep for a reliable checking of the stopping criterioncanstop = true;for (std::size_t i=0; i<ell; i++) pref[i] = 1.0;prefsum = ell;}}else{if (verbose) std::cout << "." << std::flush;canstop = false;}}timer.stop();// compute solution statisticsstd::size_t free_SV = 0; // 不在决策边界上的支持向量个数std::size_t bounded_SV = 0; // 在决策边界上的支持向量的个数double objective = -0.5 * shark::blas::inner_prod(w, w); //计算最终的目标函数值，但计算的值有些诡异，既不是原问题的目标值也不是对偶问题的目标值for (std::size_t i=0; i<ell; i++){double a = alpha(i);if (a > 0.0){objective += a;objective -= reg/2.0 * a * a;if (a == bound) bounded_SV++;else free_SV++;}}// return solution statisticsif (prop != NULL){prop->accuracy = max_violation;       // this is approximate, but a good guessprop->iterations = ell * epoch;prop->value = objective;prop->seconds = timer.lastLap();}// output solution statistics// 这里verbose只是一个是否需要输出信息的标志if (verbose){std::cout << std::endl;std::cout << "training time (seconds): " << timer.lastLap() << std::endl;std::cout << "number of epochs: " << epoch << std::endl;std::cout << "number of iterations: " << (ell * epoch) << std::endl;std::cout << "number of non-zero steps: " << steps << std::endl;std::cout << "dual accuracy: " << max_violation << std::endl;std::cout << "dual objective value: " << objective << std::endl;std::cout << "number of free support vectors: " << free_SV << std::endl;std::cout << "number of bounded support vectors: " << bounded_SV << std::endl;}// return the solutionreturn w;}protected:DataView<const DatasetType> m_data; // 训练数据            RealVector m_xSquared; //m_data^T m_data                      std::size_t m_dim; // 输入数据的维度              
};