svm.cs

支持向量机程序,非常有用,可以供大家实验使用,改进.希望能多大家工作有帮助.

		
		//UPGRADE_NOTE: 方法“Solve”的访问修饰符被更改为“internal”。 'ms-help://MS.VSCC.2003/commoner/redir/redirect.htm?keyword="jlca1204"'
		internal virtual void  Solve(int l, QMatrix Q, double[] b_, sbyte[] y_, double[] alpha_, double Cp, double Cn, double eps, SolutionInfo si, int shrinking)
		{
			this.l = l;
			this.Q = Q;
			QD = Q.get_QD();
			b = new double[b_.Length];
			b_.CopyTo(b, 0);
			y = new sbyte[y_.Length];
			y_.CopyTo(y, 0);
			alpha = new double[alpha_.Length];
			alpha_.CopyTo(alpha, 0);
			this.Cp = Cp;
			this.Cn = Cn;
			this.eps = eps;
			this.unshrinked = false;
			
			// initialize alpha_status
			{
				alpha_status = new sbyte[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;
			}
			
			// 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))
					{
						float[] Q_i = Q.get_Q(i, l);
						double alpha_i = alpha[i];
						int j;
						for (j = 0; j < l; j++)
							G[j] += alpha_i * Q_i[j];
						if (is_upper_bound(i))
							for (j = 0; j < l; j++)
								G_bar[j] += get_C(i) * Q_i[j];
					}
			}
			
			// optimization step
			
			int iter = 0;
			int counter = System.Math.Min(l, 1000) + 1;
			int[] working_set = new int[2];
			
			while (true)
			{
				// show progress and do shrinking
				
				if (--counter == 0)
				{
					counter = System.Math.Min(l, 1000);
					if (shrinking != 0)
						do_shrinking();
					System.Console.Error.Write(".");
				}
				
				if (select_working_set(working_set) != 0)
				{
					// reconstruct the whole gradient
					reconstruct_gradient();
					// reset active set size and check
					active_size = l;
					System.Console.Error.Write("*");
					if (select_working_set(working_set) != 0)
						break;
					else
						counter = 1; // do shrinking next iteration
				}
				
				int i = working_set[0];
				int j = working_set[1];
				
				++iter;
				
				// update alpha[i] and alpha[j], handle bounds carefully
				
				float[] Q_i = Q.get_Q(i, active_size);
				float[] Q_j = Q.get_Q(j, active_size);
				
				double C_i = get_C(i);
				double C_j = get_C(j);
				
				double old_alpha_i = alpha[i];
				double old_alpha_j = alpha[j];
				
				if (y[i] != y[j])
				{
					double quad_coef = Q_i[i] + Q_j[j] + 2 * Q_i[j];
					if (quad_coef <= 0)
						quad_coef = 1e-12;
					double delta = (- G[i] - G[j]) / quad_coef;
					double diff = alpha[i] - alpha[j];
					alpha[i] += delta;
					alpha[j] += delta;
					
					if (diff > 0)
					{
						if (alpha[j] < 0)
						{
							alpha[j] = 0;
							alpha[i] = diff;
						}
					}
					else
					{
						if (alpha[i] < 0)
						{
							alpha[i] = 0;
							alpha[j] = - diff;
						}
					}
					if (diff > C_i - C_j)
					{
						if (alpha[i] > C_i)
						{
							alpha[i] = C_i;
							alpha[j] = C_i - diff;
						}
					}
					else
					{
						if (alpha[j] > C_j)
						{
							alpha[j] = C_j;
							alpha[i] = C_j + diff;
						}
					}
				}
				else
				{
					double quad_coef = Q_i[i] + Q_j[j] - 2 * Q_i[j];
					if (quad_coef <= 0)
						quad_coef = 1e-12;
					double delta = (G[i] - G[j]) / quad_coef;
					double sum = alpha[i] + alpha[j];
					alpha[i] -= delta;
					alpha[j] += delta;
					
					if (sum > C_i)
					{
						if (alpha[i] > C_i)
						{
							alpha[i] = C_i;
							alpha[j] = sum - C_i;
						}
					}
					else
					{
						if (alpha[j] < 0)
						{
							alpha[j] = 0;
							alpha[i] = sum;
						}
					}
					if (sum > C_j)
					{
						if (alpha[j] > C_j)
						{
							alpha[j] = C_j;
							alpha[i] = sum - C_j;
						}
					}
					else
					{
						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
				
				{
					bool ui = is_upper_bound(i);
					bool 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.Console.Out.Write("\noptimization finished, #iter = " + iter + "\r\n");
		}
		
		// return 1 if already optimal, return 0 otherwise
		//UPGRADE_NOTE: 方法“select_working_set”的访问修饰符被更改为“internal”。 'ms-help://MS.VSCC.2003/commoner/redir/redirect.htm?keyword="jlca1204"'
		internal virtual 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;
			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 (grad_diff >= eps)
						{
							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 (grad_diff >= eps)
						{
							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 (Gmin_idx == - 1)
				return 1;
			
			working_set[0] = Gmax_idx;
			working_set[1] = Gmin_idx;
			return 0;
		}
		
		// return 1 if already optimal, return 0 otherwise
		internal virtual 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;
						}
					}
				}
				// y = -1
				else
				{
					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;
		}
		
		//UPGRADE_NOTE: 方法“do_shrinking”的访问修饰符被更改为“internal”。 'ms-help://MS.VSCC.2003/commoner/redir/redirect.htm?keyword="jlca1204"'
		internal virtual 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))