Algorithm-LibLinear

 view release on metacpan or  search on metacpan

src/liblinear/linear.cpp  view on Meta::CPAN

	double v = 0;
	int nnz = 0;
	for(j=0; j<w_size; j++)
		if(w[j] != 0)
		{
			v += fabs(w[j]);
			nnz++;
		}
	if (regularize_bias == 0)
		v -= fabs(w[w_size-1]);
	for(j=0; j<l; j++)
		if(y[j] == 1)
			v += C[GETI(j)]*log(1+1/exp_wTx[j]);
		else
			v += C[GETI(j)]*log(1+exp_wTx[j]);

	info("Objective value = %lf\n", v);
	info("#nonzeros/#features = %d/%d\n", nnz, w_size);

	delete [] index;
	delete [] y;
	delete [] Hdiag;
	delete [] Grad;
	delete [] wpd;
	delete [] xjneg_sum;
	delete [] xTd;
	delete [] exp_wTx;
	delete [] exp_wTx_new;
	delete [] tau;
	delete [] D;

	return newton_iter;
}

static int compare_feature_node(const void *a, const void *b)
{
	double a_value = (*(feature_node *)a).value;
	double b_value = (*(feature_node *)b).value;
	int a_index = (*(feature_node *)a).index;
	int b_index = (*(feature_node *)b).index;

	if(a_value < b_value)
		return -1;
	else if(a_value == b_value)
	{
		if(a_index < b_index)
			return -1;
		else if(a_index == b_index)
			return 0;
	}
	return 1;
}

// elements before the returned index are < pivot, while those after are >= pivot
static int partition(feature_node *nodes, int low, int high)
{
	int i;
	int index;

	swap(nodes[low + rand()%(high-low+1)], nodes[high]); // select and move pivot to the end

	index = low;
	for(i = low; i < high; i++)
		if (compare_feature_node(&nodes[i], &nodes[high]) == -1)
		{
			swap(nodes[index], nodes[i]);
			index++;
		}

	swap(nodes[high], nodes[index]);
	return index;
}

// rearrange nodes so that nodes[:k] contains nodes with the k smallest values.
static void quick_select_min_k(feature_node *nodes, int low, int high, int k)
{
	int pivot;
	if(low == high)
		return;
	pivot = partition(nodes, low, high);
	if(pivot == k)
		return;
	else if(k-1 < pivot)
		return quick_select_min_k(nodes, low, pivot-1, k);
	else
		return quick_select_min_k(nodes, pivot+1, high, k);
}

// A two-level coordinate descent algorithm for
// a scaled one-class SVM dual problem
//
//  min_\alpha  0.5(\alpha^T Q \alpha),
//    s.t.      0 <= \alpha_i <= 1 and
//              e^T \alpha = \nu l
//
//  where Qij = xi^T xj
//
// Given:
// x, nu
// eps is the stopping tolerance
//
// solution will be put in w and rho
//
// this function returns the number of iterations
//
// See Algorithm 7 in supplementary materials of Chou et al., SDM 2020.

static int solve_oneclass_svm(const problem *prob, const parameter *param, double *w, double *rho)
{
	int l = prob->l;
	int w_size = prob->n;
	double eps = param->eps;
	double nu = param->nu;
	int i, j, s, iter = 0;
	double Gi, Gj;
	double Qij, quad_coef, delta, sum;
	double old_alpha_i;
	double *QD = new double[l];
	double *G = new double[l];
	int *index = new int[l];
	double *alpha = new double[l];
	int max_inner_iter;
	int max_iter = 1000;
	int active_size = l;

	double negGmax;                 // max { -grad(f)_i | i in Iup }
	double negGmin;                 // min { -grad(f)_i | i in Ilow }
	// Iup = { i | alpha_i < 1 }, Ilow = { i | alpha_i > 0 }
	feature_node *max_negG_of_Iup = new feature_node[l];
	feature_node *min_negG_of_Ilow = new feature_node[l];
	feature_node node;

	int n = (int)(nu*l);            // # of alpha's at upper bound
	for(i=0; i<n; i++)
		alpha[i] = 1;
	if (n<l)
		alpha[i] = nu*l-n;
	for(i=n+1; i<l; i++)
		alpha[i] = 0;

	for(i=0; i<w_size; i++)
		w[i] = 0;
	for(i=0; i<l; i++)
	{
		feature_node * const xi = prob->x[i];
		QD[i] = sparse_operator::nrm2_sq(xi);

src/liblinear/linear.cpp  view on Meta::CPAN

		negGmin = INF;

		for (s=0; s<active_size; s++)
		{
			i = index[s];
			feature_node * const xi = prob->x[i];
			G[i] = sparse_operator::dot(w, xi);
			if (alpha[i] < 1)
				negGmax = max(negGmax, -G[i]);
			if (alpha[i] > 0)
				negGmin = min(negGmin, -G[i]);
		}

		if (negGmax - negGmin < eps)
		{
			if (active_size == l)
				break;
			else
			{
				active_size = l;
				info("*");
				continue;
			}
		}

		for(s=0; s<active_size; s++)
		{
			i = index[s];
			if ((alpha[i] == 1 && -G[i] > negGmax) ||
			    (alpha[i] == 0 && -G[i] < negGmin))
			{
				active_size--;
				swap(index[s], index[active_size]);
				s--;
			}
		}

		max_inner_iter = max(active_size/10, 1);
		int len_Iup = 0;
		int len_Ilow = 0;
		for(s=0; s<active_size; s++)
		{
			i = index[s];
			node.index = i;
			node.value = -G[i];

			if (alpha[i] < 1)
			{
				max_negG_of_Iup[len_Iup] = node;
				len_Iup++;
			}

			if (alpha[i] > 0)
			{
				min_negG_of_Ilow[len_Ilow] = node;
				len_Ilow++;
			}
		}
		max_inner_iter = min(max_inner_iter, min(len_Iup, len_Ilow));

		quick_select_min_k(max_negG_of_Iup, 0, len_Iup-1, len_Iup-max_inner_iter);
		qsort(&(max_negG_of_Iup[len_Iup-max_inner_iter]), max_inner_iter, sizeof(struct feature_node), compare_feature_node);

		quick_select_min_k(min_negG_of_Ilow, 0, len_Ilow-1, max_inner_iter);
		qsort(min_negG_of_Ilow, max_inner_iter, sizeof(struct feature_node), compare_feature_node);

		for (s=0; s<max_inner_iter; s++)
		{
			i = max_negG_of_Iup[len_Iup-s-1].index;
			j = min_negG_of_Ilow[s].index;

			if ((alpha[i] == 0 && alpha[j] == 0) ||
			    (alpha[i] == 1 && alpha[j] == 1))
				continue;

			feature_node const * xi = prob->x[i];
			feature_node const * xj = prob->x[j];

			Gi = sparse_operator::dot(w, xi);
			Gj = sparse_operator::dot(w, xj);

			int violating_pair = 0;
			if (alpha[i] < 1 && alpha[j] > 0 && -Gj + 1e-12 < -Gi)
				violating_pair = 1;
			else
				if (alpha[i] > 0 && alpha[j] < 1 && -Gi + 1e-12 < -Gj)
					violating_pair = 1;
			if (violating_pair == 0)
				continue;

			Qij = sparse_operator::sparse_dot(xi, xj);
			quad_coef = QD[i] + QD[j] - 2*Qij;
			if(quad_coef <= 0)
				quad_coef = 1e-12;
			delta = (Gi - Gj) / quad_coef;
			old_alpha_i = alpha[i];
			sum = alpha[i] + alpha[j];
			alpha[i] = alpha[i] - delta;
			alpha[j] = alpha[j] + delta;
			if (sum > 1)
			{
				if (alpha[i] > 1)
				{
					alpha[i] = 1;
					alpha[j] = sum - 1;
				}
			}
			else
			{
				if (alpha[j] < 0)
				{
					alpha[j] = 0;
					alpha[i] = sum;
				}
			}
			if (sum > 1)
			{
				if (alpha[j] > 1)
				{
					alpha[j] = 1;
					alpha[i] = sum - 1;
				}
			}
			else

src/liblinear/linear.cpp  view on Meta::CPAN

			solve_l1r_lr(&prob_col, param, w, Cp, Cn, primal_solver_tol);
			delete [] prob_col.y;
			delete [] prob_col.x;
			delete [] x_space;
			break;
		}
		case L2R_LR_DUAL:
		{
			iter = solve_l2r_lr_dual(prob, param, w, Cp, Cn, dual_solver_max_iter);
			if(iter >= dual_solver_max_iter)
			{
				info("\nWARNING: reaching max number of iterations\nSwitching to use -s 0\n\n");
				// primal_solver_tol obtained from eps for dual may be too loose
				primal_solver_tol *= 0.1;
				l2r_lr_fun fun_obj(prob, param, C);
				NEWTON newton_obj(&fun_obj, primal_solver_tol);
				newton_obj.set_print_string(liblinear_print_string);
				newton_obj.newton(w);
			}
			break;
		}
		case L2R_L2LOSS_SVR:
		{
			l2r_l2_svr_fun fun_obj(prob, param, C);
			NEWTON newton_obj(&fun_obj, primal_solver_tol);
			newton_obj.set_print_string(liblinear_print_string);
			newton_obj.newton(w);
			break;
		}
		case L2R_L1LOSS_SVR_DUAL:
		{
			iter = solve_l2r_l1l2_svr(prob, param, w, dual_solver_max_iter);
			if(iter >= dual_solver_max_iter)
				info("\nWARNING: reaching max number of iterations\nUsing -s 11 may be faster (also see FAQ)\n\n");			

			break;
		}
		case L2R_L2LOSS_SVR_DUAL:
		{
			iter = solve_l2r_l1l2_svr(prob, param, w, dual_solver_max_iter);
			if(iter >= dual_solver_max_iter)
			{
				info("\nWARNING: reaching max number of iterations\nSwitching to use -s 11\n\n");
				// primal_solver_tol obtained from eps for dual may be too loose
				primal_solver_tol *= 0.001;
				l2r_l2_svr_fun fun_obj(prob, param, C);
				NEWTON newton_obj(&fun_obj, primal_solver_tol);
				newton_obj.set_print_string(liblinear_print_string);
				newton_obj.newton(w);
			}
			break;
		}
		default:
			fprintf(stderr, "ERROR: unknown solver_type\n");
			break;
	}

	delete[] C;
}

// Calculate the initial C for parameter selection
static double calc_start_C(const problem *prob, const parameter *param)
{
	int i;
	double xTx, max_xTx;
	max_xTx = 0;
	for(i=0; i<prob->l; i++)
	{
		xTx = 0;
		feature_node *xi=prob->x[i];
		while(xi->index != -1)
		{
			double val = xi->value;
			xTx += val*val;
			xi++;
		}
		if(xTx > max_xTx)
			max_xTx = xTx;
	}

	double min_C = 1.0;
	if(param->solver_type == L2R_LR)
		min_C = 1.0 / (prob->l * max_xTx);
	else if(param->solver_type == L2R_L2LOSS_SVC)
		min_C = 1.0 / (2 * prob->l * max_xTx);
	else if(param->solver_type == L2R_L2LOSS_SVR)
	{
		double sum_y, loss, y_abs;
		double delta2 = 0.1;
		sum_y = 0, loss = 0;
		for(i=0; i<prob->l; i++)
		{
			y_abs = fabs(prob->y[i]);
			sum_y += y_abs;
			loss += max(y_abs - param->p, 0.0) * max(y_abs - param->p, 0.0);
		}
		if(loss > 0)
			min_C = delta2 * delta2 * loss / (8 * sum_y * sum_y * max_xTx);
		else
			min_C = INF;
	}

	return pow( 2, floor(log(min_C) / log(2.0)) );
}

static double calc_max_p(const problem *prob)
{
	int i;
	double max_p = 0.0;
	for(i = 0; i < prob->l; i++)
		max_p = max(max_p, fabs(prob->y[i]));

	return max_p;
}

static void find_parameter_C(const problem *prob, parameter *param_tmp, double start_C, double max_C, double *best_C, double *best_score, const int *fold_start, const int *perm, const problem *subprob, int nr_fold)
{
	// variables for CV
	int i;
	double *target = Malloc(double, prob->l);