bsvm.cpp

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

		for (j=0;j<nr_class;j++)
		{
			start1[i*nr_class+j] = p;
			start2[i*nr_class+j] = l;
			if (i != j)
				for (k=0;k<count[j];k++)
				{
					yy[p] = i;
					real_i[p] = q;
					active_set[p] = p;
					p++;
					q++;
				}
			else
				q += count[j];
		}
	}
	start1[nr_class*nr_class] = start2[nr_class*nr_class] = l;
}

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

	if(active_size == l) return;

	int i, j;
	for(i=active_size;i<l;i++)
		G[i] = G_bar[i] + lin;
	
	for(i=0;i<active_size;i++)
		if(is_free(i))
		{
			const Qfloat *Q_i = Q->get_Q(real_i[i],real_l);
			double alpha_i = alpha[i], t;
			int y_i = y[i], yy_i = yy[i], ub, k;
			
			t = 2*alpha_i;
			ub = start2[yy_i*nr_class+y_i+1];
			for (j=start2[yy_i*nr_class+y_i];j<ub;j++)
				G[j] += t*Q_i[real_i[j]];
			ub = start2[y_i*nr_class+yy_i+1];	
			for (j=start2[y_i*nr_class+yy_i];j<ub;j++)
				G[j] -= t*Q_i[real_i[j]];
					
			for (k=0;k<nr_class;k++)
				if (k != y_i && k != yy_i)
				{
					ub = start2[k*nr_class+y_i+1];
					for (j=start2[k*nr_class+y_i];j<ub;j++)
						G[j] += alpha_i*Q_i[real_i[j]];
					ub = start2[yy_i*nr_class+k+1];
					for (j=start2[yy_i*nr_class+k];j<ub;j++)
						G[j] += alpha_i*Q_i[real_i[j]];
							
					ub = start2[y_i*nr_class+k+1];
					for (j=start2[y_i*nr_class+k];j<ub;j++)
						G[j] -= alpha_i*Q_i[real_i[j]];
					ub = start2[k*nr_class+yy_i+1];
					for (j=start2[k*nr_class+yy_i];j<ub;j++)
						G[j] -= alpha_i*Q_i[real_i[j]];	
				}			
		}
}

void Solver_MB::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)
{
	this->l = l;
	this->nr_class = nr_class;
	this->real_l = l/(nr_class - 1);
	this->Q = &Q;
	this->lin = lin;
	clone(y,y_,l);
	clone(alpha,alpha_,l);
	C = C_;
	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];
	active_size = l;
	yy = new short[l];
	real_i = new int[l];
	start1 = new int[nr_class*nr_class+1];
	start2 = new int[nr_class*nr_class+1];

	initial_index_table(count);

	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] = lin;
			G_bar[i] = 0;
		}
		
		for (i=0;i<l;i++)
			if (!is_lower_bound(i))
			{
				Qfloat *Q_i = Q.get_Q(real_i[i], real_l);
				double alpha_i = alpha[i];
				double C_i = get_C(i);
				int y_i = y[i], yy_i = yy[i], ub, j, k;
				
				ub = start1[yy_i*nr_class+y_i+1];
				for (j=start1[yy_i*nr_class+y_i];j<ub;j++)
					G[j] += alpha_i*Q_i[real_i[j]];
				if (shrinking && is_upper_bound(i)) 
					for (j=start1[yy_i*nr_class+y_i];j<ub;j++)
						G_bar[j] += C_i*Q_i[real_i[j]];				
					
				ub = start1[y_i*nr_class+yy_i+1];	
				for (j=start1[y_i*nr_class+yy_i];j<ub;j++)
					G[j] -= alpha_i*Q_i[real_i[j]];
				if (shrinking && is_upper_bound(i)) 
					for (j=start1[y_i*nr_class+yy_i];j<ub;j++)
						G_bar[j] += C_i*Q_i[real_i[j]];					
					
				for (k=0;k<nr_class;k++)
					if (k != y_i && k != yy_i)
					{
						ub = start1[k*nr_class+y_i+1];
						for (j=start1[k*nr_class+y_i];j<ub;j++)
							G[j] += alpha_i*Q_i[real_i[j]];
						if (shrinking && is_upper_bound(i)) 
							for (j=start1[k*nr_class+y_i];j<ub;j++)
								G_bar[j] += C_i*Q_i[real_i[j]];							
						ub = start1[yy_i*nr_class+k+1];
						for (j=start1[yy_i*nr_class+k];j<ub;j++)
							G[j] += alpha_i*Q_i[real_i[j]];
						if (shrinking && is_upper_bound(i)) 
							for (j=start1[yy_i*nr_class+k];j<ub;j++)
								G_bar[j] += C_i*Q_i[real_i[j]];							
							
						ub = start1[y_i*nr_class+k+1];
						for (j=start1[y_i*nr_class+k];j<ub;j++)
							G[j] -= alpha_i*Q_i[real_i[j]];
						if (shrinking && is_upper_bound(i)) 
							for (j=start1[y_i*nr_class+k];j<ub;j++)
								G_bar[j] += C_i*Q_i[real_i[j]];							
						ub = start1[k*nr_class+yy_i+1];
						for (j=start1[k*nr_class+yy_i];j<ub;j++)
							G[j] -= alpha_i*Q_i[real_i[j]];
						if (shrinking && is_upper_bound(i)) 
							for (j=start1[k*nr_class+yy_i];j<ub;j++)
								G_bar[j] += C_i*Q_i[real_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
			
			short *y0;
			clone(y0,y,l);
			for (i=0;i<l;i++)
				y[active_set[i]] = y0[i];
			delete[] y0;

			char *alpha_status0;
			clone(alpha_status0,alpha_status,l);
			for (i=0;i<l;i++)
				alpha_status[active_set[i]] = alpha_status0[i];
			delete[] alpha_status0;

			double *alpha0;
			clone(alpha0,alpha,l);
			for (i=0;i<l;i++)
				alpha[active_set[i]] = alpha0[i];
			delete[] alpha0;

			double *G0;
			clone(G0,G,l);
			for (i=0;i<l;i++)
				G[active_set[i]] = G0[i];
			delete[] G0;

			double *G_bar0;
			clone(G_bar0,G_bar,l);
			for (i=0;i<l;i++)
				G_bar[active_set[i]] = G_bar0[i];
			delete[] G_bar0;

			initial_index_table(count);
		}

		++iter;	

		// construct subproblem
		Qfloat **QB;
		QB = new Qfloat *[q];
		for (i=0;i<q;i++)
			QB[i] = Q.get_Q(real_i[working_set[i]], real_l);
		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], y_Bi = y[Bi], yy_Bi = yy[Bi];
			qp.x[i] = alpha[Bi];
			qp.C[i] = get_C(Bi);
			qp.Q[i*qp.n+i] = yyy(y_Bi, yy_Bi, y_Bi, yy_Bi)*
			QB[i][real_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] = 
				yyy(y_Bi, yy_Bi, y[Bj], yy[Bj])*QB[i][real_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++)
		{
			int Bi = working_set[i];
			double d = qp.x[i] - alpha[working_set[i]];
			if(fabs(d) > 1e-12)
			{
				alpha[Bi] = qp.x[i];
				Qfloat *QB_i = QB[i];
				int y_Bi = y[Bi], yy_Bi = yy[Bi], ub, k;

				double t = 2*d;
				ub = start1[yy_Bi*nr_class+y_Bi+1];
				for (j=start1[yy_Bi*nr_class+y_Bi];j<ub;j++)
					G[j] += t*QB_i[real_i[j]];
				ub = start1[y_Bi*nr_class+yy_Bi+1];	
				for (j=start1[y_Bi*nr_class+yy_Bi];j<ub;j++)
					G[j] -= t*QB_i[real_i[j]];
					
				for (k=0;k<nr_class;k++)
					if (k != y_Bi && k != yy_Bi)
					{
						ub = start1[k*nr_class+y_Bi+1];
						for (j=start1[k*nr_class+y_Bi];j<ub;j++)
							G[j] += d*QB_i[real_i[j]];
						ub = start1[yy_Bi*nr_class+k+1];
						for (j=start1[yy_Bi*nr_class+k];j<ub;j++)
							G[j] += d*QB_i[real_i[j]];
							
						ub = start1[y_Bi*nr_class+k+1];
						for (j=start1[y_Bi*nr_class+k];j<ub;j++)
							G[j] -= d*QB_i[real_i[j]];
						ub = start1[k*nr_class+yy_Bi+1];
						for (j=start1[k*nr_class+yy_Bi];j<ub;j++)
							G[j] -= d*QB_i[real_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 = QB[i];
				double C_i = qp.C[i], t = 2*C_i;
				int ub, y_Bi = y[Bi], yy_Bi = yy[Bi], k;
				if (u)
				{
					ub = start1[yy_Bi*nr_class+y_Bi+1];
					for (j=start1[yy_Bi*nr_class+y_Bi];j<ub;j++)
						G_bar[j] -= t*QB_i[real_i[j]];
					ub = start1[y_Bi*nr_class+yy_Bi+1];	
					for (j=start1[y_Bi*nr_class+yy_Bi];j<ub;j++)
						G_bar[j] += t*QB_i[real_i[j]];

					ub = start2[yy_Bi*nr_class+y_Bi+1];
					for (j=start2[yy_Bi*nr_class+y_Bi];j<ub;j++)
						G_bar[j] -= t*QB_i[real_i[j]];
					ub = start2[y_Bi*nr_class+yy_Bi+1];	
					for (j=start2[y_Bi*nr_class+yy_Bi];j<ub;j++)
						G_bar[j] += t*QB_i[real_i[j]];
					
					for (k=0;k<nr_class;k++)
						if (k != y_Bi && k != yy_Bi)
						{
							ub = start1[k*nr_class+y_Bi+1];
							for (j=start1[k*nr_class+y_Bi];j<ub;j++)
								G_bar[j] -= C_i*QB_i[real_i[j]];
							ub = start1[yy_Bi*nr_class+k+1];
							for (j=start1[yy_Bi*nr_class+k];j<ub;j++)
								G_bar[j] -= C_i*QB_i[real_i[j]];
							
							ub = start1[y_Bi*nr_class+k+1];
							for (j=start1[y_Bi*nr_class+k];j<ub;j++)
								G_bar[j] += C_i*QB_i[real_i[j]];
							ub = start1[k*nr_class+yy_Bi+1];
							for (j=start1[k*nr_class+yy_Bi];j<ub;j++)
								G_bar[j] += C_i*QB_i[real_i[j]];

							ub = start2[k*nr_class+y_Bi+1];
							for (j=start2[k*nr_class+y_Bi];j<ub;j++)
								G_bar[j] -= C_i*QB_i[real_i[j]];
							ub = start2[yy_Bi*nr_class+k+1];
							for (j=start2[yy_Bi*nr_class+k];j<ub;j++)
								G_bar[j] -= C_i*QB_i[real_i[j]];
							
							ub = start2[y_Bi*nr_class+k+1];
							for (j=start2[y_Bi*nr_class+k];j<ub;j++)
								G_bar[j] += C_i*QB_i[real_i[j]];
							ub = start2[k*nr_class+yy_Bi+1];
							for (j=start2[k*nr_class+yy_Bi];j<ub;j++)
								G_bar[j] += C_i*QB_i[real_i[j]];
						}
				}
				else
				{
					ub = start1[yy_Bi*nr_class+y_Bi+1];
					for (j=start1[yy_Bi*nr_class+y_Bi];j<ub;j++)
						G_bar[j] += t*QB_i[real_i[j]];
					ub = start1[y_Bi*nr_class+yy_Bi+1];	
					for (j=start1[y_Bi*nr_class+yy_Bi];j<ub;j++)
						G_bar[j] -= t*QB_i[real_i[j]];

					ub = start2[yy_Bi*nr_class+y_Bi+1];
					for (j=start2[yy_Bi*nr_class+y_Bi];j<ub;j++)
						G_bar[j] += t*QB_i[real_i[j]];
					ub = start2[y_Bi*nr_class+yy_Bi+1];	
					for (j=start2[y_Bi*nr_class+yy_Bi];j<ub;j++)
						G_bar[j] -= t*QB_i[real_i[j]];
					
					for (k=0;k<nr_class;k++)
						if (k != y_Bi && k != yy_Bi)
						{
							ub = start1[k*nr_class+y_Bi+1];
							for (j=start1[k*nr_class+y_Bi];j<ub;j++)
								G_bar[j] += C_i*QB_i[real_i[j]];
							ub = start1[yy_Bi*nr_class+k+1];
							for (j=start1[yy_Bi*nr_class+k];j<ub;j++)
								G_bar[j] += C_i*QB_i[real_i[j]];
							
							ub = start1[y_Bi*nr_class+k+1];
							for (j=start1[y_Bi*nr_class+k];j<ub;j++)
								G_bar[j] -= C_i*QB_i[real_i[j]];
							ub = start1[k*nr_class+yy_Bi+1];
							for (j=start1[k*nr_class+yy_Bi];j<ub;j++)
								G_bar[j] -= C_i*QB_i[real_i[j]];

							ub = start2[k*nr_class+y_Bi+1];
							for (j=start2[k*nr_class+y_Bi];j<ub;j++)
								G_bar[j] += C_i*QB_i[real_i[j]];
							ub = start2[yy_Bi*nr_class+k+1];
							for (j=start2[yy_Bi*nr_class+k];j<ub;j++)
								G_bar[j] += C_i*QB_i[real_i[j]];
							
							ub = start2[y_Bi*nr_class+k+1];
							for (j=start2[y_Bi*nr_class+k];j<ub;j++)
								G_bar[j] -= C_i*QB_i[real_i[j]];
							ub = start2[k*nr_class+yy_Bi+1];
							for (j=start2[k*nr_class+yy_Bi];j<ub;j++)
								G_bar[j] -= C_i*QB_i[real_i[j]];
						}				
				}
			}
		}

		delete[] QB;
	}

	// calculate objective value
	{
		double v = 0;
		int i;
		for(i=0;i<l;i++)
			v += alpha[i] * (G[i] + lin);
		si->obj = v/4;
	}

	clone(si->upper_bound,C,nr_class);
	//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[] start1;
	delete[] start2;
	delete[] y;
	delete[] yy;
	delete[] real_i;
	delete[] active_set;

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

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

void Solver_MB::shrink_one(int k)
{
	int i, s = yy[k]*nr_class+y[k], t;
	t = nr_class*nr_class;
	for (i=s+1;i<=t;i++)
		start1[i]--;
	for (i=0;i<=s;i++)
		start2[i]--;
	swap_index(k, start1[s+1]);
	for (i=s+1;i<t;i++)
		swap_index(start1[i], start1[i+1]);
	for (i=0;i<s;i++)
		swap_index(start2[i], start2[i+1]);
}

void Solver_MB::unshrink_one(int k)
{
	int i, s = yy[k]*nr_class+y[k], t;
	swap_index(k, start2[s]);
	for (i=s;i>0;i--)
		swap_index(start2[i], start2[i-1]);
	t = s + 1;
	for (i=nr_class*nr_class;i>t;i--)
		swap_index(start1[i], start1[i-1]);
	t = nr_class*nr_class;
	for (i=s+1;i<=t;i++)
		start1[i]++;
	for (i=0;i<=s;i++)
		start2[i]++;
}

//
// Q matrices for various formulations
//
class BSVC_Q: public Kernel
{ 
public:
	BSVC_Q(const svm_problem& prob, const svm_parameter& param, const schar *y_)
	:Kernel(prob.l, prob.x, param)
	{
		clone(y,y_,prob.l);
		cache = new Cache(prob.l,(int)(param.cache_size*(1<<20)));
		ktype = param.kernel_type;
		expr = param.expr;
		rho = param.rho;
	}
	
	Qfloat *get_Q(int i, int len) const
	{
		Qfloat *data;
		int start;
		if(((start = cache->get_data(i,&data,len)) < len) && ktype != 4)
		{
			for(int j=start;j<len;j++)
				data[j] = (Qfloat)y[i]*y[j]*((this->*kernel_function)(i,j) + 1);
		}
		else if(ktype == 4)
		  {
		    SEXP R_fcall, dummy, ans;
		    PROTECT(R_fcall = lang2(expr, R_NilValue));
		    PROTECT(ans= allocSExp(REALSXP));
		    PROTECT(dummy = allocVector(INTSXP, 1));
		    INTEGER(dummy)[0] = i+1;
		    SETCADR(R_fcall,dummy);
		    ans = eval(R_fcall,rho);
		    for(int j=start;j<len;j++)
		      data[j] = (Qfloat)(y[i]*y[j]*REAL(ans)[j]);
		    UNPROTECT(3);
		  }
		return data;
	}

	void swap_index(int i, int j) const
	{
		cache->swap_index(i,j);
		Kernel::swap_index(i,j);
		swap(y[i],y[j]);
	}

	~BSVC_Q()
	{
		delete[] y;
		delete cache;
	}
private:
        SEXP expr;
        SEXP rho;
        int ktype;
	schar *y;
	Cache *cache;
};

class ONE_CLASS_Q: public Kernel
{
public:
	ONE_CLASS_Q(const svm_problem& prob, const svm_parameter& param)
	:Kernel(prob.l, prob.x, param)
	{
		cache = new Cache(prob.l,(int)(param.cache_size*(1<<20)));
	}
	
	Qfloat *get_Q(int i, int len) const
	{
		Qfloat *data;
		int start;
		if((start = cache->get_data(i,&data,len)) < len)
		{
			for(int j=start;j<len;j++)
				data[j] = (Qfloat)(this->*kernel_function)(i,j);
		}
		return data;
	}

	void swap_index(int i, int j) const
	{
		cache->swap_index(i,j);
		Kernel::swap_index(i,j);
	}

	~ONE_CLASS_Q()
	{
		delete cache;
	}
private:
	Cache *cache;
};

class BONE_CLASS_Q: public Kernel
{
public:
	BONE_CLASS_Q(const svm_problem& prob, const svm_parameter& param)
	:Kernel(prob.l, prob.x, param)
	{
		cache = new Cache(prob.l,(int)(param.cache_size*(1<<20)));
	}
	
	Qfloat *get_Q(int i, int len) const
	{
		Qfloat *data;
		int start;
		if((start = cache->get_data(i,&data,len)) < len)
		{
			for(int j=start;j<len;j++)
				data[j] = (Qfloat)(this->*kernel_function)(i,j) + 1;
		}
		return data;
	}

	~BONE_CLASS_Q()
	{
		delete cache;
	}
private:
	Cache *cache;
};

class BSVR_Q: public Kernel
{ 
public:
	BSVR_Q(const svm_problem& prob, const svm_parameter& param)
	:Kernel(prob.l, prob.x, param)
	{
		l = prob.l;
		cache = new Cache(l,(int)(param.cache_size*(1<<20)));
		sign = new schar[2*l];
		index = new int[2*l];
		for(int k=0;k<l;k++)
		{
			sign[k] = 1;
			sign[k+l] = -1;
			index[k] = k;
			index[k+l] = k;
		}
		q = param.qpsize;
		buffer = new Qfloat*[q];
		for (int i=0;i<q;i++)
			buffer[i] = new Qfloat[2*l];	
		next_buffer = 0;
	}

	void swap_index(int i, int j) const
	{
		swap(sign[i],sign[j]);
		swap(index[i],index[j]);
	}
	
	Qfloat *get_Q(int i, int len) const
	{
		Qfloat *data;
		int real_i = index[i];
		if(cache->get_data(real_i,&data,l) < l)
		{
			for(int j=0;j<l;j++)
				data[j] = (Qfloat)(this->*kernel_function)(real_i,j) + 1;
		}

		// reorder and copy
		Qfloat *buf = buffer[next_buffer];
		next_buffer = (next_buffer+1)%q;
		schar si = sign[i];
		for(int j=0;j<len;j++)
			buf[j] = si * sign[j] * data[index[j]];
		return buf;
	}

	~BSVR_Q()
	{
		delete cache;
		delete[] sign;
		delete[] index;
		for (int i=0;i<q;i++)
			delete[] buffer[i];
		delete[] buffer;
	}
private:
	int l, q;
	Cache *cache;
	schar *sign;
	int *index;
	mutable int next_buffer;
	Qfloat** buffer;
};

//
// construct and solve various formulations
//
static void solve_c_svc(
	const svm_problem *prob, const svm_parameter* param,
	double *alpha, Solver_B::SolutionInfo* si, double Cp, double Cn)
{
	int l = prob->l;
	double *minus_ones = new double[l];
	schar *y = new schar[l];

	int i;

	for(i=0;i<l;i++)
	{
		alpha[i] = 0;
		minus_ones[i] = -1;
		if(prob->y[i] > 0) y[i] = +1; else y[i]=-1;
	}
	
	if (param->kernel_type == LINEAR)
	{
		double *w = new double[prob->n+1];
		for (i=0;i<=prob->n;i++)
			w[i] = 0;
		Solver_B_linear s;
		int totaliter = 0;
		double Cpj = param->Cbegin, Cnj = param->Cbegin*Cn/Cp;
		while (Cpj < Cp)
		{
			totaliter += s.Solve(l, prob->x, minus_ones, y, alpha, w, 
			Cpj, Cnj, param->eps, si, param->shrinking, param->qpsize);
			if (Cpj*param->Cstep >= Cp)
			{