bsvm.cpp

svm的实现源码

			if (G_i[m] >= th)
				goto out1;
		for (m++;m<nr_class;m++)
			if (G_i[m] >= th)
				goto out1;
		
		swap_index(i, active_size);
		++active_size;
		++i;
	out1:	;
	}
}

int compar(const void *a, const void *b)
{
	if (*(double *)a > *(double *)b)
		return -1;
	else
		if (*(double *)a < *(double *)b)
			return 1;
	return 0;
}
void Solver_SPOC::solve_sub_problem(double A, double *B, double C, double *nu)
{
	int r;
	double *D;
	
	clone(D, B, nr_class+1);
	qsort(D, nr_class, sizeof(double), compar);
	D[nr_class] = -INF;
	
	double phi = D[0] - A*C;
	for (r=0;phi<(r+1)*D[r+1];r++)
		phi += D[r+1];
	delete[] D;
		
	phi /= (r+1);
	for (r=0;r<nr_class;r++)
		nu[r] = min((double) 0, phi - B[r])/A;
}

#ifdef __cplusplus
extern "C" {
#endif
void solvebqp(struct BQP*);
#ifdef __cplusplus
}
#endif

class Solver_B {
public:
	Solver_B() {};
	virtual ~Solver_B() {};

	struct SolutionInfo {
		double obj;
		double *upper_bound;
	};

	virtual void Solve(int l, const Kernel& Q, double *b_, schar *y_,
	double *alpha_, double Cp, double Cn, double eps, SolutionInfo* si, 
	int shrinking, int qpsize);
protected:
	int active_size;
	double *G;		// gradient of objective function
	enum { LOWER_BOUND, UPPER_BOUND, FREE };
	char *alpha_status;	// LOWER_BOUND, UPPER_BOUND, FREE
	double *alpha;
	const Kernel *Q;
	double eps;

	int *active_set;
	double *G_bar;		// gradient, if we treat free variables as 0
	int l;
	bool unshrinked;	// XXX

	int qpsize;
	int *working_set;
	int *old_working_set;

	virtual double get_C(int i)
	{
		return (y[i] > 0)? Cp : Cn;
	}
	void update_alpha_status(int i)
	{
		if(alpha[i] >= get_C(i))
			alpha_status[i] = UPPER_BOUND;
		else if(alpha[i] <= 0)
			alpha_status[i] = LOWER_BOUND;
		else alpha_status[i] = FREE;
	}
	bool is_upper_bound(int i) { return alpha_status[i] == UPPER_BOUND; }
	bool is_lower_bound(int i) { return alpha_status[i] == LOWER_BOUND; }
	bool is_free(int i) { return alpha_status[i] == FREE; }
	virtual void swap_index(int i, int j);
	virtual void reconstruct_gradient();
	virtual void shrink_one(int k);
	virtual void unshrink_one(int k);
	double select_working_set(int &q);
	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];
	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] = 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;

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