apaam/mlpack
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Initial release of "regularized risk minimization" project

Browse Source
This commit is contained in:
houyang
2010-03-17 20:10:13 +00:00
parent 7b44b8b51b
commit c7a38c968a
8 changed files with 4207 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files
fastlib/trunk/contrib/houyang/regmin/build.py
+23
View File
@@ -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"]
)
fastlib/trunk/contrib/houyang/regmin/opt_md.h
+420
View File
@@ -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<typename TKernel>
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<int> y_; /* list that stores "labels" */
ArrayList<NZ_entry> w_; /* the slope of the decision hyperplane, including bias: [w, b] */
ArrayList<NZ_entry> w_p_; /* coefficients for positive w_t+ */
ArrayList<NZ_entry> 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<index_t> old_from_new_; // for generating a random sequence of training data
public:
MD() {}
~MD() {}
/**
* Initialization for parameters
*/
void InitPara(int learner_typeid, ArrayList<double> &param_) {
// 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<NZ_entry> &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<wp_size; i++) {
if (fabs(w_p_[i].value) >= 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<index_t> &dataset_index, ArrayList<double> &coef, ArrayList<bool> &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<typename TKernel>
void MD<TKernel>::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; i<n_features_; i++) {
w_p_[i].index = i;
w_p_[i].value = C_iv_2f; // TODO
w_n_[i].index = i;
w_n_[i].value = C_iv_2f; // TODO
}
y_.Init(n_data_);
for (i = 0; i < n_data_; i++) {
y_[i] = (dataset_->y)[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<typename TKernel>
void MD<TKernel>::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; epo<n_epochs_; epo++) {
/* To mimic the online learning senario, in each epoch,
we randomly permutate the training set, indexed by old_from_new_ */
for (i=0; i<n_data_; i++) {
old_from_new_[i] = i;
}
for (i=0; i<n_data_; i++) {
j = rand() % n_data_;
swap(old_from_new_[i], old_from_new_[j]);
}
ct = 0;
while (ct <= n_iter_) {
work_idx_old = old_from_new_[ct % n_data_];
//eta_ = 1.0 /(lambda_ * sqrt(t_));
//eta_ = 1.0 / (lambda_ * t_);
//eta_ = 1.0 / t_;
eta_ = 1.0 /sqrt(t_);
NZ_entry *xt;
xt = (dataset_->x)[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<w_p_.size(); i++) {
w_sum += w_p_[i].value;
}
for (i=0; i<w_n_.size(); i++) {
w_sum += w_n_[i].value;
}
if (w_sum > 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<w_p_.size(); i++)
printf("w_p[%d]=%f\n", i, w_p_[i].value);
for (i=0; i<w_n_.size(); i++)
printf("w_n[%d]=%f\n", i, w_n_[i].value);
for (i=0; i<w_.size(); i++)
printf("w[%d]=%f\n", i, w_[i].value);
*/
t_ += 1.0;
ct ++;
}
}// for epo
SparseSubOverwrite(w_n_, w_p_); // w_p <= w_p - w_n
// rounding w
index_t wp_size = w_p_.size();
round_thd_ = fx_param_double(NULL, "thd", 1.0e-5);
w_nnz_ = 0;
for (ct=0; ct<wp_size; ct++) {
if (fabs(w_p_[ct].value) >= 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<w_.size(); i++) {
wi_abs = fabs(w_[i].value);
if (wi_abs > 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<w_.size(); i++) {
//printf("w[%d]=%f\n", i, w_[i].value);
if ( fabs(w_[i].value) > 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
fastlib/trunk/contrib/houyang/regmin/opt_sgd.h
+564
View File
@@ -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<typename TKernel>
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<int> y_; /* list that stores "labels" */
ArrayList<NZ_entry> 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<index_t> old_from_new_; // for generating a random sequence of training data
//ArrayList<index_t> 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<double> &param_) {
// 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<NZ_entry> &w_out) {
w_out.InitCopy(w_, w_.size());
}
double ScaleW() const {
return scale_w_;
}
void GetSV(ArrayList<index_t> &dataset_index, ArrayList<double> &coef, ArrayList<bool> &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<typename TKernel>
void SGD<TKernel>::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<typename TKernel>
void SGD<TKernel>::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; epo<n_epochs_; epo++) {
/* To mimic the online learning senario, in each epoch,
we randomly permutate the training set, indexed by old_from_new_ */
for (i=0; i<n_data_; i++) {
old_from_new_[i] = i;
}
for (i=0; i<n_data_; i++) {
j = rand() % n_data_;
swap(old_from_new_[i], old_from_new_[j]);
}
ct = 0;
while (ct <= n_iter_) {
work_idx_old = old_from_new_[ct % n_data_];
eta_ = 1.0 / (lambda_ * t_); // update step size
scale_w_ = scale_w_ - scale_w_ / t_; // update scale of w
//la::Scale(scale_w, &w_); // Note: moving w's scaling calculation to the testing session is faster
if (scale_w_ < SGD_SCALE_W_TOLERANCE) {
SparseScale(scale_w_, w_);
scale_w_ = 1.0;
//printf("epo %d: scale_w tolerance reached.\n", epo);
}
NZ_entry *xt;
xt = (dataset_->x)[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<w_.size(); j++) {
v += w_[j].value * w_[j].value;
}
v = v * scale_w_ * scale_w_ / 2.0 + loss_sum;
printf("Primal objective value: %lf\n", v);
// find max(abs(w_i))
double wi_abs_max = -INFINITY;
double wi_abs;
for (i=0; i<w_.size(); i++) {
wi_abs = fabs(w_[i].value);
if (wi_abs > 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<w_.size(); i++) {
if ( fabs(w_[i].value) > 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<n_data_; i++) {
old_from_new_[i] = i;
}
for (i=0; i<n_data_; i++) {
j = rand() % n_data_;
swap(old_from_new_[i], old_from_new_[j]);
}
t_ = 1.0;
ct = 0;
printf("SGD_kernel training begins...\n");
while (ct < n_iter_) {
work_idx_old = old_from_new_[ct % n_data_];
yt = y_[work_idx_old];
yt_hat = 0.0;
n_real_epo = index_t(ceil(ct/n_data_));
n_data_res = ct - n_real_epo * n_data_;
if (n_real_epo > 0) {
for (i=0; i<ct; i++) {
if (fabs(coef_long[i]) >= SGD_ALPHA_ZERO) {
b_calc_kernel[i%n_data_] = 1;
}
}
for (i=0; i<n_data_res; i++) {
if (b_calc_kernel[i]>0) {
kernel_value = CalcKernelValue_(old_from_new_[i], work_idx_old);
for (j=0; j<n_real_epo+1; j++) {
yt_hat += coef_long[j*n_data_+i] * kernel_value;
}
}
}
for (i=n_data_res; i<n_data_; i++) {
if (b_calc_kernel[i]>0) {
kernel_value = CalcKernelValue_(old_from_new_[i], work_idx_old);
for (j=0; j<n_real_epo; j++) {
yt_hat += coef_long[j*n_data_+i] * kernel_value;
}
}
}
}
else { // ct < n_data_
for (i=0; i<ct; i++) {
if (fabs(coef_long[i]) >= 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<t)
one_minus_eta_lambda = 1.0 - eta_ * lambda_;
//printf("%d: %f\n", ct, one_minus_eta_lambda);
for (i=0; i<ct; i++) {
coef_long[i] = coef_long[i] * one_minus_eta_lambda;
}
// update current coef (for i==t)
//coef_long[ct] = eta_ * LossFunctionGradient_(learner_typeid, yy_hat) * yt; // Hard margin SVM
//printf("%f, %f, %f\n", eta_, LossFunctionGradient_(learner_typeid, yy_hat), yt);
// soft margin svm
if (yy_hat <= rho_) {
//printf("%d: %f, %f\n", ct, yy_hat, yt);
delta = 1.0;
}
else {
delta = 0.0;
}
coef_long[ct] = eta_ * delta * yt;
// update bias
//bias_ += coef_long[ct];
t_ += 1.0;
ct ++;
}
// convert coef_long to coef_
for (i=0; i<n_iter_; i++) {
work_idx_old = old_from_new_[i % n_data_];
coef_[work_idx_old] = coef_[work_idx_old] + coef_long[i];
}
} // else
}
/* Get results for nonlinear SGD: coefficients(alpha*y), 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<typename TKernel>
void SGD<TKernel>::GetSV(ArrayList<index_t> &dataset_index, ArrayList<double> &coef, ArrayList<bool> &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
fastlib/trunk/contrib/houyang/regmin/opt_smo.h
+929
View File
@@ -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 <class T> inline void swap(T& x, T& y) { T t=x; x=y; y=t; }
template<typename TKernel>
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<index_t> 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<int> 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<double> &param_) {
// 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<index_t> &dataset_index, ArrayList<double> &coef, ArrayList<bool> &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<typename TKernel>
void SMO<TKernel>::ReconstructGradient_() {
index_t i, j;
if (n_active_ == n_alpha_)
return;
if (learner_typeid_ == 0) { // SVM_C
for (i=n_active_; i<n_alpha_; i++) {
grad_[i] = 1 - grad_bar_[i];
}
}
else if (learner_typeid_ == 1) { // SVM_R
for (i=n_active_; i<n_alpha_; i++) {
j = 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<n_active_; i++) {
if (alpha_status_[i] == SMO_ID_FREE) {
for (j=n_active_; j<n_alpha_; j++) {
grad_[j] = grad_[j] - y_[j] * alpha_[i] * y_[i] * CalcKernelValue_(i,j);
}
}
}
}
/**
* Test whether need to do shrinking for provided index and y_grad_max, y_grad_min
*
*/
template<typename TKernel>
bool SMO<TKernel>::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_min
}
}
else if (IsLowerBounded(i)) {
if (y_[i] == 1) {
return (grad_[i] < y_grad_min);
}
else { // y_[i] == -1
return (grad_[i] + y_grad_max < 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<typename TKernel>
void SMO<TKernel>::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<n_active_; t++) { // find argmax(y*grad), t\in I_up
if (y_[t] == 1) {
if (!IsUpperBounded(t)) // t\in I_up, y==1: y[t]alpha[t] < C
if (grad_[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<n_active_; t++) { // find argmin(y*grad), t\in I_down
if (y_[t] == 1) {
if (!IsLowerBounded(t)) // t\in I_down, y==1: y[t]alpha[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<n_active_; t++) {
// Shrinking: put inactive alphas behind the active set
if (TestShrink_(t, y_grad_max, y_grad_min)) {
n_active_ --;
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_ --;
}
}
}
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<typename TKernel>
void SMO<TKernel>::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<typename TKernel>
void SMO<TKernel>::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; i<n_alpha_; i++) {
active_set_[i] = i;
}
alpha_status_.Init(n_alpha_);
for (i=0; i<n_alpha_; i++)
UpdateAlphaStatus_(i);
// initialize gradient (already set to init values)
/*
for (i=0; i<n_alpha_; i++) {
for(j=0; j<n_alpha_; j++) {
if (!IsLowerbounded(j)) { // alpha_j >0
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<n_alpha_; i++) {
for(j=0; j<n_alpha_; j++) {
if(IsUpperBounded(j)) // alpha_j >= 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<n_alpha_; i++)
// printf("%f.\n", y_[i]*alpha_[i]);
//printf("\n\n");
// for every min(n_data_, 1000) iterations, do shrinking
if (do_shrinking_) {
if ( --ct_shrinking_ == 0) {
Shrinking_();
ct_shrinking_ = min(n_data_, SMO_NUM_FOR_SHRINKING);
}
}
// Find working set, check stopping criterion, update gradient and alphas
stop_condition = SMOIterations_();
// Termination check, if stop_condition==1 or ==2 => 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<typename TKernel>
int SMO<TKernel>::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<typename TKernel>
bool SMO<TKernel>::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<n_active_; t++) { // find argmax(y*grad), t\in I_up
if (y_[t] == 1) {
if (!IsUpperBounded(t)) // t\in I_up, y==1: y[t]alpha[t] < C
if (grad_[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<n_active_; t++) { // find argmin(y*grad), t\in I_down
if (y_[t] == 1) {
if (!IsLowerBounded(t)) // t\in I_down, y==1: y[t]alpha[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<n_active_; t++) {
double K_it = CalcKernelValue_(out_i, t);
double K_tt = CalcKernelValue_(t, t);
if (y_[t] == 1) {
if (!IsLowerBounded(t)) { // t\in I_down, y==1: y[t]alpha[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<typename TKernel>
void SMO<TKernel>::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<n_active_; t++) {
grad_[t] = grad_[t] - y_[t] * (y_[i] * diff_i * CalcKernelValue_(i, t) + y_[j] * diff_j * CalcKernelValue_(j, t));
}
bool ub_i = IsUpperBounded(i);
bool ub_j = IsUpperBounded(j);
// Update alpha active status
UpdateAlphaStatus_(i);
UpdateAlphaStatus_(j);
if (do_shrinking_) {
// Update gradient_bar
if( ub_i != IsUpperBounded(i) ) { // updated_alpha_i >= C
if(ub_i) // old_alpha_i >= C, new_alpha_i < C
for(t=0; t<n_alpha_; t++)
grad_bar_[t] = grad_bar_[t] - C_i * y_[i] * y_[t] * CalcKernelValue_(i, t);
else // old_alpha_i < C, new_alpha_i >= C
for(t=0; t<n_alpha_; t++)
grad_bar_[t] = grad_bar_[t] + C_i * y_[i] * y_[t] * CalcKernelValue_(i, t);
}
if( ub_j != IsUpperBounded(j) ) {
if(ub_j) // old_alpha_j >= C, new_alpha_j < C
for(t=0; t<n_alpha_; t++)
grad_bar_[t] = grad_bar_[t] - C_j * y_[j] * y_[t] * CalcKernelValue_(j, t);
else // old_alpha_j < C, new_alpha_j >= C
for(t=0; t<n_alpha_; t++)
grad_bar_[t] = grad_bar_[t] + C_j * y_[j] * y_[t] * CalcKernelValue_(j, t);
}
}
}
/**
* Calcualte bias term
*
* @return: the bias
*
*/
template<typename TKernel>
void SMO<TKernel>::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<n_active_; i++){
double yg = y_[i] * grad_[i];
if (IsUpperBounded(i)) { // bounded: alpha_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 <C
n_free_alpha++;
sum_free_yg += yg;
}
}
if(n_free_alpha>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<typename TKernel>
void SMO<TKernel>::GetSV(ArrayList<index_t> &dataset_index, ArrayList<double> &coef, ArrayList<bool> &sv_indicator) {
ArrayList<index_t> 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
fastlib/trunk/contrib/houyang/regmin/opt_tgd.h
+390
View File
@@ -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<typename TKernel>
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<int> y_; /* list that stores "labels" */
ArrayList<NZ_entry> 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<index_t> old_from_new_; // for generating a random sequence of training data
public:
TGD() {}
~TGD() {}
/**
* Initialization for parameters
*/
void InitPara(int learner_typeid, ArrayList<double> &param_) {
// 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<NZ_entry> &w_out) {
index_t w_size = w_.size();
index_t ct_nz = 0;
w_out.Init(w_nnz_);
for (index_t i=0; i<w_size; i++) {
if (fabs(w_[i].value) >= 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<index_t> &dataset_index, ArrayList<double> &coef, ArrayList<bool> &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<typename TKernel>
void TGD<TKernel>::LearnersInit_(int learner_typeid) {
index_t i;
learner_typeid_ = learner_typeid;
// init w, w+, w-
w_.Init(n_features_);
for (i=0; i<n_features_; i++) {
w_[i].index = i;
w_[i].value = 0.0; // TODO
}
y_.Init(n_data_);
for (i = 0; i < n_data_; i++) {
y_[i] = (dataset_->y)[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<typename TKernel>
void TGD<TKernel>::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; epo<n_epochs_; epo++) {
/* To mimic the online learning senario, in each epoch,
we randomly permutate the training set, indexed by old_from_new_ */
for (i=0; i<n_data_; i++) {
old_from_new_[i] = i;
}
for (i=0; i<n_data_; i++) {
j = rand() % n_data_;
swap(old_from_new_[i], old_from_new_[j]);
}
ct = 0;
while (ct <= n_iter_) {
work_idx_old = old_from_new_[ct % n_data_];
// eta_ = 1.0 / t_;
eta_ = 1.0 / sqrt(t_);
//scale_w_ = scale_w_ / w_sum;
NZ_entry *xt;
xt = (dataset_->x)[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; i<n_features_; i++) {
if (i < xt->index) {
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<w_.size(); i++)
printf("w[%d]=%f\n", i, w_[i].value);
*/
t_ += 1.0;
ct ++;
ct_k ++;
}
}// for epo
/*
// find max(abs(w_i))
double wi_abs_max = -INFINITY;
double wi_abs;
for (i=0; i<w_.size(); i++) {
wi_abs = fabs(w_[i].value);
if (wi_abs > 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<w_.size(); i++) {
if ( fabs(w_[i].value) > 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<w_size; ct++) {
if (fabs(w_[ct].value) >= 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; i<n_data_; i++) {
hinge_loss = 1- y_[i] * SparseDot(w_, (dataset_->x)[i]);
if (hinge_loss > 0) {
loss_sum += hinge_loss;
}
}
for (i=0; i<n_features_; i++) {
w_sum += fabs(w_[i].value);
}
double obj_value = loss_sum + g_ * w_sum;
printf("Primal objective value: %lf\n", obj_value);
}
}
#endif
fastlib/trunk/contrib/houyang/regmin/regmin.cc
+341
View File
@@ -0,0 +1,341 @@
/**
* @author Hua Ouyang
*
* @file regmin.cc
*
* This file contains main routines for performing Regularized Risk Minimizations.
* Sparse matrix vector manipulations are implemented.
*
* It provides four modes:
* "cv": cross validation;
* "train": model training
* "train_test": training and then online batch testing;
* "test": offline batch testing.
*
* Please refer to README for detail description of usage and examples.
*
* @see regmin.h
* @see opt_smo.h
* @see opt_sgd.h
*/
#include <errno.h>
#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<dataset.n_points; i++) {
inst_max_index = -1;
ReadLine(fp);
dataset.x[i] = &nz_pool[j];
label = strtok(line, " \t");
dataset.y[i] = strtod(label, &endptr);
if (endptr == label) {
printf("No values found for a data point at line %d\n", i+1);
return 0;
}
while (1) {
index = strtok(NULL, ":");
value = strtok(NULL, " \t");
if (value == NULL) {
break;
}
errno = 0;
// in svmlight's data format, feature index begins from 1, not 0
nz_pool[j].index = (index_t) strtol(index, &endptr, 10) - 1;
if (endptr == index || errno !=0 || *endptr !='\0' || nz_pool[j].index<=inst_max_index) {
printf("No values found for a data point at line %d\n", i+1);
return 0;
}
else {
inst_max_index = nz_pool[j].index;
}
errno = 0;
nz_pool[j].value = strtod(value, &endptr);
if ( endptr == value || errno !=0 || (*endptr!='\0' && !isspace(*endptr)) ) {
printf("No values found for a data point at line %d\n", i+1);
return 0;
}
j++;
}
if (inst_max_index > 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<SVMLinearKernel> 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<SVMRBFKernel> 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<SVMLinearKernel> 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<SVMRBFKernel> 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);
}
fastlib/trunk/contrib/houyang/regmin/regmin.h
+1100
View File
File diff suppressed because it is too large Load Diff
fastlib/trunk/contrib/houyang/regmin/regmin_data.h
+440
View File
@@ -0,0 +1,440 @@
/**
* @author Hua Ouyang
*
* @file regmin_data.h
*
* This head file contains sparse data related functions
*
*/
#ifndef U_REGMIN_DATA_H
#define U_REGMIN_DATA_H
#include "fastlib/fastlib.h"
#define Malloc(type,n) (type *)malloc((n)*sizeof(type))
#define ID_LINEAR 0
#define ID_GAUSSIAN 1
/**
* An nonzero entry(dimension) of a data point
*/
struct NZ_entry {
index_t index;
double value;
};
/**
* Sparse labeled dataset
*/
struct Dataset_sl {
index_t n_features;
index_t n_points;
index_t n_classes;
double *y; // labels of data points
struct NZ_entry **x; // data points
//ArrayList<NZ_entry> *x; // TODO
};
/**
* Class for Linear Kernel
*/
class SVMLinearKernel {
public:
// Init of kernel parameters
ArrayList<double> 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<double> 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<NZ_entry> &w) {
if (scale == 0) {
w.ShrinkTo(0);
}
else {
for (index_t i=0; i<w.size(); i++) {
w[i].value = w[i].value * scale;
}
}
}
/**
* Sparse dot product: return w^T x
*/
double SparseDot(ArrayList<NZ_entry> &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 (ct<w_nz && x->index !=-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<NZ_entry> &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 (ct<w.size() || x->index !=-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<w.size(); ct++) {
if (fabs(w[ct].value) < 1.0e-5) { // TODO: thresholding
w.Remove(ct);
ct--;
}
}
}
}
/**
* Sparse vector subtraction: z <= y-x
*/
void SparseSub(ArrayList<NZ_entry> &x,ArrayList<NZ_entry> &y, ArrayList<NZ_entry> &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<y.size(); i++) {
printf("y[%d].index=%d,.value=%f\n", i, y[i].index, y[i].value);
}
for (index_t i=0; i<x.size(); i++) {
printf("x[%d].index=%d,.value=%f\n", i, x[i].index, x[i].value);
}
*/
z.ShrinkTo(0);
if (y_size == 0) { // y: all zeros
z.GrowTo(x_size);
for (i=0; i<x_size; i++) {
z[i].index = x[i].index;
z[i].value = - x[i].value;
}
}
else if (x_size == 0) { // x: all zeros
z.GrowTo(y_size);
for (i=0; i<y_size; i++) {
z[i].index = y[i].index;
z[i].value = y[i].value;
}
}
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
nz_tmp.index = y[ct_y].index;
nz_tmp.value = y[ct_y].value;
z.PushBackCopy(nz_tmp);
++ct_y;
}
else if (ct_y == y_size) { // y reaches end, while x still not
nz_tmp.index = x[ct_x].index;
nz_tmp.value = - x[ct_x].value;
z.PushBackCopy(nz_tmp);
++ct_x;
}
else { // neither x nor y reaches end
if (x[ct_x].index == y[ct_y].index) {
nz_tmp.index = x[ct_x].index;
nz_tmp.value = y[ct_y].value - x[ct_x].value;
z.PushBackCopy(nz_tmp);
++ct_x;
++ct_y;
}
else if (y[ct_y].index > 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<z.size(); i++) {
printf("z[%d].index=%d,.value=%f\n", i, z[i].index, z[i].value);
}
*/
}
/**
* Sparse vector subtraction: y <= y-x
*/
void SparseSubOverwrite(ArrayList<NZ_entry> &x,ArrayList<NZ_entry> &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_size; i++) {
y[i].index = x[i].index;
y[i].value = - x[i].value;
}
}
else if (x_size == 0) { // x: all zeros
// y remains unchanged
}
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
break;
}
else if (ct_y == y.size()) { // y reaches end, while x still not
nz_tmp.index = x[ct_x].index;
nz_tmp.value = - x[ct_x].value;
y.PushBackCopy(nz_tmp);
++ct_x;
++ct_y;
}
else { // neither x nor y reaches end
if (y[ct_y].index == x[ct_x].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) {
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<NZ_entry> &y, ArrayList<NZ_entry> &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