bsvm.cpp

这是核学习的一个基础软件包

	void do_shrinking();
private:
	double Cp, Cn;
	double *b;
	schar *y;
};

void Solver_B::swap_index(int i, int j)
{
	Q->swap_index(i,j);
	swap(y[i],y[j]);
	swap(G[i],G[j]);
	swap(alpha_status[i],alpha_status[j]);
	swap(alpha[i],alpha[j]);
	swap(b[i],b[j]);
	swap(active_set[i],active_set[j]);
	swap(G_bar[i],G_bar[j]);
}

void Solver_B::reconstruct_gradient()
{
	// reconstruct inactive elements of G from G_bar and free variables

	if(active_size == l) return;

	int i;
	for(i=active_size;i<l;i++)
		G[i] = G_bar[i] + b[i];
	
	for(i=0;i<active_size;i++)
		if(is_free(i))
		{
			const Qfloat *Q_i = Q->get_Q(i,l);
			double alpha_i = alpha[i];
			for(int j=active_size;j<l;j++)
				G[j] += alpha_i * Q_i[j];
		}
}

void Solver_B::Solve(int l, const Kernel& Q, double *b_, schar *y_,
	double *alpha_, double Cp, double Cn, double eps, SolutionInfo* si,
	int shrinking, int qpsize)
{
	this->l = l;
	this->Q = &Q;
	b = b_;
	clone(y, y_, l);
	clone(alpha,alpha_,l);
	this->Cp = Cp;
	this->Cn = Cn;
	this->eps = eps;
	this->qpsize = qpsize;
	unshrinked = false;

	// initialize alpha_status
	{
		alpha_status = new char[l];
		for(int i=0;i<l;i++)
			update_alpha_status(i);
	}

	// initialize active set (for shrinking)
	{
		active_set = new int[l];
		for(int i=0;i<l;i++)
			active_set[i] = i;
		active_size = l;
	}

	BQP qp;
	working_set = new int[qpsize];
	qp.eps = eps/10;
	qp.C = new double[qpsize];
	qp.x = new double[qpsize];
	qp.p = new double[qpsize];
	qp.Q = new double[qpsize*qpsize];

	// initialize gradient
	{
		G = new double[l];
		G_bar = new double[l];
		int i;
		for(i=0;i<l;i++)
		{
			G[i] = b[i];
			G_bar[i] = 0;
		}
		for(i=0;i<l;i++)
			if(!is_lower_bound(i))
			{
				Qfloat *Q_i = Q.get_Q(i,l);
				double C_i = get_C(i);
				double alpha_i = alpha[i];
				int j;
				for(j=0;j<l;j++)
					G[j] += alpha_i*Q_i[j];
				if (shrinking)
					if(is_upper_bound(i))
						for(j=0;j<l;j++)
							G_bar[j] += C_i*Q_i[j];
			}
	}

	// optimization step

	int iter = 0;
	int counter = min(l*2/qpsize,2000/qpsize)+1;

	while(1)
	{
		// show progress and do shrinking

		if(--counter == 0)
		{
			counter = min(l*2/qpsize, 2000/qpsize);
			if(shrinking) do_shrinking();
			//info(".");
		}

		int i,j,q;
		if (select_working_set(q) < eps)
		{
			// reconstruct the whole gradient
			reconstruct_gradient();
			// reset active set size and check
			active_size = l;
			//info("*");info_flush();
			if (select_working_set(q) < eps)
				break;
			else
				counter = 1;	// do shrinking next iteration
		}
		
		++iter;

		// construct subproblem
		Qfloat **QB;
		QB = new Qfloat *[q];
		for (i=0;i<q;i++)
			QB[i] = Q.get_Q(working_set[i], active_size);
		qp.n = q;
		for (i=0;i<qp.n;i++)
			qp.p[i] = G[working_set[i]];
		for (i=0;i<qp.n;i++)
		{
			int Bi = working_set[i];
			qp.x[i] = alpha[Bi];
			qp.C[i] = get_C(Bi);
			qp.Q[i*qp.n+i] = QB[i][Bi];
			qp.p[i] -= qp.Q[i*qp.n+i]*alpha[Bi];
			for (j=i+1;j<qp.n;j++)
			{			
				int Bj = working_set[j];
				qp.Q[i*qp.n+j] = qp.Q[j*qp.n+i] = QB[i][Bj];
				qp.p[i] -= qp.Q[i*qp.n+j]*alpha[Bj];
				qp.p[j] -= qp.Q[j*qp.n+i]*alpha[Bi];
			}
		}

		solvebqp(&qp);

		// update G

		for(i=0;i<q;i++)
		{
			double d = qp.x[i] - alpha[working_set[i]];
			if(fabs(d)>1e-12)
			{
				alpha[working_set[i]] = qp.x[i];
				Qfloat *QB_i = QB[i];
				for(j=0;j<active_size;j++)
					G[j] += d*QB_i[j];
			}
		}

		// update alpha_status and G_bar

		for (i=0;i<q;i++)
		{
			int Bi = working_set[i];
			bool u = is_upper_bound(Bi);
			update_alpha_status(Bi);
			if (!shrinking)
				continue;
			if (u != is_upper_bound(Bi))
			{
				Qfloat *QB_i = Q.get_Q(Bi, l);
				double C_i = qp.C[i];
				if (u)
					for (j=0;j<l;j++)
						G_bar[j] -= C_i*QB_i[j];
				else
					for (j=0;j<l;j++)
						G_bar[j] += C_i*QB_i[j];
			}
		}

		delete[] QB;
	}

	// calculate objective value
	{
		double v = 0;
		int i;
		for(i=0;i<l;i++)
			v += alpha[i] * (G[i] + b[i]);

		si->obj = v/2;
	}

	// juggle everything back
	/*{
		for(int i=0;i<l;i++)
			while(active_set[i] != i)
				swap_index(i,active_set[i]);
				// or Q.swap_index(i,active_set[i]);
	}*/

	si->upper_bound = new double[2];
	si->upper_bound[0] = Cp;
	si->upper_bound[1] = Cn;

	//	info("\noptimization finished, #iter = %d\n",iter);

	// put back the solution
	{
		for(int i=0;i<l;i++)
			alpha_[active_set[i]] = alpha[i];
	}

	delete[] active_set;
	delete[] alpha;
	delete[] alpha_status;
	delete[] G;
	delete[] G_bar;
	delete[] y;

	delete[] working_set;
	delete[] qp.p;
	delete[] qp.C;
	delete[] qp.x;
	delete[] qp.Q;
}

// return maximal violation
double Solver_B::select_working_set(int &q)
{
	int i, j, q_2 = qpsize/2;
	double maxvio = 0, max0;
	double *positive_max;
	int *positive_set;

	positive_max = new double[qpsize];
	positive_set = new int[qpsize];
	q = 0;

	for (i=0;i<q_2;i++)
		positive_max[i] = INF/2;
	for (i=0;i<active_size;i++)
	{	
		if(!is_free(i)) continue;
		double v = fabs(G[i]);
		if(v < positive_max[0])
		{
			for (j=1;j<q_2;j++)
			{
				if (v >= positive_max[j])
					break;
				positive_max[j-1] = positive_max[j];
				positive_set[j-1] = positive_set[j];
			}
			positive_max[j-1] = v;
			positive_set[j-1] = i;
		}
	} 
	for (i=0;i<q_2;i++)
		if (positive_max[i] != INF/2)
			working_set[q++] = positive_set[i];
	max0 = q ? positive_max[0] : 0;
	q_2 = qpsize - q;
				
	for (i=0;i<q_2;i++)
		positive_max[i] = -INF;
	for (i=0;i<active_size;i++)
	{
		double v = fabs(G[i]);
		if (is_free(i) && v <= max0) continue;
		if(is_upper_bound(i))
		{
			if(G[i]<0) continue;
		}
		else if(is_lower_bound(i))
		{
			if(G[i]>0) continue;
		}
		if (v > positive_max[0])
		{
			for (j=1;j<q_2;j++)
			{
				if (v <= positive_max[j])
					break;
				positive_max[j-1] = positive_max[j];
				positive_set[j-1] = positive_set[j];
			}
			positive_max[j-1] = v;
			positive_set[j-1] = i;
		}
	}
	for (i=0;i<q_2;i++)
		if (positive_max[i] != -INF)
		{
			working_set[q++] = positive_set[i];
			maxvio = max(maxvio,positive_max[i]);
		}

	delete[] positive_set;
	delete[] positive_max;
	return maxvio;
}

void Solver_B::shrink_one(int k)
{
	swap_index(k, active_size);
}

void Solver_B::unshrink_one(int k)
{
	swap_index(k, active_size);
}

void Solver_B::do_shrinking()
{
	int k;

	double Gm = select_working_set(k);
	if (Gm < eps)
		return;

	// shrink
	
	for(k=0;k<active_size;k++)
	{
		if (is_lower_bound(k))
		{
			if (G[k] <= Gm)
				continue;
		}
		else
			if (is_upper_bound(k))
			{
				if (G[k] >= -Gm)
					continue;
			}
			else
				continue;

		--active_size;
		shrink_one(k);
		--k;	// look at the newcomer
	}

	// unshrink, check all variables again before final iterations

	if (unshrinked || Gm > eps*10)
		return;
	
	unshrinked = true;
	reconstruct_gradient();

	for(k=l-1;k>=active_size;k--)
	{
		if (is_lower_bound(k))
		{
			if (G[k] > Gm)
				continue;
		}
		else
			if (is_upper_bound(k))
			{
				if (G[k] < -Gm)
					continue;
			}
			else
				continue;

		unshrink_one(k);
		active_size++;
		++k;	// look at the newcomer
	}
}

class Solver_B_linear : public Solver_B
{
public:
	Solver_B_linear() {};
	~Solver_B_linear() {};
	int Solve(int l, svm_node * const * x_, double *b_, schar *y_,
	double *alpha_, double *w, double Cp, double Cn, double eps, SolutionInfo* si, 
	int shrinking, int qpsize);
private:
	double get_C(int i)
	{
		return (y[i] > 0)? Cp : Cn;
	}
	void swap_index(int i, int j);
	void reconstruct_gradient();
	double dot(int i, int j);
	double Cp, Cn;
	double *b;
	schar *y;
	double *w;
	const svm_node **x;
};

double Solver_B_linear::dot(int i, int j)
{
	const svm_node *px = x[i], *py = x[j];	
	double sum = 0;
	while(px->index != -1 && py->index != -1)
	{
		if(px->index == py->index)
		{
			sum += px->value * py->value;
			++px;
			++py;
		}
		else
		{
			if(px->index > py->index)
				++py;
			else
				++px;
		}			
	}
	return sum;
}

void Solver_B_linear::swap_index(int i, int j)
{
	swap(y[i],y[j]);
	swap(G[i],G[j]);
	swap(alpha_status[i],alpha_status[j]);
	swap(alpha[i],alpha[j]);
	swap(b[i],b[j]);
	swap(active_set[i],active_set[j]);
	swap(x[i], x[j]);
}

void Solver_B_linear::reconstruct_gradient()
{
	int i;
	for(i=active_size;i<l;i++)
	{
		double sum = 0;
		for (const svm_node *px = x[i];px->index != -1;px++)
			sum += w[px->index]*px->value;
		sum += w[0];
		G[i] = y[i]*sum + b[i];
	}
}

int Solver_B_linear::Solve(int l, svm_node * const * x_, double *b_, schar *y_,
	double *alpha_, double *w, double Cp, double Cn, double eps, SolutionInfo* si,
	int shrinking, int qpsize)
{
	this->l = l;
	clone(x, x_, l);	
	clone(b, b_, l);
	clone(y, y_, l);
	clone(alpha,alpha_,l);
	this->Cp = Cp;
	this->Cn = Cn;
	this->eps = eps;
	this->qpsize = qpsize;
	this->w = w;
	unshrinked = false;

	// initialize alpha_status
	{
		alpha_status = new char[l];
		for(int i=0;i<l;i++)
			update_alpha_status(i);
	}

	// initialize active set (for shrinking)
	{
		active_set = new int[l];
		for(int i=0;i<l;i++)
			active_set[i] = i;
		active_size = l;
	}

	BQP qp;
	working_set = new int[qpsize];
	qp.eps = eps/100;
	qp.C = new double[qpsize];
	qp.x = new double[qpsize];
	qp.p = new double[qpsize];
	qp.Q = new double[qpsize*qpsize];

	// initialize gradient
	{
		G = new double[l];
		int i;
		bool allzero = true;
		for(i=0;i<l;i++)
		{
			G[i] = b[i];
			if(!is_lower_bound(i))
				allzero = false;
		}
		if (!allzero)
			for(i=0;i<l;i++)
			{
				double sum = 0;
				for (const svm_node *px = x[i];px->index != -1;px++)
					sum += w[px->index]*px->value;
				sum += w[0];
				G[i] += y[i]*sum;
			}			
	}

	// optimization step

	int iter = 0;
	int counter = min(l*2/qpsize,2000/qpsize)+1;

	while(1)
	{
		// show progress and do shrinking

		if(--counter == 0)
		{
			counter = min(l*2/qpsize, 2000/qpsize);
			if(shrinking) do_shrinking();
			//	info(".");
		}

		int i,j,q;
		if (select_working_set(q) < eps)
		{
			// reconstruct the whole gradient
			reconstruct_gradient();
			// reset active set size and check
			active_size = l;
			//	info("*");info_flush();
			if (select_working_set(q) < eps)
				break;
			else
				counter = 1;	// do shrinking next iteration
		}
		
		++iter;

		// construct subproblem
		qp.n = q;
		for (i=0;i<qp.n;i++)
			qp.p[i] = G[working_set[i]];
		for (i=0;i<qp.n;i++)
		{
			int Bi = working_set[i];
			qp.x[i] = alpha[Bi];
			qp.C[i] = get_C(Bi);
			qp.Q[i*qp.n+i] = dot(Bi, Bi) + 1;
			qp.p[i] -= qp.Q[i*qp.n+i]*alpha[Bi];
			for (j=i+1;j<qp.n;j++)
			{			
				int Bj = working_set[j];
				qp.Q[i*qp.n+j] = qp.Q[j*qp.n+i] = y[Bi]*y[Bj]*(dot(Bi, Bj) + 1);
				qp.p[i] -= qp.Q[i*qp.n+j]*alpha[Bj];
				qp.p[j] -= qp.Q[j*qp.n+i]*alpha[Bi];
			}
		}

		solvebqp(&qp);

		// update G

		for(i=0;i<q;i++)
		{
			int Bi = working_set[i];
			double d = qp.x[i] - alpha[Bi];
			if(fabs(d)>1e-12)
			{
				alpha[Bi] = qp.x[i];
				update_alpha_status(Bi);
				double yalpha = y[Bi]*d;
				for (const svm_node *px = x[Bi];px->index != -1;px++)
					w[px->index] += yalpha*px->value;
				w[0] += yalpha;
			}
		}
		for(j=0;j<active_size;j++)
		{
			double sum = 0;
			for (const svm_node *px = x[j];px->index != -1;px++)
				sum += w[px->index]*px->value;
			sum += w[0];
			G[j] = y[j]*sum + b[j];
		}

	}

	// calculate objective value
	{
		double v = 0;
		int i;
		for(i=0;i<l;i++)
			v += alpha[i] * (G[i] + b[i]);

		si->obj = v/2;
	}

	// juggle everything back
	/*{
		for(int i=0;i<l;i++)
			while(active_set[i] != i)
				swap_index(i,active_set[i]);
				// or Q.swap_index(i,active_set[i]);
	}*/

	si->upper_bound = new double[2];
	si->upper_bound[0] = Cp;
	si->upper_bound[1] = Cn;

	// info("\noptimization finished, #iter = %d\n",iter);

	// put back the solution
	{
		for(int i=0;i<l;i++)
			alpha_[active_set[i]] = alpha[i];
	}

	delete[] active_set;
	delete[] alpha;
	delete[] alpha_status;
	delete[] G;
	delete[] y;
	delete[] b;

	delete[] working_set;
	delete[] qp.p;
	delete[] qp.C;
	delete[] qp.x;
	delete[] qp.Q;

	return iter;
}


class Solver_MB : public Solver_B
{
public:
	Solver_MB() {};
	~Solver_MB() {};
	void Solve(int l, const Kernel& Q, double lin, double *alpha_,
	short *y_, double *C, double eps, SolutionInfo* si,
	int shrinking, int qpsize, int nr_class, int *count);
private:
	short *y, *yy;
	double *C;
	double lin;
	int *real_i;
	int real_l;

	int nr_class;
	int *start1, *start2;

	double get_C(int i)
	{
		return C[y[i]];
	}
	void swap_index(int i, int j);
	void reconstruct_gradient();
	void shrink_one(int k);
	void unshrink_one(int k);
	void initial_index_table(int *);
	int yyy(int yi, int yyi, int yj, int yyj) const
	{
		int xx = 0;
		if (yi == yj)
			xx++;
		if (yyi == yyj)
			xx++;
		if (yi == yyj)
			xx--;
		if (yj == yyi)
			xx--;
		return xx;
	}
};

void Solver_MB::swap_index(int i, int j)
{
	if (i == j)
		return;
	swap(y[i],y[j]);
	swap(yy[i],yy[j]);
	swap(G[i],G[j]);
	swap(alpha_status[i],alpha_status[j]);
	swap(alpha[i],alpha[j]);
	swap(active_set[i],active_set[j]);
	swap(real_i[i], real_i[j]);
	swap(G_bar[i],G_bar[j]);
}

void Solver_MB::initial_index_table(int *count)
{
	int i, j, k, p, q;

	p = 0;
	for (i=0;i<nr_class;i++)
	{
		q = 0;