diff --git a/fastlib/trunk/contrib/houyang/regmin/build.py b/fastlib/trunk/contrib/houyang/regmin/build.py new file mode 100644 index 0000000000..e9073e7f62 --- /dev/null +++ b/fastlib/trunk/contrib/houyang/regmin/build.py @@ -0,0 +1,23 @@ +librule( + name = "svm", + + sources = ["regmin.cc"], + + headers = ["opt_smo.h", "opt_sgd.h", "opt_md.h", "opt_tgd.h", "regmin.h", "regmin_data.h"], + + deplibs = ["fastlib:fastlib"], + +) + +binrule( + # The name of the executable. + name = "regmin", + + sources = ["regmin.cc"], + + headers = ["opt_smo.h", "opt_sgd.h", "opt_md.h", "opt_tgd.h", "regmin.h", "regmin_data.h"], + + deplibs = [":svm"] + + ) + diff --git a/fastlib/trunk/contrib/houyang/regmin/opt_md.h b/fastlib/trunk/contrib/houyang/regmin/opt_md.h new file mode 100644 index 0000000000..4fc3c957de --- /dev/null +++ b/fastlib/trunk/contrib/houyang/regmin/opt_md.h @@ -0,0 +1,420 @@ +/** + * @author Hua Ouyang + * + * @file opt_md.h + * + * This head file contains functions for performing L1-regularized linear loss optimization, using Mirror Descent + * + * + * @see svm.h + */ + +#ifndef U_SVM_OPT_MD_H +#define U_SVM_OPT_MD_H + +#include "fastlib/fastlib.h" + +const double MD_ZERO = 1.0e-30; + +template +class MD { + FORBID_ACCIDENTAL_COPIES(MD); + + public: + typedef TKernel Kernel; + + private: + int learner_typeid_; + + Kernel kernel_; + const Dataset_sl *dataset_; + index_t n_data_; /* number of data samples */ + index_t n_features_; /* # of features == # of row - 1, exclude the last row (for labels) */ + // index_t n_features_bias_; /* # of features + 1 , [x, 1], for the bias term */ + + index_t n_sv_; /* number of support vectors */ + + index_t w_nnz_; + + double round_thd_; + + ArrayList y_; /* list that stores "labels" */ + + ArrayList w_; /* the slope of the decision hyperplane, including bias: [w, b] */ + ArrayList w_p_; /* coefficients for positive w_t+ */ + ArrayList w_n_; /* coefficients for negative w_t- */ + + double scale_w_; // the scale for w + + // parameters + double C_; // \|w\|_1^1 \leq C + double lambda_; // regularization parameter. lambda = 1/(C*n_data) + index_t n_iter_; // number of iterations + index_t n_epochs_; // number of epochs + double accuracy_; // accuracy for stopping creterion + double eta_; // step length. eta = 1/(lambda*t) + double t_; + + //bool is_constant_step_size_; // whether use constant step size (default) or not + + ArrayList old_from_new_; // for generating a random sequence of training data + + public: + MD() {} + ~MD() {} + + /** + * Initialization for parameters + */ + void InitPara(int learner_typeid, ArrayList ¶m_) { + // init parameters + if (learner_typeid == 0) { // SVM_C + C_ = param_[0]; + n_epochs_ = (index_t)param_[2]; + n_iter_ = (index_t)param_[3]; + accuracy_ = param_[4]; + } + else if (learner_typeid == 1) { // SVM_R + } + } + + void Train(int learner_typeid, Dataset_sl &dataset_in); + + Kernel& kernel() { + return kernel_; + } + + void GetW(ArrayList &w_out) { + index_t wp_size = w_p_.size(); + index_t ct_nz = 0; + w_out.Init(w_nnz_); + for (index_t i=0; i= round_thd_ ) { + w_out[ct_nz].index = w_p_[i].index; + w_out[ct_nz].value = w_p_[i].value; + ct_nz ++; + } + } + } + + double ScaleW() const { + return scale_w_; + } + + //void GetSV(ArrayList &dataset_index, ArrayList &coef, ArrayList &sv_indicator); + + private: + /** + * Loss functions + */ + double LossFunction_(int learner_typeid, double yy_hat) { + if (learner_typeid_ == 0) { // SVM_C + return HingeLoss_(yy_hat); + } + else if (learner_typeid_ == 1) { // SVM_R + return 0.0; // TODO + } + else + return HingeLoss_(yy_hat); + } + + /** + * Gradient of loss functions + */ + double LossFunctionGradient_(int learner_typeid, double yy_hat) { + if (learner_typeid_ == 0) { // SVM_C + return HingeLossGradient_(yy_hat); + //return LogisticLossGradient_(yy_hat); + } + else if (learner_typeid_ == 1) { // SVM_R + return 0.0; // TODO + } + else { + if (yy_hat < 1.0) + return 1.0; + else + return 0.0; + } + } + + /** + * Hinge Loss function + */ + double HingeLoss_(double yy_hat) { + if (yy_hat < 1.0) + return 1.0 - yy_hat; + else + return 0.0; + } + + /** + * Gradient of the Hinge Loss function + */ + double HingeLossGradient_(double yy_hat) { + if (yy_hat < 1.0) + return -1.0; + else + return 0.0; + } + + /** + * Gradient of the Logistic Loss function + */ + double LogisticLossGradient_(double yy_hat) { + double tmp = exp(-yy_hat); + return -tmp/(1+tmp); + } + + void LearnersInit_(int learner_typeid); + + int TrainIteration_(); + + double GetC_(index_t i) { + return C_; + } + +}; + + +/** + * Initialization according to different SVM learner types + * + * @param: learner type id + */ +template +void MD::LearnersInit_(int learner_typeid) { + index_t i; + double C_iv_f = C_/(n_features_); + double C_iv_2f = C_iv_f/2.0; + + learner_typeid_ = learner_typeid; + + // init w, w+, w- + w_.Init(0); + w_p_.Init(n_features_); + w_n_.Init(n_features_); + for (i=0; iy)[i] > 0 ? 1 : -1; + } +} + + +/** +* L1-regularization training for 2-classes +* +* @param: input 2-classes data matrix with labels (1,-1) in the last row +*/ +template +void MD::Train(int learner_typeid, Dataset_sl &dataset_in) { + index_t i, j, epo, ct; + + index_t total_n_iter; + //double DX, M, x_sq_sup; + double cons_step, w_sum = 1.0; + + /* general learner-independent initializations */ + dataset_ = &dataset_in; + n_data_ = dataset_->n_points; + n_features_ = dataset_->n_features; + + printf("n_data=%d,n_feature=%d\n", n_data_, n_features_); + + if (n_epochs_ > 0) { // # of epochs provided, use it + n_iter_ = n_data_; + total_n_iter = n_iter_ * n_epochs_; + } + else { // # of epochs not provided, use n_iter_ to count iterations + n_epochs_ = 1; // not exactly one epoch, just use it for one loop + total_n_iter = n_iter_; + } + + DEBUG_ASSERT(C_ != 0); + lambda_ = 1.0/(C_*n_data_); + + /* learners initialization */ + LearnersInit_(learner_typeid); + old_from_new_.Init(n_data_); + + // determine step sizes + double x_sq_sp = 0.0; + NZ_entry *xs; + xs = (dataset_->x)[0]; + while (xs->index != -1) { + x_sq_sp += math::Sqr(xs->value); + ++xs; + } + //eta_ = 1/(4.0*x_sq_sp); + t_ = 4.0 * x_sq_sp; + + cons_step = 1; + + index_t work_idx_old = 0; + + /* Begin training iterations */ + double yt, yt_hat, yy_hat; + scale_w_ = 1.0; + + printf("MD training begins...\n"); + for (epo = 0; epox)[work_idx_old]; + yt = y_[work_idx_old]; + yt_hat = SparseDot(w_p_, xt) - SparseDot(w_n_, xt); + yy_hat = yt * yt_hat; + + double eta_grad_tmp = eta_* yt * LossFunctionGradient_(learner_typeid, yy_hat); + //printf("epo:%d, ct:%d, eta_grad:%f\n", epo, ct, eta_grad_tmp); + if (eta_grad_tmp != 0) { + xt = (dataset_->x)[work_idx_old]; + index_t ct_w = 0; + + while (ct_w < w_p_.size() && xt->index != -1) { + if (w_p_[ct_w].index < xt->index) { + ++ct_w; + } + else if (w_p_[ct_w].index == xt->index) { + w_p_[ct_w].value = w_p_[ct_w].value * exp(-eta_grad_tmp * xt->value); + ++ct_w; + ++xt; + } + else { // w_p_[ct_w].index > xt->index + ++xt; + } + } + ct_w = 0; + xt = (dataset_->x)[work_idx_old]; + while (ct_w < w_n_.size() && xt->index != -1) { + if (w_n_[ct_w].index < xt->index) { + ++ct_w; + } + else if (w_n_[ct_w].index == xt->index) { + w_n_[ct_w].value = w_n_[ct_w].value * exp(eta_grad_tmp * xt->value); + ++ct_w; + ++xt; + } + else { // w_n_[ct_w].index > xt->index + ++xt; + } + } + + // calc sum_i(w_p_i+w_n_i) + w_sum = 0.0; + for (i=0; i C_) { + // printf("epo:%d,iter:%d, w_sum=%f\n", epo, ct, w_sum); + SparseScale(C_/w_sum, w_p_); + SparseScale(C_/w_sum, w_n_); + } + + } + + //printf("epo:%d, ct:%d-- sacle_w:%f\n", epo, ct, scale_w_); + + /* + for (i=0; i= round_thd_) { // TODO: thresholding + w_nnz_++; + } + } + + printf("%d out of %d features are non zero. NZ rate:%f\n", w_nnz_, n_features_, (double)(w_nnz_)/(double)n_features_); + + /* + // find max(abs(w_i)) + double wi_abs_max = -INFINITY; + double wi_abs; + for (i=0; i wi_abs_max) { + wi_abs_max = wi_abs; + } + } + // round small w_i to 0 + index_t w_ct = 0; + double round_factor = fx_param_double(NULL, "round_factor", 1.0e32); + double round_thd = wi_abs_max / round_factor; + for (i=0; i round_thd ) { + w_ct ++; + //printf("w_dim:%d, w_value:%lf\n", i, w_[i]); + } + else { + w_.Remove(i); + } + } + printf("%d out of %d features are non zero\n", w_ct, n_features_); + */ + + // Calculate objective value; default: no calculation to save time + int objvalue = fx_param_int(NULL, "objvalue", 0); + if (objvalue > 0) { + double hinge_loss = 0.0, loss_sum= 0.0; + + // primal objective value + for (i=0; i< n_data_; i++) { + hinge_loss = 1- y_[i] * SparseDot(w_p_, (dataset_->x)[i]); + if (hinge_loss > 0) { + loss_sum += hinge_loss; + } + } + + printf("Primal objective value: %lf\n", loss_sum); + } + +} + +#endif + + diff --git a/fastlib/trunk/contrib/houyang/regmin/opt_sgd.h b/fastlib/trunk/contrib/houyang/regmin/opt_sgd.h new file mode 100644 index 0000000000..4477cbac3d --- /dev/null +++ b/fastlib/trunk/contrib/houyang/regmin/opt_sgd.h @@ -0,0 +1,564 @@ +/** + * @author Hua Ouyang + * + * @file opt_sgd.h + * + * This head file contains functions for performing Stochastic Gradient Descent (SGD) based optimization for SVM + * + * The algorithms in the following papers are implemented: + * + * 1. SGD for linear SVM + * @ARTICLE{Zhang_SGD, + * author = "Tong Zhang", + * title = "{Solving Large Scale Linear Prediction Problems Using Stochastic Gradient Descent Algorithms}", + * booktitle = "{International Conference on Machine Learning}", + * year = 2004, + * } + * + * 2. SGD for nonlinear SVM + * @ARTICLE{Kivinen_SGD, + * author = "Jyrki Kivinen", + * title = "{Online Learning with Kernels}", + * booktitle = NIPS, + * number = 14, + * year = 2001, + * } + * + * @see svm.h + */ + +#ifndef U_SVM_OPT_SGD_H +#define U_SVM_OPT_SGD_H + +#include "fastlib/fastlib.h" + +// tolerance of sacale_w +const double SGD_SCALE_W_TOLERANCE = 1.0e-9; +// threshold that determines whether an alpha is a SV or not +const double SGD_ALPHA_ZERO = 1.0e-7; + + +template +class SGD { + FORBID_ACCIDENTAL_COPIES(SGD); + + public: + typedef TKernel Kernel; + + private: + int learner_typeid_; + + Kernel kernel_; + const Dataset_sl *dataset_; + index_t n_data_; /* number of data samples */ + index_t n_features_; /* # of features == # of row - 1, exclude the last row (for labels) */ + + Vector coef_; /* alpha*y, to be optimized */ + index_t n_alpha_; /* number of lagrangian multipliers in the dual */ + index_t n_sv_; /* number of support vectors */ + + index_t i_cache_, j_cache_; /* indices for the most recently cached kernel value */ + double cached_kernel_value_; /* cache */ + + ArrayList y_; /* list that stores "labels" */ + ArrayList w_; + double bias_; + double scale_w_; // the scale for w + + // parameters + double C_; // for SVM_C + double epsilon_; // for SVM_R + bool b_linear_; // whether it's a linear SVM + double lambda_; // regularization parameter. lambda = 1/(C*n_data) + index_t n_iter_; // number of iterations + index_t n_epochs_; // number of epochs + double accuracy_; // accuracy for stopping creterion + double eta_; // step length. eta = 1/(lambda*t) + double t_; + + ArrayList old_from_new_; // for generating a random sequence of training data + //ArrayList new_from_old_; // for generating a random sequence of training data + + double rho_;// for soft margin nonlinear SGD SVM + + public: + SGD() {} + ~SGD() {} + + /** + * Initialization for parameters + */ + void InitPara(int learner_typeid, ArrayList ¶m_) { + // init parameters + if (learner_typeid == 0) { // SVM_C + C_ = param_[0]; + b_linear_ = param_[2]>0.0 ? false: true; // whether it's a linear learner + n_epochs_ = (index_t)param_[3]; + n_iter_ = (index_t)param_[4]; + accuracy_ = param_[5]; + } + else if (learner_typeid == 1) { // SVM_R + } + } + + void Train(int learner_typeid, const Dataset_sl &dataset_in); + + Kernel& kernel() { + return kernel_; + } + + double Bias() const { + return bias_; + } + + void GetW(ArrayList &w_out) { + w_out.InitCopy(w_, w_.size()); + } + + double ScaleW() const { + return scale_w_; + } + + void GetSV(ArrayList &dataset_index, ArrayList &coef, ArrayList &sv_indicator); + + private: + + /** + * Loss functions + */ + double LossFunction_(int learner_typeid, double yy_hat) { + if (learner_typeid_ == 0) { // SVM_C + return HingeLoss_(yy_hat); + } + else if (learner_typeid_ == 1) { // SVM_R + return 0.0; // TODO + } + else + return HingeLoss_(yy_hat); + } + + /** + * Gradient of loss functions + */ + double LossFunctionGradient_(int learner_typeid, double yy_hat) { + if (learner_typeid_ == 0) { // SVM_C + return HingeLossGradient_(yy_hat); + } + else if (learner_typeid_ == 1) { // SVM_R + return 0.0; // TODO + } + else { + if (yy_hat < 1.0) + return 1.0; + else + return 0.0; + } + } + + /** + * Hinge Loss function + */ + double HingeLoss_(double yy_hat) { + if (yy_hat < 1.0) + return 1.0 - yy_hat; + else + return 0.0; + } + + /** + * Gradient of the Hinge Loss function + */ + double HingeLossGradient_(double yy_hat) { + if (yy_hat < 1.0) + return -1.0; + else + return 0.0; + } + + void LearnersInit_(int learner_typeid); + + int TrainIteration_(); + + double GetC_(index_t i) { + return C_; + } + + /** + * Calculate kernel values + */ + double CalcKernelValue_(index_t i, index_t j) { + // for SVM_R where n_alpha_==2*n_data_ + if (learner_typeid_ == 1) { + i = i >= n_data_ ? (i-n_data_) : i; + j = j >= n_data_ ? (j-n_data_) : j; + } + + // Check cache + //if (i == i_cache_ && j == j_cache_) { + // return cached_kernel_value_; + //} + + // Do Caching. Store the recently caculated kernel values. + //i_cache_ = i; + //j_cache_ = j; + + cached_kernel_value_ = kernel_.Eval((dataset_->x)[i], (dataset_->x)[j]); + + return cached_kernel_value_; + } +}; + + +/** + * Initialization according to different SVM learner types + * + * @param: learner type id + */ +template +void SGD::LearnersInit_(int learner_typeid) { + index_t i; + learner_typeid_ = learner_typeid; + rho_ = fx_param_double(NULL, "rho", 1.0); // specify the soft margin. default value 1.0: hard margin + + w_.Init(0); // all coefs of w_ are init as 0 + + if (learner_typeid_ == 0) { // SVM_C + if (b_linear_) { // linear SVM + coef_.Init(0); // not used, plain init + } + else { // nonlinear SVM + n_alpha_ = n_data_; + coef_.Init(n_alpha_); + coef_.SetZero(); + + //w_.Init(0); // not used, plain init + } + + y_.Init(n_data_); + for (i = 0; i < n_data_; i++) { + y_[i] = (dataset_->y)[i] > 0 ? 1 : -1; + } + } + else if (learner_typeid_ == 1) { // SVM_R + // TODO + n_alpha_ = 2 * n_data_; + + coef_.Init(n_alpha_); + coef_.SetZero(); + + y_.Init(n_alpha_); + for (i = 0; i < n_data_; i++) { + y_[i] = 1; // -> alpha_i + y_[i + n_data_] = -1; // -> alpha_i^* + } + } + else if (learner_typeid_ == 2) { // SVM_DE + // TODO + } +} + + +/** +* Steepest descent based SGD training for 2-classes +* +* @param: input 2-classes data matrix with labels (1,-1) in the last row +*/ +template +void SGD::Train(int learner_typeid, const Dataset_sl &dataset_in) { + index_t i, j, epo, ct; + // Load data + dataset_ = &dataset_in; + n_data_ = dataset_->n_points; + n_features_ = dataset_->n_features; + + /* general learner-independent initializations */ + if (n_epochs_ > 0) { // # of epochs provided, use it + n_iter_ = n_data_; + } + else { // # of epochs not provided, use n_iter_ to count iterations + n_epochs_ = 1; // not exactly one epoch, just use it for one loop + } + + DEBUG_ASSERT(C_ != 0); + lambda_ = 1.0/(C_*n_data_); + bias_ = 0.0; + + /* learners initialization */ + LearnersInit_(learner_typeid); + old_from_new_.Init(n_data_); + + index_t work_idx_old = 0; + + /* Begin SGD iterations */ + if (b_linear_) { // linear SVM, output: w, bias + double yt, yt_hat, yy_hat; + double sqrt_n = sqrt(n_data_); + double eta0 = sqrt_n / max(1.0, LossFunctionGradient_(learner_typeid, -sqrt_n)); // initial step length + t_ = 1.0 / (eta0 * lambda_); + scale_w_ = 1.0; + + printf("SGD_linear training begins...\n"); + for (epo = 0; epox)[work_idx_old]; + yt = y_[work_idx_old]; + //yt_hat = la::Dot(w_, xt) * scale_w_ + bias_; + yt_hat = SparseDot(w_, xt) * scale_w_ + bias_; + yy_hat = yt * yt_hat; + + // update w by Stochastic Gradient Descent: w_{t+1} = w_t - eta_t * grad_t = (1-eta*lambda) * w_t + eta * [yt*xt]^+ + double eta_grad = eta_ * LossFunctionGradient_(learner_typeid, yy_hat) * yt; + if (eta_grad != 0) { + // Note: moving w's scaling calculation w_t*(1-1/t) to the testing session is faster + SparseAddExpert(-eta_grad / scale_w_, xt, w_); + // update bias + bias_ -= eta_grad * 0.01; + } + + t_ += 1.0; + ct ++; + } + } // for epo + + // Calculate objective value; default: no calculation to save time + int objvalue = fx_param_int(NULL, "objvalue", 0); + if (objvalue > 0) { + double v = 0.0, hinge_loss = 0.0, loss_sum= 0.0; + + // primal objective value + for (i=0; i< n_data_; i++) { + NZ_entry *xt; + xt = (dataset_->x)[work_idx_old]; + //hinge_loss = 1- y_[i] * (scale_w_ * la::Dot(w_, xt) + bias_); + hinge_loss = 1- y_[i] * (scale_w_ * SparseDot(w_, xt) + bias_); + if (hinge_loss > 0) { + loss_sum += hinge_loss * C_; + } + } + + for (j=0; j wi_abs_max) { + wi_abs_max = wi_abs; + } + } + // round small w_i to 0 + index_t w_ct = 0; + double round_factor = fx_param_double(NULL, "round_factor", 1.0e32); + double round_thd = wi_abs_max / round_factor; + for (i=0; i round_thd ) { + w_ct ++; + } + else { + w_.Remove(i); + } + } + printf("%d out of %d features are non zero\n", w_ct, n_features_); + } + + } + else { // nonlinear SVM, output: coefs(i.e. alpha*y), bias + // it's more expensive to calc the accuracy than linear SVM, so we just use n_iter_ as stop criterion + double delta; + double yt, yt_hat, yy_hat; + double one_minus_eta_lambda; + index_t n_real_epo, n_data_res; + double kernel_value; + + Vector b_calc_kernel; + b_calc_kernel.Init(n_data_); + b_calc_kernel.SetAll(0); + + Vector coef_long; + n_iter_ = n_iter_ * n_epochs_; + coef_long.Init(n_iter_); + + // initial step length + //double sqrt_n = sqrt(n_data_); + //double eta0 = sqrt_n / max(1.0, LossFunctionGradient_(learner_typeid, -sqrt_n)); // initial step length + //t_ = 1.0 / (eta0 * lambda_); + //double eta0 = 1.0 / (2*lambda_); + + /* To mimic the online learning senario, we randomly permutate the training set, indexed by old_from_new_ */ + for (i=0; i 0) { + for (i=0; i= SGD_ALPHA_ZERO) { + b_calc_kernel[i%n_data_] = 1; + } + } + for (i=0; i0) { + kernel_value = CalcKernelValue_(old_from_new_[i], work_idx_old); + for (j=0; j0) { + kernel_value = CalcKernelValue_(old_from_new_[i], work_idx_old); + for (j=0; j= SGD_ALPHA_ZERO) { + yt_hat += coef_long[i] * CalcKernelValue_(old_from_new_[i], work_idx_old); + } + } + } + + yt_hat += bias_; + yy_hat = yt * yt_hat; + + // update step length + eta_ = 1.0 / (lambda_ * t_); + //eta_ = eta0 / sqrt(t_); + + if (ct >= 1) { + bias_ += coef_long[ct-1]; + } + + // update old coefs (for i +void SGD::GetSV(ArrayList &dataset_index, ArrayList &coef, ArrayList &sv_indicator) { + n_sv_ = 0; + if (learner_typeid_ == 0) {// SVM_C + for (index_t i = 0; i < n_data_; i++) { + if (fabs(coef_[i]) >= SGD_ALPHA_ZERO) { // support vectors found + coef.PushBack() = coef_[i]; + sv_indicator[dataset_index[i]] = true; + n_sv_++; + } + else { + coef.PushBack() = 0; + } + } + printf("Number of support vectors: %d.\n", n_sv_); + } + else if (learner_typeid_ == 1) {// SVM_R + // TODO + /* + for (index_t i = 0; i < n_data_; i++) { + double alpha_diff = -alpha_[i] + alpha_[i+n_data_]; // alpha_i^* - alpha_i + if (fabs(alpha_diff) >= SGD_ALPHA_ZERO) { // support vectors found + coef.PushBack() = alpha_diff; + sv_indicator[dataset_index[i]] = true; + n_sv_++; + } + else { + coef.PushBack() = 0; + } + } + */ + } +} + +#endif + +// ./regmin --learner_name=svm_c --mode=train_test --train_data=heart_scale --test_data=heart_scale --opt=sgd --kernel=linear --normalize=0 --objvalue=1 --c=0.1 --n_iter=4000 + +// ./regmin --learner_name=svm_c --mode=train_test --train_data=heart_scale --test_data=heart_scale --opt=sgd --kernel=gaussian --normalize=0 --objvalue=1 --sigma=1 --c=10 --n_iter=5000 + +// ./regmin --learner_name=svm_c --mode=train_test --train_data=real_sim --test_data=real_sim --opt=sgd --kernel=linear --normalize=0 --objvalue=0 --c=1 --n_iter=140000 diff --git a/fastlib/trunk/contrib/houyang/regmin/opt_smo.h b/fastlib/trunk/contrib/houyang/regmin/opt_smo.h new file mode 100644 index 0000000000..0dafb262dc --- /dev/null +++ b/fastlib/trunk/contrib/houyang/regmin/opt_smo.h @@ -0,0 +1,929 @@ +/** + * @author Hua Ouyang + * + * @file opt_smo.h + * + * This head file contains functions for performing Sequential Minimal Optimization (SMO) + * + * The algorithms in the following papers are implemented: + * + * 1. SMO and Working set selecting using 1st order expansion + * @ARTICLE{Platt_SMO, + * author = "J. C. Platt", + * title = "{Fast Training of Support Vector Machines using Sequential Minimal Optimization}", + * booktitle = "{Advances in Kernel Methods - Support Vector Learning}", + * year = 1999, + * publisher = "MIT Press" + * } + * + * 2. Shrinkng and Caching for SMO + * @ARTICLE{Joachims_SVMLIGHT, + * author = "T. Joachims", + * title = "{Making large-Scale SVM Learning Practical}", + * booktitle = "{Advances in Kernel Methods - Support Vector Learning}", + * year = 1999, + * publisher = "MIT Press" + * } + * + * 3. Working set selecting using 2nd order expansion + * @ARTICLE{Fan_JMLR, + * author = "R. Fan, P. Chen, C. Lin", + * title = "{Working Set Selection using Second Order Information for Training Support Vector Machines}", + * journal = "{Jornal of Machine Learning Research}", + * year = 2005 + * } + * + * @see svm.h + */ + +#ifndef U_SVM_OPT_SMO_H +#define U_SVM_OPT_SMO_H + +#include "fastlib/fastlib.h" +#include "regmin_data.h" + + +// maximum # of interations for SMO training +const index_t MAX_NUM_ITER_SMO = 10000000; +// after # of iterations to do shrinking +const index_t SMO_NUM_FOR_SHRINKING = 1000; +// threshold that determines whether need to do unshrinking +const double SMO_UNSHRINKING_FACTOR = 10; +// threshold that determines whether an alpha is a SV or not +const double SMO_ALPHA_ZERO = 1.0e-7; +// for indefinite kernels +const double TAU = 1e-12; + +const double SMO_ID_LOWER_BOUNDED = -1; +const double SMO_ID_UPPER_BOUNDED = 1; +const double SMO_ID_FREE = 0; + +template inline void swap(T& x, T& y) { T t=x; x=y; y=t; } + +template +class SMO { + FORBID_ACCIDENTAL_COPIES(SMO); + + public: + typedef TKernel Kernel; + + private: + int learner_typeid_; + int hinge_; // do L2-SVM or L1-SVM, default: L1 + + index_t ct_iter_; /* counter for the number of iterations */ + index_t ct_shrinking_; /* counter for doing shrinking */ + bool do_shrinking_; // 1(default): do shrinking after 1000 iterations; 0: don't do shrinking + + Kernel kernel_; + index_t n_data_; /* number of data samples */ + Dataset_sl *dataset_; /* alias for the input dataste */ + + Vector alpha_; /* the alphas, to be optimized */ + Vector alpha_status_; /* ID_LOWER_BOUND (-1), ID_UPPER_BOUND (1), ID_FREE (0) */ + index_t n_sv_; /* number of support vectors */ + + index_t n_alpha_; /* number of variables to be optimized */ + index_t n_active_; /* number of samples in the active set */ + ArrayList active_set_; /* list that stores the old indices of active alphas followed by inactive alphas. == old_from_new*/ + bool reconstructed_; /* indicator: where unshrinking has been carried out */ + index_t i_cache_, j_cache_; /* indices for the most recently cached kernel value */ + double cached_kernel_value_; /* cache */ + + ArrayList y_; /* list that stores "labels" */ + + double bias_; + + Vector grad_; /* gradient value */ + Vector grad_bar_; /* gradient value when treat un-upperbounded variables as 0: grad_bar_i==C\sum_{j:a_j=C} y_i y_j K_ij */ + + // parameters + int budget_; + double Cp_; // C_+, for SVM_C, y==1 + double Cn_; // C_-, for SVM_C, y==-1 + double C_; + double inv_two_C_; // 1/2C + double epsilon_; // for SVM_R + int wss_; // working set selection scheme, 1 for 1st order expansion; 2 for 2nd order expansion + index_t n_iter_; // number of iterations + double accuracy_; // accuracy for stopping creterion + double gap_; // for stopping criterion + + public: + SMO() {} + ~SMO() {} + + /** + * Initialization for parameters + */ + void InitPara(int learner_typeid, ArrayList ¶m_) { + // init parameters + budget_ = (int)param_[0]; + wss_ = (int) param_[4]; + hinge_ = (int) param_[3]; + n_iter_ = (index_t) param_[5]; + n_iter_ = n_iter_ < MAX_NUM_ITER_SMO ? n_iter_: MAX_NUM_ITER_SMO; + accuracy_ = param_[6]; + if (learner_typeid == 0) { // SVM_C + if (hinge_==2) { // L2-SVM: squated hinge loss + Cp_ = INFINITY; + Cn_ = INFINITY; + C_ = param_[1]; + inv_two_C_ = 1/(2*C_); + } + else { // L1-SVM + Cp_ = param_[1]; + Cn_ = param_[2]; + } + } + else if (learner_typeid == 1) { // SVM_R + Cp_ = param_[1]; + Cn_ = Cp_; + epsilon_ = param_[2]; + } + } + + void Train(int learner_typeid, Dataset_sl &dataset_in); + + Kernel& kernel() { + return kernel_; + } + + double Bias() const { + return bias_; + } + + void GetSV(ArrayList &dataset_index, ArrayList &coef, ArrayList &sv_indicator); + + private: + void LearnersInit_(int learner_typeid); + + int SMOIterations_(); + + void ReconstructGradient_(); + + bool TestShrink_(index_t i, double y_grad_max, double y_grad_min); + + void Shrinking_(); + + bool WorkingSetSelection_(index_t &i, index_t &j); + + void UpdateGradientAlpha_(index_t i, index_t j); + + void CalcBias_(); + + /* void GetVector_(index_t i, Vector *v) const { + datamatrix_.MakeColumnSubvector(i, 0, datamatrix_.n_rows()-1, v); + } + */ + + /** + * Instead of C, we use C_+ and C_- to handle unbalanced data + */ + double GetC_(index_t i) { + return (y_[i] > 0 ? Cp_ : Cn_); + } + + void UpdateAlphaStatus_(index_t i) { + if (alpha_[i] >= GetC_(i)) { + alpha_status_[i] = SMO_ID_UPPER_BOUNDED; + } + else if (alpha_[i] <= 0) { + alpha_status_[i] = SMO_ID_LOWER_BOUNDED; + } + else { // 0 < alpha_[i] < C + alpha_status_[i] = SMO_ID_FREE; + } + } + + bool IsUpperBounded(index_t i) { + return alpha_status_[i] == SMO_ID_UPPER_BOUNDED; + } + bool IsLowerBounded(index_t i) { + return alpha_status_[i] == SMO_ID_LOWER_BOUNDED; + } + + /** + * Calculate kernel values + */ + double CalcKernelValue_(index_t ii, index_t jj) { + // the indices have been swaped in the shrinking processes + index_t i = active_set_[ii]; // ii/jj: index in the new permuted set + index_t j = active_set_[jj]; // i/j: index in the old set + + // for SVM_R where n_alpha_==2*n_data_ + if (learner_typeid_ == 1) { + i = i >= n_data_ ? (i-n_data_) : i; + j = j >= n_data_ ? (j-n_data_) : j; + } + + // Check cache + //if (i == i_cache_ && j == j_cache_) { + // return cached_kernel_value_; + //} + + cached_kernel_value_ = kernel_.Eval((dataset_->x)[i], (dataset_->x)[j]); + + if (hinge_ == 2) { // for L2-SVM + if (i == j) { + cached_kernel_value_ = cached_kernel_value_ + inv_two_C_; + } + } + + return cached_kernel_value_; + } +}; + + +/** +* Reconstruct inactive elements of G from G_bar and free variables +* +* @param: learner type id +*/ +template +void SMO::ReconstructGradient_() { + index_t i, j; + if (n_active_ == n_alpha_) + return; + if (learner_typeid_ == 0) { // SVM_C + for (i=n_active_; i= n_data_ ? (i-n_data_) : i; + //grad_[j] = grad_bar_[j] + datamatrix_.get(datamatrix_.n_rows()-1, active_set_[j]) - epsilon_; // TODO + grad_[j] = grad_bar_[j] + (dataset_->y)[j] - epsilon_; // TODO + } + } + + for (i=0; i +bool SMO::TestShrink_(index_t i, double y_grad_max, double y_grad_min) { + if (IsUpperBounded(i)) { // alpha_[i] = C + if (y_[i] == 1) { + return (grad_[i] > y_grad_max); + } + else { // y_[i] == -1 + return (grad_[i] + y_grad_min > 0); // -grad_[i]y_grad_max + } + } + else + return false; +} + +/** + * Do Shrinking. Temporarily remove alphas (from the active set) that are + * unlikely to be selected in the working set, since they have reached their + * lower/upper bound. + * + */ +template +void SMO::Shrinking_() { + index_t t; + + // Find m(a) == y_grad_max(i\in I_up) and M(a) == y_grad_min(j\in I_down) + double y_grad_max = -INFINITY; + double y_grad_min = INFINITY; + for (t=0; t y_grad_max) { // y==1 + y_grad_max = grad_[t]; + } + } + else { // y[t] == -1 + if (!IsLowerBounded(t)) // t\in I_up, y==-1: y[t]alpha[t] < 0 + if (grad_[t] + y_grad_max < 0) { // y==-1... <=> -grad_[t] > y_grad_max + y_grad_max = -grad_[t]; + } + } + } + for (t=0; t 0 + if (grad_[t] < y_grad_min) { // y==1 + y_grad_min = grad_[t]; + } + } + else { // y[t] == -1 + if (!IsUpperBounded(t)) // t\in I_down, y==-1: y[t]alpha[t] > -C + if (grad_[t] + y_grad_min > 0) { // y==-1...<=> -grad_[t] < y_grad_min + y_grad_min = -grad_[t]; + } + } + } + + // Find the alpha to be shrunk + printf("Shrinking...\n"); + for (t=0; t t) { + if (!TestShrink_(n_active_, y_grad_max, y_grad_min)) { + swap(active_set_[t], active_set_[n_active_]); + swap(alpha_[t], alpha_[n_active_]); + swap(alpha_status_[t], alpha_status_[n_active_]); + swap(y_[t], y_[n_active_]); + swap(grad_[t], grad_[n_active_]); + swap(grad_bar_[t], grad_bar_[n_active_]); + break; + } + n_active_ --; + } + } + } + + double gap = y_grad_max - y_grad_min; + //printf("%d: gap:%f, n_active:%d\n", ct_iter_, gap, n_active_); + // do unshrinking for the first time when y_grad_max - y_grad_min <= SMO_UNSHRINKING_FACTOR * accuracy_ + if ( reconstructed_==false && gap <= SMO_UNSHRINKING_FACTOR * accuracy_ ) { + printf("Unshrinking...\n"); + // Unshrinking: put shrinked alphas back to active set + // 1.recover gradient + ReconstructGradient_(); + // 2.recover active status + for (t=n_alpha_-1; t>n_active_; t--) { + if (!TestShrink_(t, y_grad_max, y_grad_min)) { + while (n_active_ < t) { + if (TestShrink_(n_active_, y_grad_max, y_grad_min)) { + swap(active_set_[t], active_set_[n_active_]); + swap(alpha_[t], alpha_[n_active_]); + swap(alpha_status_[t], alpha_status_[n_active_]); + swap(y_[t], y_[n_active_]); + swap(grad_[t], grad_[n_active_]); + swap(grad_bar_[t], grad_bar_[n_active_]); + break; + } + n_active_ ++; + } + n_active_ ++; + } + } + reconstructed_ = true; // indicator: unshrinking has been carried out in this round + } + +} + + +/** + * Initialization according to different SVM learner types + * + * @param: learner type id + */ +template +void SMO::LearnersInit_(int learner_typeid) { + index_t i; + learner_typeid_ = learner_typeid; + + if (learner_typeid_ == 0) { // SVM_C + n_alpha_ = n_data_; + + alpha_.Init(n_alpha_); + alpha_.SetZero(); + + // initialize gradient + grad_.Init(n_alpha_); + grad_.SetAll(1.0); + + y_.Init(n_alpha_); + for (i = 0; i < n_alpha_; i++) { + y_[i] = (dataset_->y)[i] > 0 ? 1 : -1; + } + } + else if (learner_typeid_ == 1) { // SVM_R + n_alpha_ = 2 * n_data_; + + alpha_.Init(2 * n_alpha_); // TODO + alpha_.SetZero(); + + // initialize gradient + grad_.Init(n_alpha_); + y_.Init(n_alpha_); + for (i = 0; i < n_data_; i++) { + y_[i] = 1; // -> alpha_i + y_[i + n_data_] = -1; // -> alpha_i^* + grad_[i] = epsilon_ - (dataset_->y)[i]; + grad_[i + n_data_] = epsilon_ + (dataset_->y)[i]; + } + } + else if (learner_typeid_ == 2) { // SVM_DE + // TODO + } + +} + + +/** +* SMO training for 2-classes +* +* @param: input 2-classes data matrix with labels (1,-1) in the last row +*/ +template +void SMO::Train(int learner_typeid, Dataset_sl &dataset_in) { + index_t i; + // Load data + dataset_ = &dataset_in; + n_data_ = dataset_->n_points; + + // Learners initialization + LearnersInit_(learner_typeid); + + // General learner-independent initializations + budget_ = min(budget_, n_data_); + bias_ = 0.0; + n_sv_ = 0; + reconstructed_ = false; + i_cache_ = -1; j_cache_ = -1; + cached_kernel_value_ = INFINITY; + + n_active_ = n_alpha_; + active_set_.Init(n_alpha_); + for (i=0; i0 + grad_[i] = grad_[i] - y_[i] * y_[j] * alpha_[j] * CalcKernelValue_(i,j); + } + } + } + */ + + // initialize gradient_bar + grad_bar_.Init(n_alpha_); + grad_bar_.SetZero(); + + do_shrinking_ = fx_param_int(NULL, "shrink", 0); + ct_shrinking_ = min(n_data_, SMO_NUM_FOR_SHRINKING); + /* + if (do_shrinking_) { + for (i=0; i= C + grad_bar_[i] = grad_bar_[i] + GetC_(j) * y_[j] * CalcKernelValue_(i,j); + } + grad_bar_[i] = y_[i] * grad_bar_[i]; + } + } + */ + + printf("SMO initialization done!\n"); + + // Begin SMO iterations + ct_iter_ = 0; + + int stop_condition = 0; + while (1) { + //for(index_t i=0; i SMO terminates + if (stop_condition == 1) {// optimality reached + // Calculate the bias term + CalcBias_(); + printf("SMO terminates since the accuracy %f achieved!!! Number of iterations: %d.\n", accuracy_, ct_iter_); + break; + } + else if (stop_condition == 2) {// max num of iterations exceeded + // Calculate the bias term + CalcBias_(); + fprintf(stderr, "SMO terminates since the number of iterations %d exceeded !!! Gap: %f.\n", n_iter_, gap_); + break; + } + } +} + +/** +* SMO training iterations +* +* @return: stopping condition id +*/ +template +int SMO::SMOIterations_() { + ct_iter_ ++; + index_t i,j; + if (WorkingSetSelection_(i,j) == true) { + if (!do_shrinking_) { // no shrinking, optimality reached + return 1; + } + else { // shrinking, need to check whether optimality really reached + ReconstructGradient_(); // restore the inactive alphas and reconstruct gradients + n_active_ = n_alpha_; + if (WorkingSetSelection_(i,j) == true) { // optimality reached + return 1; + } + else { + ct_shrinking_ = 1; // do shrinking in the next iteration + return 0; + } + } + } + else if (ct_iter_ >= n_iter_) { // number of iterations exceeded + if (!do_shrinking_) { // no shrinking, optimality reached + return 2; + } + else if ( ct_iter_ >= min(n_data_, SMO_NUM_FOR_SHRINKING) ) { // shrinking has been carried out, need to calculate the true gap + ReconstructGradient_(); // restore the inactive alphas and reconstruct gradients + n_active_ = n_alpha_; + WorkingSetSelection_(i,j); + return 2; + } + else { + return 2; + } + } + else{ // update gradient, alphas and bias term, and continue iterations + UpdateGradientAlpha_(i, j); + return 0; + } +} + +/** +* Try to find a working set (i,j). Both 1st(default) and 2nd order approximations of +* the objective function Z(\alpha+\lambda u_ij)-Z(\alpha) are implemented. +* +* @param: reference to working set (i, j) +* +* @return: working set (i, j); indicator of whether the optimal solution is reached (true:reached) +*/ +template +bool SMO::WorkingSetSelection_(index_t &out_i, index_t &out_j) { + double y_grad_max = -INFINITY; + double y_grad_min = INFINITY; + int idx_i = -1; + int idx_j = -1; + + // Find i using maximal violating pair scheme + index_t t; + for (t=0; t y_grad_max) { // y==1 + y_grad_max = grad_[t]; + idx_i = t; + } + } + else { // y[t] == -1 + if (!IsLowerBounded(t)) // t\in I_up, y==-1: y[t]alpha[t] < 0 + if (grad_[t] + y_grad_max < 0) { // y==-1... <=> -grad_[t] > y_grad_max + y_grad_max = -grad_[t]; + idx_i = t; + } + } + } + out_i = idx_i; // i found + + /* Find j using maximal violating pair scheme (1st order approximation) */ + if (wss_ == 1) { + for (t=0; t 0 + if (grad_[t] < y_grad_min) { // y==1 + y_grad_min = grad_[t]; + idx_j = t; + } + } + else { // y[t] == -1 + if (!IsUpperBounded(t)) // t\in I_down, y==-1: y[t]alpha[t] > -C + if (grad_[t] + y_grad_min > 0) { // y==-1...<=> -grad_[t] < y_grad_min + y_grad_min = -grad_[t]; + idx_j = t; + } + } + } + out_j = idx_j; // j found + } + /* Find j using 2nd order working set selection scheme; need to calc kernels, but faster convergence */ + else if (wss_ == 2) { + double K_ii = CalcKernelValue_(out_i, out_i); + double opt_gain_max = -INFINITY; + double grad_diff; + double quad_kernel; + double opt_gain = -INFINITY; + for (t=0; t 0 + // calculate y_grad_min for Stopping Criterion + if (grad_[t] < y_grad_min) // y==1 + y_grad_min = grad_[t]; + // find j + grad_diff = y_grad_max - grad_[t]; // max(y_i*grad_i) - y_t*grad_t + if (grad_diff > 0) { + quad_kernel = K_ii + K_tt - 2 * K_it; + if (quad_kernel > 0) // for positive definite kernels + opt_gain = ( grad_diff * grad_diff ) / quad_kernel; // actually ../2*quad_kernel + else // handle non-positive definite kernels + opt_gain = ( grad_diff * grad_diff ) / TAU; + // find max(opt_gain) + if (opt_gain > opt_gain_max) { + idx_j = t; + opt_gain_max = opt_gain; + } + } + } + } + else { // y[t] == -1 + if (!IsUpperBounded(t)) {// t\in I_down, y==-1: y[t]alpha[t] > -C + // calculate y_grad_min for Stopping Criterion + if (grad_[t] + y_grad_min > 0) // y==-1, -grad_[t] < y_grad_min + y_grad_min = -grad_[t]; + // find j + grad_diff = y_grad_max + grad_[t]; // max(y_i*grad_i) - y_t*grad_t + if (grad_diff > 0) { + quad_kernel = K_ii + K_tt - 2 * K_it; + if (quad_kernel > 0) // for positive definite kernels + opt_gain = ( grad_diff * grad_diff ) / quad_kernel; // actually ../2*quad_kernel + else // handle non-positive definite kernels + opt_gain = ( grad_diff * grad_diff ) / TAU; + // find max(opt_gain) + if (opt_gain > opt_gain_max) { + idx_j = t; + opt_gain_max = opt_gain; + } + } + } + } + } + } + out_j = idx_j; // j found + + //printf("y_i=%d, y_j=%d\n", y_[out_i], y_[out_j]); + //printf("a_i=%f, a_j=%f\n", alpha_[out_i], alpha_[out_j]); + + // Stopping Criterion check + //printf("ct_iter:%d, accu:%f\n", ct_iter_, y_grad_max - y_grad_min); + gap_ = y_grad_max - y_grad_min; + //printf("%d: gap=%f\n", ct_iter_, gap_); + if (gap_ <= accuracy_) { + return true; // optimality reached + } + + return false; +} + +/** +* Search direction; Update gradient, alphas and bias term +* +* @param: a working set (i,j) found by working set selection +* +*/ +template +void SMO::UpdateGradientAlpha_(index_t i, index_t j) { + index_t t; + + double a_i = alpha_[i]; // old alphas + double a_j = alpha_[j]; + int y_i = y_[i]; + int y_j = y_[j]; + double C_i = GetC_(i); // can be Cp (for y==1) or Cn (for y==-1) + double C_j = GetC_(j); + + // cached kernel values + double K_ii, K_ij, K_jj; + K_ii = CalcKernelValue_(i, i); + K_ij = CalcKernelValue_(i, j); + K_jj = CalcKernelValue_(j, j); + + double first_order_diff = y_i * grad_[i] - y_j * grad_[j]; + double second_order_diff = K_ii + K_jj - 2 * K_ij; + if (second_order_diff <= 0) // handle non-positive definite kernels + second_order_diff = TAU; + double lambda = first_order_diff / second_order_diff; // step size + + //printf("step size=%f\n", lambda); + + // Update alphas + alpha_[i] = a_i + y_i * lambda; + alpha_[j] = a_j - y_j * lambda; + + // Handle bounds for updated alphas + if (y_i != y_j) { + double alpha_old_diff = a_i - a_j; + if (alpha_old_diff > 0) { + if (alpha_[j] < 0) { + alpha_[j] = 0; + alpha_[i] = alpha_old_diff; + } + } + else { // alpha_old_diff <= 0 + if (alpha_[i] < 0) { + alpha_[i] = 0; + alpha_[j] = - alpha_old_diff; + } + } + if (alpha_old_diff > C_i - C_j) { + if (alpha_[i] > C_i) { + alpha_[i] = C_i; + alpha_[j] = C_i - alpha_old_diff; + } + } + else { + if (alpha_[j] > C_j) { + alpha_[j] = C_j; + alpha_[i] = C_j + alpha_old_diff; + } + } + } + else { // y_i == y_j + double alpha_old_sum = a_i + a_j; + if (alpha_old_sum > C_i) { + if (alpha_[i] > C_i) { + alpha_[i] = C_i; + alpha_[j] = alpha_old_sum - C_i; + } + } + else { + if (alpha_[j] < 0) { + alpha_[j] = 0; + alpha_[i] = alpha_old_sum; + } + } + if (alpha_old_sum > C_j) { + if (alpha_[j] > C_j) { + alpha_[j] = C_j; + alpha_[i] = alpha_old_sum - C_j; + } + } + else { + if (alpha_[i] < 0) { + alpha_[i] = 0; + alpha_[j] = alpha_old_sum; + } + } + } + + // Update gradient + double diff_i = alpha_[i] - a_i; + double diff_j = alpha_[j] - a_j; + for (t=0; t= C + if(ub_i) // old_alpha_i >= C, new_alpha_i < C + for(t=0; t= C + for(t=0; t= C, new_alpha_j < C + for(t=0; t= C + for(t=0; t +void SMO::CalcBias_() { + double b; + index_t n_free_alpha = 0; + double ub = INFINITY, lb = -INFINITY, sum_free_yg = 0.0; + + for (index_t i=0; i= C + if(y_[i] == 1) + lb = max(lb, yg); + else + ub = min(ub, yg); + } + else if (IsLowerBounded(i)) { // bounded: alpha_i <= 0 + if(y_[i] == -1) + lb = max(lb, yg); + else + ub = min(ub, yg); + } + else { // free: 0< alpha_i 0) + b = sum_free_yg / n_free_alpha; + else + b = (ub + lb) / 2; + + bias_ = b; +} + +/* Get SVM results:coefficients, number and indecies of SVs +* +* @param: sample indices of the training (sub)set in the total training set +* @param: support vector coefficients: alpha*y +* @param: bool indicators FOR THE TRAINING SET: is/isn't a support vector +* +*/ +template +void SMO::GetSV(ArrayList &dataset_index, ArrayList &coef, ArrayList &sv_indicator) { + ArrayList new_from_old; // it's used to retrieve the permuted new index from old index + new_from_old.Init(n_alpha_); + for (index_t i = 0; i < n_alpha_; i++) { + new_from_old[active_set_[i]] = i; + } + if (learner_typeid_ == 0) {// SVM_C + for (index_t ii = 0; ii < n_data_; ii++) { + index_t i = new_from_old[ii]; // retrive the index of permuted vector + if (alpha_[i] >= SMO_ALPHA_ZERO) { // support vectors found + //printf("%f\n", alpha_[i] * y_[i]); + coef.PushBack() = alpha_[i] * y_[i]; + sv_indicator[dataset_index[ii]] = true; + n_sv_++; + } + else { + coef.PushBack() = 0; + } + } + printf("Number of SVs: %d\n", n_sv_); + } + else if (learner_typeid_ == 1) {// SVM_R + for (index_t ii = 0; ii < n_data_; ii++) { + index_t i = new_from_old[ii]; // retrive the index of permuted vector + index_t iplusn = new_from_old[ii+n_data_]; + double alpha_diff = -alpha_[i] + alpha_[iplusn]; // alpha_i^* - alpha_i + if (fabs(alpha_diff) >= SMO_ALPHA_ZERO) { // support vectors found + coef.PushBack() = alpha_diff; + sv_indicator[dataset_index[ii]] = true; + n_sv_++; + } + else { + coef.PushBack() = 0; + } + } + } +} + +#endif + +// ./regmin --learner_name=svm_c --mode=train_test --train_data=heart_scale --test_data=heart_scale --opt=smo --kernel=linear --normalize=0 --objvalue=1 --c=0.1 --wss=2 --shrink=1 --n_iter=300 + +// ./regmin --learner_name=svm_c --mode=train_test --train_data=heart_scale --test_data=heart_scale --opt=smo --kernel=gaussian --normalize=0 --objvalue=1 --wss=2 --shrink=1 --c=10 --sigma=1 --n_iter=300 diff --git a/fastlib/trunk/contrib/houyang/regmin/opt_tgd.h b/fastlib/trunk/contrib/houyang/regmin/opt_tgd.h new file mode 100644 index 0000000000..cc6fe04e83 --- /dev/null +++ b/fastlib/trunk/contrib/houyang/regmin/opt_tgd.h @@ -0,0 +1,390 @@ +/** + * @author Hua Ouyang + * + * @file opt_tgd.h + * + * This head file contains functions for performing L1-regularized linear loss optimization, using Truncated Gradient Descent + * + * + * @see svm.h + */ + +#ifndef U_SVM_OPT_TGD_H +#define U_SVM_OPT_TGD_H + +#include "fastlib/fastlib.h" + +const double TGD_ZERO = 1.0e-30; + +template +class TGD { + FORBID_ACCIDENTAL_COPIES(TGD); + + public: + typedef TKernel Kernel; + + private: + int learner_typeid_; + + Kernel kernel_; + const Dataset_sl *dataset_; + index_t n_data_; /* number of data samples */ + index_t n_features_; /* # of features == # of row - 1, exclude the last row (for labels) */ + // index_t n_features_bias_; /* # of features + 1 , [x, 1], for the bias term */ + + index_t n_sv_; /* number of support vectors */ + + index_t w_nnz_; + + double round_thd_; + + ArrayList y_; /* list that stores "labels" */ + + ArrayList w_; /* the slope of the decision hyperplane, including bias: [w, b] */ + + double scale_w_; // the scale for w + + // parameters + double C_; // \|w\|_1^1 \leq C + double g_; // regularization parameter in Langford's paper. g = 1/C + index_t n_iter_; // number of iterations + index_t n_epochs_; // number of epochs + double accuracy_; // accuracy for stopping creterion + double eta_; // step length. eta = 1/sqrt(t) + double t_; + index_t k_; // perform truncation every k iterations + + ArrayList old_from_new_; // for generating a random sequence of training data + + public: + TGD() {} + ~TGD() {} + + /** + * Initialization for parameters + */ + void InitPara(int learner_typeid, ArrayList ¶m_) { + // init parameters + if (learner_typeid == 0) { // SVM_C + C_ = param_[0]; + n_epochs_ = (index_t)param_[2]; + n_iter_ = (index_t)param_[3]; + accuracy_ = param_[4]; + } + else if (learner_typeid == 1) { // SVM_R + } + } + + void Train(int learner_typeid, Dataset_sl &dataset_in); + + Kernel& kernel() { + return kernel_; + } + + void GetW(ArrayList &w_out) { + index_t w_size = w_.size(); + index_t ct_nz = 0; + w_out.Init(w_nnz_); + for (index_t i=0; i= round_thd_ ) { + w_out[ct_nz].index = w_[i].index; + w_out[ct_nz].value = w_[i].value; + ct_nz ++; + } + } + } + + double ScaleW() const { + return scale_w_; + } + + //void GetSV(ArrayList &dataset_index, ArrayList &coef, ArrayList &sv_indicator); + + private: + /** + * Loss functions + */ + double LossFunction_(int learner_typeid, double yy_hat) { + if (learner_typeid_ == 0) { // SVM_C + return HingeLoss_(yy_hat); + } + else if (learner_typeid_ == 1) { // SVM_R + return 0.0; // TODO + } + else + return HingeLoss_(yy_hat); + } + + /** + * Gradient of loss functions + */ + double LossFunctionGradient_(int learner_typeid, double yy_hat) { + if (learner_typeid_ == 0) { // SVM_C + //return HingeLossGradient_(yy_hat); + return LogisticLossGradient_(yy_hat); + } + else if (learner_typeid_ == 1) { // SVM_R + return 0.0; // TODO + } + else { + if (yy_hat < 1.0) + return 1.0; + else + return 0.0; + } + } + + /** + * Hinge Loss function + */ + double HingeLoss_(double yy_hat) { + if (yy_hat < 1.0) + return 1.0 - yy_hat; + else + return 0.0; + } + + /** + * Gradient of the Hinge Loss function + */ + double HingeLossGradient_(double yy_hat) { + if (yy_hat < 1.0) + return -1.0; + else + return 0.0; + } + + /** + * Gradient of the Logistic Loss function + */ + double LogisticLossGradient_(double yy_hat) { + double tmp = exp(-yy_hat); + return -tmp/(1+tmp); + } + + void LearnersInit_(int learner_typeid); + + int TrainIteration_(); + + double GetC_(index_t i) { + return C_; + } + +}; + + +/** + * Initialization according to different SVM learner types + * + * @param: learner type id + */ +template +void TGD::LearnersInit_(int learner_typeid) { + index_t i; + + learner_typeid_ = learner_typeid; + + // init w, w+, w- + w_.Init(n_features_); + for (i=0; iy)[i] > 0 ? 1 : -1; + } +} + + +/** +* L1-regularization training for 2-classes +* +* @param: input 2-classes data matrix with labels (1,-1) in the last row +*/ +template +void TGD::Train(int learner_typeid, Dataset_sl &dataset_in) { + index_t i, j, epo, ct, ct_k; + + index_t total_n_iter; + + k_ = fx_param_int(NULL, "k", 1); + g_ = 1.0 / C_; + + /* general learner-independent initializations */ + dataset_ = &dataset_in; + n_data_ = dataset_->n_points; + n_features_ = dataset_->n_features; + + if (n_epochs_ > 0) { // # of epochs provided, use it + n_iter_ = n_data_; + total_n_iter = n_iter_ * n_epochs_; + } + else { // # of epochs not provided, use n_iter_ to count iterations + n_epochs_ = 1; // not exactly one epoch, just use it for one loop + total_n_iter = n_iter_; + } + + DEBUG_ASSERT(C_ != 0); + + /* learners initialization */ + LearnersInit_(learner_typeid); + old_from_new_.Init(n_data_); + + index_t work_idx_old = 0; + + /* Begin training iterations */ + double yt, yt_hat, yy_hat; + //double sqrt_n = sqrt(n_data_); + //double eta0 = sqrt_n / max(1.0, LossFunctionGradient_(learner_typeid, -sqrt_n)); // initial step length + //double eta_grad = 0; + //t_ = 1.0 / (eta0 * lambda_); + t_ = 1.0; + scale_w_ = 1.0; + + ct_k = 0; + printf("TGD training begins...\n"); + for (epo = 0; epox)[work_idx_old]; + yt = y_[work_idx_old]; + yt_hat = SparseDot(w_, xt) * scale_w_; + yy_hat = yt * yt_hat; + + + double grad_tmp = yt * LossFunctionGradient_(learner_typeid, yy_hat); + + //printf("epo:%d, ct:%d, w_size=%d, y_grad:%f\n", epo, ct, w_.size(), grad_tmp); + for (i=0; iindex) { + if (w_[i].value >= -g_ && w_[i].value <= g_) { + w_[i].value = 0; + } + // else w_[i].value remains unchanged + } + else { // i== xt->index + double db_tmp = w_[i].value - eta_ * grad_tmp * xt->value; + if (db_tmp >= 0.0 && db_tmp <= g_) { + if (!(ct_k % k_)) { + db_tmp -= eta_ * k_ * g_; + } + if (db_tmp >=0.0) { + w_[i].value = db_tmp; + } + else { + w_[i].value = 0.0; + } + } + else if (db_tmp >= -g_ && db_tmp <= 0.0) { + if (!(ct_k % k_)) { + db_tmp += eta_ * k_ * g_; + } + if (db_tmp <=0.0) { + w_[i].value = db_tmp; + } + else { + w_[i].value = 0.0; + } + } + else { + w_[i].value = db_tmp; + } + ++xt; + } + } + + /* + printf("epo:%d, ct:%d\n", epo, ct); + for (i=0; i wi_abs_max) { + wi_abs_max = wi_abs; + } + } + + // round small w_i to 0 + index_t w_ct = 0; + double round_factor = fx_param_double(NULL, "round_factor", 1.0e32); + double round_thd = wi_abs_max / round_factor; + for (i=0; i round_thd ) { + w_ct ++; + //printf("w_dim:%d, w_value:%lf\n", i, w_[i]); + } + else { + w_.Remove(i); + } + } + printf("%d out of %d features are non zero\n", w_ct, n_features_); + */ + + // rounding w + index_t w_size = w_.size(); + round_thd_ = fx_param_double(NULL, "thd", 1.0e-5); + w_nnz_ = 0; + for (ct=0; ct= round_thd_) { // TODO: thresholding + w_nnz_++; + } + } + + printf("%d out of %d features are non zero. NZ rate:%f\n", w_nnz_, n_features_, (double)(w_nnz_)/(double)n_features_); + + // Calculate objective value; default: no calculation to save time + int objvalue = fx_param_int(NULL, "objvalue", 0); + if (objvalue > 0) { + double hinge_loss = 0.0, loss_sum= 0.0, w_sum = 0.0; + + // primal objective value + for (i=0; ix)[i]); + if (hinge_loss > 0) { + loss_sum += hinge_loss; + } + } + for (i=0; i +#include "regmin.h" + +char *line = NULL; +index_t max_line_length; // buffer size for read a line + +Dataset_sl train_set; +NZ_entry *train_nz_pool; // a pool for all non-zero entries in the training set +Dataset_sl test_set; +NZ_entry *test_nz_pool; // a pool for all non-zero entries in the testing set + + +static char *ReadLine(FILE *fp_in) { + index_t length; + if ( fgets(line, max_line_length, fp_in)==NULL ) { + return NULL; + } + while ( strrchr(line,'\n')==NULL ) { + max_line_length *= 2; + line = (char *) realloc(line, max_line_length); + length = (index_t) strlen(line); + if ( fgets(line+length, max_line_length-length, fp_in)==NULL ) { + break; + } + } + return line; +} + +int ReadData(Dataset_sl &dataset, struct NZ_entry *nz_pool, FILE *fp) { + index_t max_index, inst_max_index, i, j; + char *endptr, *label, *index, *value; + + max_index = 0; + j = 0; + for (i=0; i max_index) { + max_index = inst_max_index; + } + nz_pool[j++].index = -1; // an indicator of the end of a data point + } + dataset.n_features = max_index + 1; + return 1; +} + +/** + * Initialize sparse training dataset from a file + * + * @param: the training set + * @param: the training filename + */ +int InitTrainsetFromFile(String param) { + if (fx_param_exists(NULL, param)) { + String train_filename = fx_param_str_req(NULL, param); + FILE *fp = fopen(train_filename, "r"); + if (fp == NULL) { + fprintf(stderr, "Cannot open the specified training file!!!\n"); + return 0; + } + else { + index_t num_nz_entries = 0; + train_set.n_points= 0; + + // count # of data points, # of non-zero entries + while ( ReadLine(fp)!= NULL ) { + // skip the lable + char *p = strtok(line, " \t"); + // count # of features + while (1) { + p = strtok(NULL, " \t"); + if (p == NULL || *p == '\n') { + break; + } + num_nz_entries ++; + } + num_nz_entries ++; // add an indicator for the end of this data point + train_set.n_points ++; + } + rewind(fp); + + train_set.y = Malloc(double, train_set.n_points); + train_set.x = Malloc(struct NZ_entry *, train_set.n_points); + train_nz_pool = Malloc(struct NZ_entry, num_nz_entries); + + if (ReadData(train_set, train_nz_pool, fp)) { + fclose(fp); + return 1; + } + else { + free(train_set.y); + free(train_set.x); + free(train_nz_pool); + fclose(fp); + fprintf(stderr, "Errors in training file format!!!\n"); + return 0; + } + } + } + else { + fprintf(stderr, "No training filename specified !!!\n"); + return 0; + } +} + +/** + * Initialize sparse testing dataset from a file + * + * @param: the testing set + * @param: the testing filename + */ +int InitTestsetFromFile(String param) { + if (fx_param_exists(NULL, param)) { + String test_filename = fx_param_str_req(NULL, param); + FILE *fp = fopen(test_filename, "r"); + if (fp == NULL) { + fprintf(stderr, "Cannot open the specified testing file!!!\n"); + return 0; + } + else { + index_t num_nz_entries = 0; + test_set.n_points= 0; + + // count # of data points, # of non-zero entries + while ( ReadLine(fp)!= NULL ) { + // skip the lable + char *p = strtok(line, " \t"); + // count # of features + while (1) { + p = strtok(NULL, " \t"); + if (p == NULL || *p == '\n') { + break; + } + num_nz_entries ++; + } + num_nz_entries ++; // add an indicator for the end of this data point + test_set.n_points ++; + } + rewind(fp); + + test_set.y = Malloc(double, test_set.n_points); + test_set.x = Malloc(struct NZ_entry *, test_set.n_points); + test_nz_pool = Malloc(struct NZ_entry, num_nz_entries); + + if (ReadData(test_set, test_nz_pool, fp)) { + fclose(fp); + return 1; + } + else { + free(test_set.y); + free(test_set.x); + free(test_nz_pool); + fclose(fp); + fprintf(stderr, "Errors in testing file format!!!\n"); + return 0; + } + } + } + else { + fprintf(stderr, "No testing filename specified!!!\n"); + return 0; + } +} + + +/** +* Multiclass SVM classification/ SVM regression - Main function +* +* @param: argc +* @param: argv +*/ +int main(int argc, char *argv[]) { + fx_init(argc, argv, NULL); + + srand(time(NULL)); + + String mode = fx_param_str_req(NULL, "mode"); + String kernel = fx_param_str_req(NULL, "kernel"); + String learner_name = fx_param_str_req(NULL,"learner_name"); + int learner_typeid; + + if (learner_name == "svm_c") { // Support Vector Classfication + learner_typeid = 0; + } + else if (learner_name == "svm_r") { // Support Vector Regression + learner_typeid = 1; + } + else if (learner_name == "svm_q") { // Support Vector Quantile Estimation + learner_typeid = 2; + } + else { + fprintf(stderr, "Unknown support vector learner name!!!\n"); + return 0; + } + + max_line_length = 1024; + line = Malloc(char, max_line_length); + + /* Training Mode, need training data | Training + Testing(online) Mode, need training data + testing data */ + if (mode=="train" || mode=="train_test"){ + fprintf(stderr, "SVM Training... \n"); + + /* Load training data */ + if (InitTrainsetFromFile("train_data") == 0) { + free(line); + exit(1); + } + + /* Begin SVM Training | Training and Testing */ + datanode *svm_module = fx_submodule(fx_root, "svm"); + if (kernel == "linear") { + SVM svm; + svm.InitTrain(learner_typeid, train_set, svm_module); + + /* training and testing, thus no need to load model from file */ + if (mode=="train_test"){ + fprintf(stderr, "SVM Predicting... \n"); + /* Load testing data */ + if (InitTestsetFromFile("test_data") == 0) { + free(line); + exit(1); + } + svm.BatchPredict(learner_typeid, test_set, "predicted_values"); + free(test_set.y); + free(test_set.x); + free(test_nz_pool); + } + + free(train_set.y); + free(train_set.x); + free(train_nz_pool); + } + else if (kernel == "gaussian") { + SVM svm; + svm.InitTrain(learner_typeid, train_set, svm_module); + + /* training and testing, thus no need to load model from file */ + if (mode=="train_test"){ + fprintf(stderr, "SVM Predicting... \n"); + /* Load testing data */ + if (InitTestsetFromFile("test_data") == 0) { + free(line); + exit(1); + } + svm.BatchPredict(learner_typeid, test_set, "predicted_values"); + free(test_set.y); + free(test_set.x); + free(test_nz_pool); + } + + free(train_set.y); + free(train_set.x); + free(train_nz_pool); + } + } + /* Testing(offline) Mode, need loading model file and testing data */ + else if (mode=="test") { + fprintf(stderr, "SVM Predicting... \n"); + + /* Load testing data */ + if (InitTestsetFromFile("test_data") == 0) { + free(line); + exit(1); + } + + /* Begin Prediction */ + datanode *svm_module = fx_submodule(fx_root, "svm"); + + if (kernel == "linear") { + SVM svm; + svm.Init(learner_typeid, test_set, svm_module); + svm.LoadModelBatchPredict(learner_typeid, test_set, "svm_model", "predicted_values"); + free(test_set.y); + free(test_set.x); + free(test_nz_pool); + } + else if (kernel == "gaussian") { + SVM svm; + svm.Init(learner_typeid, test_set, svm_module); + svm.LoadModelBatchPredict(learner_typeid, test_set, "svm_model", "predicted_values"); + free(test_set.y); + free(test_set.x); + free(test_nz_pool); + } + } + free(line); + fx_done(NULL); +} + diff --git a/fastlib/trunk/contrib/houyang/regmin/regmin.h b/fastlib/trunk/contrib/houyang/regmin/regmin.h new file mode 100644 index 0000000000..15e9bda370 --- /dev/null +++ b/fastlib/trunk/contrib/houyang/regmin/regmin.h @@ -0,0 +1,1100 @@ +/** + * @author Hua Ouyang + * + * @file regmin.h + * + * This head file contains functions for training and prediction of + * regulazied risk minimization problems. + * Supported learner type:SVM_C, SVM_R, SVM_Q + * + * @see opt_smo.h + * @see opt_sgd.h + * @see opt_md.h + * @see opt_tgd.h + */ + +#ifndef U_REGMIN_H +#define U_REGMIN_H + +#include "fastlib/fastlib.h" +#include "regmin_data.h" +#include "opt_smo.h" +#include "opt_sgd.h" +#include "opt_md.h" +#include "opt_tgd.h" + +#include + + +/** +* Class for SVM +*/ +template +class SVM { + + private: + /** + * Type id of the SVM learner: + * 0:SVM Classification (svm_c); + * 1:SVM Regression (svm_r); + * 2:SVM quantile estimation (svm_q); + * Developers may add more learner types if necessary + */ + int learner_typeid_; + // Optimization method: smo, lasvm, sgd, tgd, cd, pegasos, rivanov, hcy, fw, mfw, sfw, smd, mfw, par, sparsereg + String opt_method_; + /* array of models for storage of the 2-class(binary) classifiers + Need to train num_classes_*(num_classes_-1)/2 binary models */ + struct SVM_MODELS { + /* bias term in each binary model */ + double bias_; + /* all coefficients (alpha*y) of the binary dataset, not necessarily thoes of SVs */ + ArrayList coef_; + /* the slope w */ + //Vector w_; + ArrayList w_; + /* scale for w*/ + double scale_w_; // Use it if w's scaling is not done in training session + }; + ArrayList models_; + + /* list of labels, double type, but may be converted to integers. + e.g. [0.0,1.0,2.0] for a 3-class dataset */ + ArrayList train_labels_list_; + /* array of label indices, after grouping. e.g. [c1[0,5,6,7,10,13,17],c2[1,2,4,8,9],c3[...]]*/ + ArrayList train_labels_index_; + /* counted number of label for each class. e.g. [7,5,8]*/ + ArrayList train_labels_ct_; + /* start positions of each classes in the training label list. e.g. [0,7,12] */ + ArrayList train_labels_startpos_; + + /* total set of support vectors and their coefficients */ + struct NZ_entry **sv_entries_; // store the SVs + struct NZ_entry *sv_nz_pool_; // store the nonzero entries of SVs for testing + Matrix sv_coef_; + ArrayList trainset_sv_indicator_; + + /* total number of support vectors */ + index_t total_num_sv_; + /* support vector list to store the indices (in the training set) of support vectors */ + ArrayList sv_index_; + /* start positions of each class of support vectors, in the support vector list */ + ArrayList sv_list_startpos_; + /* counted number of support vectors for each class */ + ArrayList sv_list_ct_; + + /* SVM parameters */ + struct PARAMETERS { + TKernel kernel_; + String kernelname_; + int kerneltypeid_; + int b_; + double C_; + // for SVM_C of unbalanced data + double Cp_; // C for y==1 + double Cn_; // C for y==-1 + // for nu-SVM + double nu_; + // for SVM_R + double epsilon_; + // working set selection scheme of SMO, 1 for 1st order expansion; 2 for 2nd order expansion + double wss_; + // whether do L1-SVM (1) or L2-SVM (2) + int hinge_; + // accuracy for the optimization stopping creterion + double accuracy_; + // number of iterations + index_t n_iter_; + // number of epochs for stochastic algorithms + index_t n_epochs_; + }; + PARAMETERS param_; + + /* number of data samples */ + index_t n_data_; + /* number of classes in the training set */ + int num_classes_; + /* number of binary models to be trained, i.e. num_classes_*(num_classes_-1)/2 */ + int num_models_; + int num_features_; + + index_t max_line_length_;; + char *line_; + + public: + typedef TKernel Kernel; + class SMO; + class SGD; + class MD; + class TGD; + + void Init(int learner_typeid, Dataset_sl& dataset, datanode *module); + void InitTrain(int learner_typeid, Dataset_sl& dataset, datanode *module); + void GetLabels(Dataset_sl& dataset, + ArrayList &labels_list, + ArrayList &labels_index, + ArrayList &labels_ct, + ArrayList &labels_startpos); + + double Predict(int learner_typeid, NZ_entry *test_vec); + void BatchPredict(int learner_typeid, Dataset_sl& testset, String predictedvalue_filename); + void LoadModelBatchPredict(int learner_typeid, Dataset_sl& testset, String model_filename, String predictedvalue_filename); + + private: + void SVM_C_Train_(int learner_typeid, Dataset_sl& dataset, datanode *module); + void SVM_R_Train_(int learner_typeid, Dataset_sl& dataset, datanode *module); + void SVM_Q_Train_(int learner_typeid, Dataset_sl& dataset, datanode *module); + double SVM_C_Predict_(NZ_entry *test_vec); + double SVM_R_Predict_(NZ_entry *test_vec); + double SVM_Q_Predict_(NZ_entry *test_vec); + + void SaveModel_(int learner_typeid, String model_filename); + void LoadModel_(int learner_typeid, String model_filename); + char *ReadLine(FILE *fp); +}; + +template +char * SVM::ReadLine(FILE *fp_in) { + index_t length; + if ( fgets(line_, max_line_length_, fp_in)==NULL ) { + return NULL; + } + while ( strrchr(line_,'\n')==NULL ) { + max_line_length_ *= 2; + line_ = (char *) realloc(line_, max_line_length_); + length = (index_t) strlen(line_); + if ( fgets(line_+length, max_line_length_-length, fp_in)==NULL ) { + break; + } + } + return line_; +} + +template +void SVM::GetLabels(Dataset_sl &dataset, + ArrayList &labels_list, + ArrayList &labels_index, + ArrayList &labels_ct, + ArrayList &labels_startpos) { + index_t i = 0; + //index_t label_row_idx = matrix_.n_rows() - 1; // the last row is for labels + index_t n_points = dataset.n_points; + index_t n_labels = 0; + + double current_label; + + // these Arraylists need initialization before-hand + labels_list.Renew(); + labels_index.Renew(); + labels_ct.Renew(); + labels_startpos.Renew(); + + labels_index.Init(n_points); + labels_list.Init(); + labels_ct.Init(); + labels_startpos.Init(); + + ArrayList labels_temp; + labels_temp.Init(n_points); + labels_temp[0] = 0; + + labels_list.PushBack() = dataset.y[0]; + labels_ct.PushBack() = 1; + n_labels++; + + for (i = 1; i < n_points; i++) { + current_label = dataset.y[i]; + index_t j = 0; + for (j = 0; j < n_labels; j++) { + if (current_label == labels_list[j]) { + labels_ct[j]++; + break; + } + } + labels_temp[i] = j; + if (j == n_labels) { // new label + labels_list.PushBack() = current_label; // add new label to list + labels_ct.PushBack() = 1; + n_labels++; + } + } + + labels_startpos.PushBack() = 0; + for(i = 1; i < n_labels; i++){ + labels_startpos.PushBack() = labels_startpos[i-1] + labels_ct[i-1]; + } + + for(i = 0; i < n_points; i++) { + labels_index[labels_startpos[labels_temp[i]]] = i; + labels_startpos[labels_temp[i]]++; + } + + labels_startpos[0] = 0; + for(i = 1; i < n_labels; i++) { + labels_startpos[i] = labels_startpos[i-1] + labels_ct[i-1]; + } + + dataset.n_classes = n_labels; + + labels_temp.Clear(); +} + + +/** +* SVM initialization +* +* @param: labeled training set or testing set +* @param: number of classes (different labels) in the data set +* @param: module name +*/ +template +void SVM::Init(int learner_typeid, Dataset_sl& dataset, datanode *module){ + learner_typeid_ = learner_typeid; + + opt_method_ = fx_param_str(NULL, "opt", "smo"); // optimization method: default using SMO + + n_data_ = dataset.n_points; + + train_labels_list_.Init(); + train_labels_index_.Init(); + train_labels_ct_.Init(); + train_labels_startpos_.Init(); + + /* 1.Find the # of classes of the training set; + 2.Group labels, split the training dataset for training bi-class SVM classifiers */ + GetLabels(dataset, train_labels_list_, train_labels_index_, train_labels_ct_, train_labels_startpos_); + num_classes_ = dataset.n_classes; + + if (learner_typeid == 0) { /* for multiclass SVM classificatioin*/ + num_models_ = num_classes_ * (num_classes_-1) / 2; + sv_list_startpos_.Init(num_classes_); + sv_list_ct_.Init(num_classes_); + } + else { /* for other SVM learners */ + num_classes_ = 2; // dummy #, only meaningful in SaveModel and LoadModel + + num_models_ = 1; + sv_list_startpos_.Init(); + sv_list_ct_.Init(); + } + + models_.Init(); + sv_index_.Init(); + total_num_sv_ = 0; + + /* bool indicators FOR THE TRAINING SET: is/isn't a support vector */ + /* Note: it has the same index as the training !!! */ + trainset_sv_indicator_.Init(n_data_); + for (index_t i=0; i +void SVM::InitTrain(int learner_typeid, Dataset_sl& dataset, datanode *module) { + Init(learner_typeid, dataset, module); + if (learner_typeid == 0) { // Multiclass SVM Clssification + SVM_C_Train_(learner_typeid, dataset, module); + } + else if (learner_typeid == 1) { // SVM Regression + SVM_R_Train_(learner_typeid, dataset, module); + } + else if (learner_typeid == 2) { // SVM Quantile Estimation + SVM_Q_Train_(learner_typeid, dataset, module); + } + + /* Save models to file "svm_model" */ + SaveModel_(learner_typeid, "svm_model"); // TODO: param_req, and for CV mode + // TODO: calculate training error +} + + +/** +* Training for Multiclass SVM Clssification, using One-vs-One method +* +* @param: type id of the learner +* @param: training set +* @param: number of classes of the training set +* @param: module name +*/ +template +void SVM::SVM_C_Train_(int learner_typeid, Dataset_sl &dataset, datanode *module) { + /* Train num_classes*(num_classes-1)/2 binary class(labels:-1, 1) models */ + index_t ct = 0; + index_t i, j; + for (i = 0; i < num_classes_; i++) { + for (j = i+1; j < num_classes_; j++) { + models_.PushBack(); + /* Construct dataset consists of two classes i and j (reassign labels 1 and -1) */ + Dataset_sl dataset_bi; + dataset_bi.n_points = train_labels_ct_[i]+train_labels_ct_[j]; + dataset_bi.n_features = dataset.n_features; + dataset_bi.x = Malloc(struct NZ_entry *, dataset_bi.n_points); + dataset_bi.y = Malloc(double, dataset_bi.n_points); + ArrayList dataset_bi_index; + dataset_bi_index.Init(dataset_bi.n_points); + // for class 1 + for (index_t m = 0; m < train_labels_ct_[i]; m++) { + dataset_bi.x[m] = dataset.x[train_labels_index_[train_labels_startpos_[i]+m]]; + dataset_bi.y[m] = 1; + dataset_bi_index[m] = train_labels_index_[train_labels_startpos_[i]+m]; + } + // flor class -1 + for (index_t n = 0; n < train_labels_ct_[j]; n++) { + dataset_bi.x[n+train_labels_ct_[i]] = dataset.x[train_labels_index_[train_labels_startpos_[j]+n]]; + dataset_bi.y[n+train_labels_ct_[i]] = -1; + dataset_bi_index[n+train_labels_ct_[i]] = train_labels_index_[train_labels_startpos_[j]+n]; + } + + if (opt_method_== "smo") { + // Initialize SMO parameters + ArrayList param_feed_db; + param_feed_db.Init(); + param_feed_db.PushBack() = param_.b_; + param_feed_db.PushBack() = param_.Cp_; + param_feed_db.PushBack() = param_.Cn_; + param_feed_db.PushBack() = param_.hinge_; + param_feed_db.PushBack() = param_.wss_; + param_feed_db.PushBack() = param_.n_iter_; + param_feed_db.PushBack() = param_.accuracy_; + SMO smo; + smo.InitPara(learner_typeid, param_feed_db); + + // Initialize kernel + smo.kernel().Init(fx_submodule(module, "kernel")); + + // 2-classes SVM training using SMO + fx_timer_start(NULL, "train_smo"); + smo.Train(learner_typeid, dataset_bi); + fx_timer_stop(NULL, "train_smo"); + + // Get the trained bi-class model + models_[ct].coef_.Init(); // alpha*y + models_[ct].bias_ = smo.Bias(); // bias + //models_[ct].w_.Init(0); // for linear classifiers only. not used here + smo.GetSV(dataset_bi_index, models_[ct].coef_, trainset_sv_indicator_); // get support vectors + } + else if (opt_method_== "sgd") { + // Initialize SGD parameters + ArrayList param_feed_db; + param_feed_db.Init(); + param_feed_db.PushBack() = param_.Cp_; + param_feed_db.PushBack() = param_.Cn_; + param_feed_db.PushBack() = param_.kerneltypeid_== 0 ? 0.0: 1.0; + param_feed_db.PushBack() = param_.n_epochs_; + param_feed_db.PushBack() = param_.n_iter_; + param_feed_db.PushBack() = param_.accuracy_; + SGD sgd; + sgd.InitPara(learner_typeid, param_feed_db); + + // Initialize kernel + sgd.kernel().Init(fx_submodule(module, "kernel")); + + // 2-classes SVM training using SGD + fx_timer_start(NULL, "train_sgd"); + sgd.Train(learner_typeid, dataset_bi); + fx_timer_stop(NULL, "train_sgd"); + // Get the trained bi-class model + models_[ct].coef_.Init(); // alpha*y, used for nonlinear SVM only + if (param_.kerneltypeid_== 0) { // linear SVM + sgd.GetW(models_[ct].w_); + models_[ct].scale_w_ = sgd.ScaleW(); // scale of w for linear SVM. Use it if w's scaling is not done in training session + } + else { // nonlinear SVM + sgd.GetSV(dataset_bi_index, models_[ct].coef_, trainset_sv_indicator_); // get support vectors + //models_[ct].w_.Init(0); // for linear SVM only. not used here + } + models_[ct].bias_ = sgd.Bias(); // bias + } + else if (opt_method_== "md") { + // Initialize MD parameters + ArrayList param_feed_db; + param_feed_db.Init(); + param_feed_db.PushBack() = param_.Cp_; + param_feed_db.PushBack() = param_.Cn_; + param_feed_db.PushBack() = param_.n_epochs_; + param_feed_db.PushBack() = param_.n_iter_; + param_feed_db.PushBack() = param_.accuracy_; + MD md; + md.InitPara(learner_typeid, param_feed_db); + + // Initialize kernel + md.kernel().Init(fx_submodule(module, "kernel")); + + // 2-classes SVM training using MD + fx_timer_start(NULL, "train_md"); + md.Train(learner_typeid, dataset_bi); + fx_timer_stop(NULL, "train_md"); + // Get the trained bi-class model + models_[ct].coef_.Init(); // alpha*y, not used here + md.GetW(models_[ct].w_); + models_[ct].scale_w_ = md.ScaleW(); // scale of w for linear SVM. Use it if w's scaling is not done in training session + } + else if (opt_method_== "tgd") { + // Initialize TGD parameters + ArrayList param_feed_db; + param_feed_db.Init(); + param_feed_db.PushBack() = param_.Cp_; + param_feed_db.PushBack() = param_.Cn_; + param_feed_db.PushBack() = param_.n_epochs_; + param_feed_db.PushBack() = param_.n_iter_; + param_feed_db.PushBack() = param_.accuracy_; + TGD tgd; + tgd.InitPara(learner_typeid, param_feed_db); + + // Initialize kernel + tgd.kernel().Init(fx_submodule(module, "kernel")); + + // 2-classes SVM training using TGD + fx_timer_start(NULL, "train_tgd"); + tgd.Train(learner_typeid, dataset_bi); + fx_timer_stop(NULL, "train_tgd"); + // Get the trained bi-class model + models_[ct].coef_.Init(); // alpha*y, not used here + tgd.GetW(models_[ct].w_); + models_[ct].scale_w_ = tgd.ScaleW(); // scale of w for linear SVM. Use it if w's scaling is not done in training session + } + else { + fprintf(stderr, "ERROR!!! Unknown optimization method!\n"); + } + + ct++; + } + } + + /* Get total set of SVs from all the binary models */ + index_t k; + sv_list_startpos_[0] = 0; + + for (i = 0; i < num_classes_; i++) { + ct = 0; + for (j = 0; j < train_labels_ct_[i]; j++) { + if (trainset_sv_indicator_[ train_labels_index_[train_labels_startpos_[i]+j] ]) { + sv_index_.PushBack() = train_labels_index_[train_labels_startpos_[i]+j]; + total_num_sv_++; + ct++; + } + } + sv_list_ct_[i] = ct; + if (i >= 1) + sv_list_startpos_[i] = sv_list_startpos_[i-1] + sv_list_ct_[i-1]; + } + //sv_.Init(num_features_, total_num_sv_); + sv_entries_ = Malloc(struct NZ_entry *, total_num_sv_); + for (i = 0; i < total_num_sv_; i++) { + sv_entries_[i] = (dataset.x)[sv_index_[i]]; + } + + /* Get the matrix sv_coef_ which stores the coefficients of all sets of SVs */ + /* i.e. models_[x].coef_ -> sv_coef_ */ + index_t ct_model = 0; + index_t p; + sv_coef_.Init(num_classes_-1, total_num_sv_); + sv_coef_.SetZero(); + for (i = 0; i < num_classes_; i++) { + for (j = i+1; j < num_classes_; j++) { + p = sv_list_startpos_[i]; + for (k = 0; k < train_labels_ct_[i]; k++) { + if (trainset_sv_indicator_[ train_labels_index_[train_labels_startpos_[i]+k] ]) { + sv_coef_.set(j-1, p++, models_[ct_model].coef_[k]); + } + } + p = sv_list_startpos_[j]; + for (k = 0; k < train_labels_ct_[j]; k++) { + if (trainset_sv_indicator_[ train_labels_index_[train_labels_startpos_[j]+k] ]) { + sv_coef_.set(i, p++, models_[ct_model].coef_[train_labels_ct_[i] + k]); + } + } + ct_model++; + } + } +} + +/** +* Training for SVM Regression +* +* @param: type id of the learner +* @param: training set +* @param: module name +*/ +template +void SVM::SVM_R_Train_(int learner_typeid, Dataset_sl &dataset, datanode *module) { + index_t i; + ArrayList dataset_index; + dataset_index.Init(n_data_); + for (i=0; i param_feed_db; + param_feed_db.Init(); + param_feed_db.PushBack() = param_.b_; + param_feed_db.PushBack() = param_.C_; + param_feed_db.PushBack() = param_.epsilon_; + param_feed_db.PushBack() = param_.wss_; + param_feed_db.PushBack() = param_.n_iter_; + param_feed_db.PushBack() = param_.accuracy_; + SMO smo; + smo.InitPara(learner_typeid, param_feed_db); + + // Initialize kernel + smo.kernel().Init(fx_submodule(module, "kernel")); + + // SVM_R Training using SMO + smo.Train(learner_typeid, dataset); + + // Get the trained model + models_[0].bias_ = smo.Bias(); // bias + models_[0].coef_.Init(); // alpha*y + //models_[0].w_.Init(0); // not using + smo.GetSV(dataset_index, models_[0].coef_, trainset_sv_indicator_); // get support vectors + } + else if (opt_method_== "sgd") { + // Initialize SGD parameters + ArrayList param_feed_db; + param_feed_db.Init(); + param_feed_db.PushBack() = param_.Cp_; + param_feed_db.PushBack() = param_.Cn_; + param_feed_db.PushBack() = param_.kerneltypeid_== 0 ? 0.0: 1.0; + SGD sgd; + sgd.InitPara(learner_typeid, param_feed_db); + + // Initialize kernel + sgd.kernel().Init(fx_submodule(module, "kernel")); + + // SVM_R Training using SGD + sgd.Train(learner_typeid, dataset); + + // Get the trained model + models_[0].bias_ = sgd.Bias(); // bias + //models_[0].w_.Copy(*(sgd.GetW())); // w + sgd.GetW(models_[0].w_); + models_[0].scale_w_ = sgd.ScaleW(); // scale of w for linear SVM. Use it if w's scaling is not done in training session + models_[0].coef_.Init(0); // not using + } + else { + fprintf(stderr, "ERROR!!! Unknown optimization method!"); + } + + /* Get index list of support vectors */ + for (i = 0; i < n_data_; i++) { + if (trainset_sv_indicator_[i]) { + sv_index_.PushBack() = i; + total_num_sv_++; + } + } + + /* Get support vecotors and coefficients */ + //sv_.Init(num_features_, total_num_sv_); + sv_entries_ = Malloc(struct NZ_entry*, total_num_sv_); + for (i = 0; i < total_num_sv_; i++) { + sv_entries_[i] = dataset.x[sv_index_[i]]; + } + + /* + for (i = 0; i < total_num_sv_; i++) { + Vector source, dest; + sv_.MakeColumnVector(i, &dest); + // last row of dataset is for labels + dataset.matrix().MakeColumnSubvector(sv_index_[i], 0, num_features_, &source); + dest.CopyValues(source); + } + */ + + sv_coef_.Init(1, total_num_sv_); + for (i = 0; i < total_num_sv_; i++) { + sv_coef_.set(0, i, models_[0].coef_[i]); + } + +} + +/** +* Training for SVM Quantile Estimation +* +* @param: type id of the learner +* @param: training set +* @param: module name +*/ +template +void SVM::SVM_Q_Train_(int learner_typeid, Dataset_sl &dataset, datanode *module) { + // TODO +} + + +/** +* SVM prediction for one testing vector +* +* @param: type id of the learner +* @param: testing vector +* +* @return: predited value +*/ +template +double SVM::Predict(int learner_typeid, struct NZ_entry *test_vec) { + double predicted_value = INFINITY; + if (learner_typeid == 0) { // Multiclass SVM Clssification + predicted_value = SVM_C_Predict_(test_vec); + } + else if (learner_typeid == 1) { // SVM Regression + predicted_value = SVM_R_Predict_(test_vec); + } + else if (learner_typeid == 2) { // SVM Quantile Estimation + predicted_value = SVM_Q_Predict_(test_vec); + } + return predicted_value; +} + +/** +* Multiclass SVM classification for one testing vector +* +* @param: testing vector +* +* @return: a label (double-type-integer, e.g. 1.0, 2.0, 3.0) +*/ +template +double SVM::SVM_C_Predict_(struct NZ_entry *test_vec) { + index_t i, j, k; + ArrayList keval; + keval.Init(total_num_sv_); + if (opt_method_!="sgd" || param_.kerneltypeid_ != 0) { + for (i = 0; i < total_num_sv_; i++) { + keval[i] = param_.kernel_.Eval(test_vec, sv_entries_[i]); + } + } + ArrayList values; + values.Init(num_models_); + index_t ct = 0; + double sum = 0.0; + for (i = 0; i < num_classes_; i++) { + for (j = i+1; j < num_classes_; j++) { + if (opt_method_== "smo") { + sum = 0.0; + for (k = 0; k < sv_list_ct_[i]; k++) { + sum += sv_coef_.get(j-1, sv_list_startpos_[i]+k) * keval[sv_list_startpos_[i]+k]; + } + for (k = 0; k < sv_list_ct_[j]; k++) { + sum += sv_coef_.get(i, sv_list_startpos_[j]+k) * keval[sv_list_startpos_[j]+k]; + } + sum += models_[ct].bias_; + } + else if (opt_method_== "sgd") { + if (param_.kerneltypeid_== 0) { // linear + sum = SparseDot(models_[ct].w_, test_vec); + sum *= models_[ct].scale_w_; // Use this if scaling of w is not done in the training session + } + else { //nonlinear + sum = 0.0; + for (k = 0; k < sv_list_ct_[i]; k++) { + sum += sv_coef_.get(j-1, sv_list_startpos_[i]+k) * keval[sv_list_startpos_[i]+k]; + } + for (k = 0; k < sv_list_ct_[j]; k++) { + sum += sv_coef_.get(i, sv_list_startpos_[j]+k) * keval[sv_list_startpos_[j]+k]; + } + } + sum += models_[ct].bias_; + } + else if (opt_method_== "md" || opt_method_== "tgd") { + sum = SparseDot(models_[ct].w_, test_vec); + //sum *= models_[ct].scale_w_; // Use this if scaling of w is not done in the training session + } + values[ct] = sum; + ct++; + } + } + + ArrayList vote; + vote.Init(num_classes_); + for (i = 0; i < num_classes_; i++) { + vote[i] = 0; + } + ct = 0; + for (i = 0; i < num_classes_; i++) { + for (j = i+1; j < num_classes_; j++) { + if(values[ct] > 0.0) { // label 1 in bi-classifiers (for i=...) + vote[i] = vote[i] + 1; + } + else { // label -1 in bi-classifiers (for j=...) + vote[j] = vote[j] + 1; + } + ct++; + } + } + index_t vote_max_idx = 0; + for (i = 1; i < num_classes_; i++) { + if (vote[i] >= vote[vote_max_idx]) { + vote_max_idx = i; + } + } + return train_labels_list_[vote_max_idx]; +} + +/** +* SVM Regression Prediction for one testing vector +* +* @param: testing vector +* +* @return: predicted regression value +*/ +template +double SVM::SVM_R_Predict_(struct NZ_entry *test_vec) { + index_t i; + double sum = 0.0; + if (opt_method_== "smo") { + for (i = 0; i < total_num_sv_; i++) { + //sum += sv_coef_.get(0, i) * param_.kernel_.Eval(datum.ptr(), sv_.GetColumnPtr(i), num_features_); + sum += sv_coef_.get(0, i) * param_.kernel_.Eval(test_vec, sv_entries_[i]); + } + } + else if (opt_method_== "sgd") { + // TODO + } + sum += models_[0].bias_; + return sum; +} + +/** +* SVM Quantile Estimation Prediction for one testing vector +* +* @param: testing vector +* +* @return: estimated quantile value (the support) +*/ +template +double SVM::SVM_Q_Predict_(struct NZ_entry *test_vec) { + // TODO + return 0.0; +} + + + +/** +* Online batch classification for multiple testing vectors. No need to load model file, +* since models are already in RAM. +* +* Note: for test set, if no true test labels provided, just put some dummy labels +* (e.g. all -1) in the last row of testset +* +* @param: type id of the learner +* @param: testing set +* @param: file name of the testing data +*/ +template +void SVM::BatchPredict(int learner_typeid, Dataset_sl& testset, String predictedvalue_filename) { + FILE *fp = fopen(predictedvalue_filename, "w"); + double predictedvalue; + if (fp == NULL) { + fprintf(stderr, "Cannot save predicted values to file!\n"); + return; + } + index_t err_ct = 0; + for (index_t i = 0; i < testset.n_points; i++) { + predictedvalue = Predict(learner_typeid, (testset.x)[i]); + if (predictedvalue != testset.y[i]) { + err_ct++; + } + // save predicted values to file + fprintf(fp, "%f\n", predictedvalue); + } + fclose(fp); + /* calculate testing error */ + printf( "\n*** %d out of %d misclassified ***\n", err_ct, testset.n_points ); + printf( "*** Testing error is %f, accuracy is %f. ***\n", double(err_ct)/double(testset.n_points), 1- double(err_ct)/double(testset.n_points) ); + //fprintf( stderr, "*** Results are save in \"%s\" ***\n\n", predictedvalue_filename.c_str()); +} + +/** +* Load models from a file, and perform offline batch classification for multiple testing vectors +* +* @param: type id of the learner +* @param: testing set +* @param: name of the model file +* @param: name of the file to store classified labels +*/ +template +void SVM::LoadModelBatchPredict(int learner_typeid, Dataset_sl& testset, String model_filename, String predictedvalue_filename) { + LoadModel_(learner_typeid, model_filename); + BatchPredict(learner_typeid, testset, predictedvalue_filename); +} + + +/** +* Save SVM model to a text file +* +* @param: type id of the learner +* @param: name of the model file +*/ +// TODO: use XML +template +void SVM::SaveModel_(int learner_typeid, String model_filename) { + FILE *fp = fopen(model_filename, "w"); + if (fp == NULL) { + fprintf(stderr, "Cannot save trained model to file!"); + return; + } + index_t i, j; + + if (learner_typeid == 0) { // for SVM_C + fprintf(fp, "svm_type SVM_C\n"); + fprintf(fp, "total_num_sv %d\n", total_num_sv_); + fprintf(fp, "num_classes %d\n", num_classes_); + // save labels + fprintf(fp, "labels "); + for (i = 0; i < num_classes_; i++) + fprintf(fp, "%f ", train_labels_list_[i]); + fprintf(fp, "\n"); + // save support vector info + fprintf(fp, "sv_list_startpos "); + for (i =0; i < num_classes_; i++) + fprintf(fp, "%d ", sv_list_startpos_[i]); + fprintf(fp, "\n"); + fprintf(fp, "sv_list_ct "); + for (i =0; i < num_classes_; i++) + fprintf(fp, "%d ", sv_list_ct_[i]); + fprintf(fp, "\n"); + } + else if (learner_typeid == 1) { // for SVM_R + fprintf(fp, "svm_type SVM_R\n"); + fprintf(fp, "total_num_sv %d\n", total_num_sv_); + fprintf(fp, "sv_index "); + for (i = 0; i < total_num_sv_; i++) + fprintf(fp, "%d ", sv_index_[i]); + fprintf(fp, "\n"); + } + else if (learner_typeid == 2) { // for SVM_Q + fprintf(fp, "svm_type SVM_Q\n"); + fprintf(fp, "total_num_sv %d\n", total_num_sv_); + fprintf(fp, "sv_index "); + for (i = 0; i < total_num_sv_; i++) + fprintf(fp, "%d ", sv_index_[i]); + fprintf(fp, "\n"); + } + + // save kernel parameters + fprintf(fp, "kernel_name %s\n", param_.kernelname_.c_str()); + fprintf(fp, "kernel_typeid %d\n", param_.kerneltypeid_); + param_.kernel_.SaveParam(fp); + + // save models: bias, coefficients and support vectors + fprintf(fp, "bias "); + for (i = 0; i < num_models_; i++) + fprintf(fp, "%.16g ", models_[i].bias_); + fprintf(fp, "\n"); + + fprintf(fp, "SV_coefs\n"); + for (i = 0; i < total_num_sv_; i++) { + for (j = 0; j < num_classes_-1; j++) { + fprintf(fp, "%.16g ", sv_coef_.get(j,i)); + } + const struct NZ_entry *sv_p = sv_entries_[i]; + while (sv_p->index != -1) { + fprintf(fp, "%d:%.8g ", sv_p->index+1, sv_p->value); // in svmlight's data format, feature index begins from 1, not 0 + sv_p ++; + } + fprintf(fp, "\n"); + } + + fclose(fp); +} + +/** +* Load SVM model file +* +* @param: type id of the learner +* @param: name of the model file +*/ +// TODO: use XML +template +void SVM::LoadModel_(int learner_typeid, String model_filename) { + if (learner_typeid == 0) {// SVM_C + train_labels_list_.Renew(); + train_labels_list_.Init(num_classes_); // get labels list from the model file + } + + /* load model file */ + FILE *fp = fopen(model_filename, "r"); + if (fp == NULL) { + fprintf(stderr, "Cannot open SVM model file!"); + return; + } + char cmd[80]; + int i, j; int temp_d; double temp_f; + for (i = 0; i < num_models_; i++) { + models_.PushBack(); + models_[i].coef_.Init(); + } + while (1) { + fscanf(fp,"%80s",cmd); + if(strcmp(cmd,"svm_type")==0) { + fscanf(fp,"%80s", cmd); + if (strcmp(cmd,"SVM_C")==0) + learner_typeid_ = 0; + else if (strcmp(cmd,"SVM_R")==0) + learner_typeid_ = 1; + else if (strcmp(cmd,"SVM_Q")==0) + learner_typeid_ = 2; + } + else if (strcmp(cmd, "total_num_sv")==0) { + fscanf(fp,"%d",&total_num_sv_); + } + // for SVM_C + else if (strcmp(cmd, "num_classes")==0) { + fscanf(fp,"%d",&num_classes_); + } + else if (strcmp(cmd, "labels")==0) { + for (i=0; i *x; // TODO +}; + + +/** +* Class for Linear Kernel +*/ +class SVMLinearKernel { + public: + // Init of kernel parameters + ArrayList kpara_; // kernel parameters + void Init(datanode *node) { //TODO: NULL->node + kpara_.Init(); + } + // Kernel name + void GetName(String* kname) { + kname->Copy("linear"); + } + // Get an type ID for kernel + int GetTypeId() { + return ID_LINEAR; + } + // Kernel value evaluation + double Eval(NZ_entry *pa, NZ_entry *pb) { + // sparse dot product + double kv = 0.0; + if (pa->index == -1 || pb->index == -1) { + return 0.0; + } + else { + while (pa->index !=-1 && pb->index !=-1) { + if (pa->index == pb->index) { + kv += pa->value * pb->value; + ++pa; + ++pb; + } + else { + if (pa->index > pb->index) { + ++pb; + } + else { + ++pa; + } + } + } + } + return kv; + } + // Save kernel parameters to file + void SaveParam(FILE* fp) { + } +}; + +/** +* Class for Gaussian RBF Kernel +*/ +class SVMRBFKernel { + public: + // Init of kernel parameters + ArrayList kpara_; // kernel parameters + void Init(datanode *node) { //TODO: NULL->node + kpara_.Init(2); + kpara_[0] = fx_param_double_req(NULL, "sigma"); // sigma + kpara_[1] = -1.0 / (2 * kpara_[0] * kpara_[0]); // -gamma = -1/(2 sigma^2) + } + // Kernel name + void GetName(String* kname) { + kname->Copy("gaussian"); + } + // Get an type ID for kernel + int GetTypeId() { + return ID_GAUSSIAN; + } + // Kernel value evaluation + double Eval(NZ_entry *pa, NZ_entry *pb) { + double kv = 0; + while (pa->index !=-1 && pb->index !=-1) { + if (pa->index == pb->index) { + double tmp = pa->value - pb->value; + kv += tmp * tmp; + ++pa; + ++pb; + } + else { + if (pa->index > pb->index) { + kv += pb->value * pb->value; + ++pb; + } + else { + kv += pa->value * pa->value; + ++pa; + } + } + } + while (pa->index != -1) { + kv += pa->value * pa->value; + ++pa; + } + while (pb->index != -1) { + kv += pb->value * pb->value; + ++ pb; + } + kv = exp( kpara_[1] * kv ); + return kv; + } + // Save kernel parameters to file + void SaveParam(FILE* fp) { + fprintf(fp, "sigma %g\n", kpara_[0]); + fprintf(fp, "gamma %g\n", kpara_[1]); + } +}; + + +/** + * Sparse vector scaling: w<= scale *w + */ +void SparseScale(double scale, ArrayList &w) { + if (scale == 0) { + w.ShrinkTo(0); + } + else { + for (index_t i=0; i &w, NZ_entry *x) { + double dv = 0.0; + index_t w_nz = w.size(); + index_t ct = 0; + if (w_nz == 0 || x->index == -1) { + return 0.0; + } + else { + while (ctindex !=-1) { + if (w[ct].index == x->index) { + dv += w[ct].value * x->value; + ++ct; + ++x; + } + else { + if (w[ct].index > x->index) { + ++x; + } + else { + ++ct; + } + } + } + return dv; + } +} + +/** + * Sparse vector scaled add: w<= w+ scale * x + */ +void SparseAddExpert(double scale, NZ_entry *x, ArrayList &w) { + index_t ct = 0; + NZ_entry nz_tmp; + + if (x->index == -1 || scale == 0) { // x: all zeros or scale ==0 + // w remains unchanged + } + else if (w.size() == 0) { // w: all zeros + while (x->index != -1) { + nz_tmp.index = x->index; + nz_tmp.value = scale * x->value; + w.PushBackCopy(nz_tmp); + ++x; + } + } + else { // neither w nor x is of all zeros + while (ctindex !=-1) { + if (ct == w.size()) { // w reaches end, while x still not + nz_tmp.index = x->index; + nz_tmp.value = scale * x->value; + w.PushBackCopy(nz_tmp); + ++ct; + ++x; + } + else if (x->index == -1) { // x reaches end, while w still not + // the succeeding w remain unchanged + break; + } + else { // neither w nor x reaches end + if (w[ct].index == x->index) { + w[ct].value = w[ct].value + scale * x->value; + ++ct; + ++x; + } + else if (w[ct].index > x->index) { + nz_tmp.index = x->index; + nz_tmp.value = scale * x->value; + w.InsertCopy(ct, nz_tmp); // w's dimension increases by 1 + ++ct; + ++x; + } + else { // w[ct].index < x->index + // w[ct] remains unchanged + ++ct; + } + } + } + // shrink w + for (ct=0; ct &x,ArrayList &y, ArrayList &z) { + index_t x_size, y_size; + index_t i, ct_x, ct_y; + x_size = x.size(); + y_size = y.size(); + ct_x = 0; ct_y = 0; + NZ_entry nz_tmp; + + /* + for (index_t i=0; i x[ct_x].index) { + nz_tmp.index = x[ct_x].index; + nz_tmp.value = - x[ct_x].value; + z.InsertCopy(ct_y, nz_tmp); // w's dimension increases by 1 + ++ct_x; + } + else { // y[ct].index < x[ct_x].index + ++ct_y; + } + } + } + + } + + /* + for (index_t i=0; i &x,ArrayList &y) { + index_t x_size; + index_t i, ct_x, ct_y; + x_size = x.size(); + ct_x = 0; ct_y = 0; + NZ_entry nz_tmp; + + if (y.size() == 0) { // y: all zeros + y.GrowTo(x_size); + for (i=0; i x[ct_x].index) { + nz_tmp.index = x[ct_x].index; + nz_tmp.value = - x[ct_x].value; + y.InsertCopy(ct_y, nz_tmp); // w's dimension increases by 1 + ++ct_x; + ++ct_y; + } + else { // y[ct].index < x[ct_x].index + ++ct_y; + } + } + } + + } + +} + + + + + +/** + * Sparse element-wise multiplication of vectors: y <= y .* x + */ +void SparseElementMulOverwrite(ArrayList &y, ArrayList &x) { + index_t ct_x = 0; + index_t ct_y = 0; + + if (y.size() == 0) { // y: all zeros + // y remains unchanged + } + else if (x.size() == 0) { // x: all zeros + y.ShrinkTo(0); + } + else { // neither x nor y is of all zeros + while (ct_x < x.size() || ct_y < y.size()) { + if (ct_x == x.size()) { // x reaches end, while y still not + y.Remove(ct_y); + //printf("ct_y=%d\n", ct_y); + } + else if (ct_y == y.size()) { // y reaches end, while x still not + break; + } + else { // neither x nor y reaches end + if (x[ct_x].index == y[ct_y].index) { + y[ct_y].value = y[ct_y].value * x[ct_x].value; + ++ct_x; + ++ct_y; + } + else if (y[ct_y].index > x[ct_x].index) { + ++ct_x; + //printf("y[%d].index=%d; x[%d].index=%d\n",ct_y, y[ct_y].index, ct_x, x[ct_x].index); + } + else { // y[ct].index < x[ct_x].index + y.Remove(ct_y); + } + } + } + } + +} + + +#endif