svm.java

SVM是一种常用的模式分类机器学习算法

				{
					if(alpha[i] < 0)
					{
						alpha[i] = 0;
						alpha[j] = sum;
					}
				}
			}

			// update G

			double delta_alpha_i = alpha[i] - old_alpha_i;
			double delta_alpha_j = alpha[j] - old_alpha_j;

			for(int k=0;k<active_size;k++)
			{
				G[k] += Q_i[k]*delta_alpha_i + Q_j[k]*delta_alpha_j;
			}

			// update alpha_status and G_bar

			{
				boolean ui = is_upper_bound(i);
				boolean uj = is_upper_bound(j);
				update_alpha_status(i);
				update_alpha_status(j);
				int k;
				if(ui != is_upper_bound(i))
				{
					Q_i = Q.get_Q(i,l);
					if(ui)
						for(k=0;k<l;k++)
							G_bar[k] -= C_i * Q_i[k];
					else
						for(k=0;k<l;k++)
							G_bar[k] += C_i * Q_i[k];
				}

				if(uj != is_upper_bound(j))
				{
					Q_j = Q.get_Q(j,l);
					if(uj)
						for(k=0;k<l;k++)
							G_bar[k] -= C_j * Q_j[k];
					else
						for(k=0;k<l;k++)
							G_bar[k] += C_j * Q_j[k];
				}
			}

		}

		// calculate rho

		si.rho = calculate_rho();

		// 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;
		}

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

		si.upper_bound_p = Cp;
		si.upper_bound_n = Cn;

		System.out.print("\noptimization finished, #iter = "+iter+"\n");
	}

	// return 1 if already optimal, return 0 otherwise
	int select_working_set(int[] working_set)
	{
		// return i,j such that
		// i: maximizes -y_i * grad(f)_i, i in I_up(\alpha)
		// j: mimimizes the decrease of obj value
		//    (if quadratic coefficeint <= 0, replace it with tau)
		//    -y_j*grad(f)_j < -y_i*grad(f)_i, j in I_low(\alpha)
		
		double Gmax = -INF;
		double Gmax2 = -INF;
		int Gmax_idx = -1;
		int Gmin_idx = -1;
		double obj_diff_min = INF;
	
		for(int t=0;t<active_size;t++)
			if(y[t]==+1)	
			{
				if(!is_upper_bound(t))
					if(-G[t] >= Gmax)
					{
						Gmax = -G[t];
						Gmax_idx = t;
					}
			}
			else
			{
				if(!is_lower_bound(t))
					if(G[t] >= Gmax)
					{
						Gmax = G[t];
						Gmax_idx = t;
					}
			}
	
		int i = Gmax_idx;
		float[] Q_i = null;
		if(i != -1) // null Q_i not accessed: Gmax=-INF if i=-1
			Q_i = Q.get_Q(i,active_size);
	
		for(int j=0;j<active_size;j++)
		{
			if(y[j]==+1)
			{
				if (!is_lower_bound(j))
				{
					double grad_diff=Gmax+G[j];
					if (G[j] >= Gmax2)
						Gmax2 = G[j];
					if (grad_diff > 0)
					{
						double obj_diff; 
						double quad_coef=Q_i[i]+QD[j]-2*y[i]*Q_i[j];
						if (quad_coef > 0)
							obj_diff = -(grad_diff*grad_diff)/quad_coef;
						else
							obj_diff = -(grad_diff*grad_diff)/1e-12;
	
						if (obj_diff <= obj_diff_min)
						{
							Gmin_idx=j;
							obj_diff_min = obj_diff;
						}
					}
				}
			}
			else
			{
				if (!is_upper_bound(j))
				{
					double grad_diff= Gmax-G[j];
					if (-G[j] >= Gmax2)
						Gmax2 = -G[j];
					if (grad_diff > 0)
					{
						double obj_diff; 
						double quad_coef=Q_i[i]+QD[j]+2*y[i]*Q_i[j];
						if (quad_coef > 0)
							obj_diff = -(grad_diff*grad_diff)/quad_coef;
						else
							obj_diff = -(grad_diff*grad_diff)/1e-12;
	
						if (obj_diff <= obj_diff_min)
						{
							Gmin_idx=j;
							obj_diff_min = obj_diff;
						}
					}
				}
			}
		}

		if(Gmax+Gmax2 < eps)
			return 1;

		working_set[0] = Gmax_idx;
		working_set[1] = Gmin_idx;
		return 0;
	}

	// return 1 if already optimal, return 0 otherwise
	int max_violating_pair(int[] working_set)
	{
		// return i,j which maximize -grad(f)^T d , under constraint
		// if alpha_i == C, d != +1
		// if alpha_i == 0, d != -1

		double Gmax1 = -INF;		// max { -y_i * grad(f)_i | i in I_up(\alpha) }
		int Gmax1_idx = -1;

		int Gmax2_idx = -1;
		double Gmax2 = -INF;		// max { y_i * grad(f)_i | i in I_low(\alpha) }

		for(int i=0;i<active_size;i++)
		{
			if(y[i]==+1)	// y = +1
			{
				if(!is_upper_bound(i))	// d = +1
				{
					if(-G[i] >= Gmax1)
					{
						Gmax1 = -G[i];
						Gmax1_idx = i;
					}
				}
				if(!is_lower_bound(i))	// d = -1
				{
					if(G[i] >= Gmax2)
					{
						Gmax2 = G[i];
						Gmax2_idx = i;
					}
				}
			}
			else		// y = -1
			{
				if(!is_upper_bound(i))	// d = +1
				{
					if(-G[i] >= Gmax2)
					{
						Gmax2 = -G[i];
						Gmax2_idx = i;
					}
				}
				if(!is_lower_bound(i))	// d = -1
				{
					if(G[i] >= Gmax1)
					{
						Gmax1 = G[i];
						Gmax1_idx = i;
					}
				}
			}
		}

		if(Gmax1+Gmax2 < eps)
	 		return 1;

		working_set[0] = Gmax1_idx;
		working_set[1] = Gmax2_idx;
		return 0;
	}

	void do_shrinking()
	{
		int i,j,k;
		int[] working_set = new int[2];
		if(max_violating_pair(working_set)!=0) return;
		i = working_set[0];
		j = working_set[1];
		double Gm1 = -y[j]*G[j];
		double Gm2 = y[i]*G[i];

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

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

		// unshrink, check all variables again before final iterations

		if(unshrinked || -(Gm1 + Gm2) > eps*10) return;

		unshrinked = true;
		reconstruct_gradient();

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

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

	double calculate_rho()
	{
		double r;
		int nr_free = 0;
		double ub = INF, lb = -INF, sum_free = 0;
		for(int i=0;i<active_size;i++)
		{
			double yG = y[i]*G[i];

			if(is_lower_bound(i))
			{
				if(y[i] > 0)
					ub = Math.min(ub,yG);
				else
					lb = Math.max(lb,yG);
			}
			else if(is_upper_bound(i))
			{
				if(y[i] < 0)
					ub = Math.min(ub,yG);
				else
					lb = Math.max(lb,yG);
			}
			else
			{
				++nr_free;
				sum_free += yG;
			}
		}

		if(nr_free>0)
			r = sum_free/nr_free;
		else
			r = (ub+lb)/2;

		return r;
	}

}

//
// Solver for nu-svm classification and regression
//
// additional constraint: e^T \alpha = constant
//
final class Solver_NU extends Solver
{
	private SolutionInfo si;

	void Solve(int l, QMatrix Q, double[] b, byte[] y,
		   double[] alpha, double Cp, double Cn, double eps,
		   SolutionInfo si, int shrinking)
	{
		this.si = si;
		super.Solve(l,Q,b,y,alpha,Cp,Cn,eps,si,shrinking);
	}

	// return 1 if already optimal, return 0 otherwise
	int select_working_set(int[] working_set)
	{
		// return i,j such that y_i = y_j and
		// i: maximizes -y_i * grad(f)_i, i in I_up(\alpha)
		// j: minimizes the decrease of obj value
		//    (if quadratic coefficeint <= 0, replace it with tau)
		//    -y_j*grad(f)_j < -y_i*grad(f)_i, j in I_low(\alpha)
	
		double Gmaxp = -INF;
		double Gmaxp2 = -INF;
		int Gmaxp_idx = -1;
	
		double Gmaxn = -INF;
		double Gmaxn2 = -INF;
		int Gmaxn_idx = -1;
	
		int Gmin_idx = -1;
		double obj_diff_min = INF;
	
		for(int t=0;t<active_size;t++)
			if(y[t]==+1)
			{
				if(!is_upper_bound(t))
					if(-G[t] >= Gmaxp)
					{
						Gmaxp = -G[t];
						Gmaxp_idx = t;
					}
			}
			else
			{
				if(!is_lower_bound(t))
					if(G[t] >= Gmaxn)
					{
						Gmaxn = G[t];
						Gmaxn_idx = t;
					}
			}
	
		int ip = Gmaxp_idx;
		int in = Gmaxn_idx;
		float[] Q_ip = null;
		float[] Q_in = null;
		if(ip != -1) // null Q_ip not accessed: Gmaxp=-INF if ip=-1
			Q_ip = Q.get_Q(ip,active_size);
		if(in != -1)
			Q_in = Q.get_Q(in,active_size);
	
		for(int j=0;j<active_size;j++)
		{
			if(y[j]==+1)
			{
				if (!is_lower_bound(j))	
				{
					double grad_diff=Gmaxp+G[j];
					if (G[j] >= Gmaxp2)
						Gmaxp2 = G[j];
					if (grad_diff > 0)
					{
						double obj_diff; 
						double quad_coef = Q_ip[ip]+QD[j]-2*Q_ip[j];
						if (quad_coef > 0)
							obj_diff = -(grad_diff*grad_diff)/quad_coef;
						else
							obj_diff = -(grad_diff*grad_diff)/1e-12;
	
						if (obj_diff <= obj_diff_min)
						{
							Gmin_idx=j;
							obj_diff_min = obj_diff;
						}
					}
				}
			}
			else
			{
				if (!is_upper_bound(j))
				{
					double grad_diff=Gmaxn-G[j];
					if (-G[j] >= Gmaxn2)
						Gmaxn2 = -G[j];
					if (grad_diff > 0)
					{
						double obj_diff; 
						double quad_coef = Q_in[in]+QD[j]-2*Q_in[j];
						if (quad_coef > 0)
							obj_diff = -(grad_diff*grad_diff)/quad_coef;
						else
							obj_diff = -(grad_diff*grad_diff)/1e-12;
	
						if (obj_diff <= obj_diff_min)
						{
							Gmin_idx=j;
							obj_diff_min = obj_diff;
						}
					}
				}
			}
		}

		if(Math.max(Gmaxp+Gmaxp2,Gmaxn+Gmaxn2) < eps)
 			return 1;
	
		if(y[Gmin_idx] == +1)
			working_set[0] = Gmaxp_idx;
		else
			working_set[0] = Gmaxn_idx;
		working_set[1] = Gmin_idx;
	
		return 0;
	}

	void do_shrinking()
	{
		double Gmax1 = -INF;	// max { -y_i * grad(f)_i | y_i = +1, i in I_up(\alpha) }
		double Gmax2 = -INF;	// max { y_i * grad(f)_i | y_i = +1, i in I_low(\alpha) }
		double Gmax3 = -INF;	// max { -y_i * grad(f)_i | y_i = -1, i in I_up(\alpha) }
		double Gmax4 = -INF;	// max { y_i * grad(f)_i | y_i = -1, i in I_low(\alpha) }
 
		// find maximal violating pair first
		int k;
		for(k=0;k<active_size;k++)
		{
			if(!is_upper_bound(k))
			{
				if(y[k]==+1)
				{
					if(-G[k] > Gmax1) Gmax1 = -G[k];
				}
				else	if(-G[k] > Gmax3) Gmax3 = -G[k];
			}
			if(!is_lower_bound(k))
			{
				if(y[k]==+1)
				{	
					if(G[k] > Gmax2) Gmax2 = G[k];
				}
				else	if(G[k] > Gmax4) Gmax4 = G[k];
			}
		}

		// shrinking

		double Gm1 = -Gmax2;
		double Gm2 = -Gmax1;
		double Gm3 = -Gmax4;
		double Gm4 = -Gmax3;

		for(k=0;k<active_size;k++)
		{
			if(is_lower_bound(k))
			{
				if(y[k]==+1)
				{
					if(-G[k] >= Gm1) continue;
				}
				else	if(-G[k] >= Gm3) continue;
			}
			else if(is_upper_bound(k))
			{
				if(y[k]==+1)
				{
					if(G[k] >= Gm2) continue;
				}
				else	if(G[k] >= Gm4) continue;
			}
			else continue;

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

		// unshrink, check all variables again before final iterations

		if(unshrinked || Math.max(-(Gm1+Gm2),-(Gm3+Gm4)) > eps*10) return;
	
		unshrinked = true;
		reconstruct_gradient();

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