bsvm.cpp

svm的实现源码

					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
		}

		if (counter == min(l*2/qpsize, 2000/qpsize))
		{
			bool same = true;
			for (i=0;i<qpsize;i++)
				if (old_working_set[i] != working_set[i]) 
				{
					same = false;
					break;
				}

			if (same)
				break;
		}

		for (i=0;i<qpsize;i++)
			old_working_set[i] = working_set[i];
		
		++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[] x;

	delete[] working_set;
	delete[] old_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;
		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];
	old_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);
		}

		if (counter == min(l*2/qpsize, 2000/qpsize))
		{
			bool same = true;
			for (i=0;i<qpsize;i++)
				if (old_working_set[i] != working_set[i]) 
				{
					same = false;
					break;
				}

			if (same)
				break;
		}

		for (i=0;i<qpsize;i++)
			old_working_set[i] = working_set[i];
		
		++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++)