?? solver.cs
字號:
private float[][] buffer;
private float[] QD;
public SVR_Q(Problem prob, Parameter param)
: base(prob.Count, prob.X, param)
{
l = prob.Count;
cache = new Cache(l, (long)(param.CacheSize * (1 << 20)));
QD = new float[2 * l];
sign = new short[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;
QD[k] = (float)kernel_function(k, k);
QD[k + l] = QD[k];
}
buffer = new float[2][];
buffer[0] = new float[2 * l];
buffer[1] = new float[2 * l];
next_buffer = 0;
}
public override void swap_index(int i, int j)
{
do { short _ = sign[i]; sign[i] = sign[j]; sign[j] = _; } while (false);
do { int _ = index[i]; index[i] = index[j]; index[j] = _; } while (false);
do { float _ = QD[i]; QD[i] = QD[j]; QD[j] = _; } while (false);
}
public override float[] get_Q(int i, int len)
{
float[][] data = new float[1][];
int real_i = index[i];
if (cache.get_data(real_i, data, l) < l)
{
for (int j = 0; j < l; j++)
data[0][j] = (float)kernel_function(real_i, j);
}
// reorder and copy
float[] buf = buffer[next_buffer];
next_buffer = 1 - next_buffer;
short si = sign[i];
for (int j = 0; j < len; j++)
buf[j] = si * sign[j] * data[0][index[j]];
return buf;
}
public override float[] get_QD()
{
return QD;
}
}
internal static class Procedures
{
//
// construct and solve various formulations
//
private static void solve_c_svc(Problem prob, Parameter param,
double[] alpha, Solver.SolutionInfo si,
double Cp, double Cn)
{
int l = prob.Count;
double[] minus_ones = new double[l];
short[] y = new short[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;
}
Solver s = new Solver();
s.Solve(l, new SVC_Q(prob, param, y), minus_ones, y,
alpha, Cp, Cn, param.EPS, si, param.Shrinking);
double sum_alpha = 0;
for (i = 0; i < l; i++)
sum_alpha += alpha[i];
if (Cp == Cn)
Debug.Write("nu = " + sum_alpha / (Cp * prob.Count) + "\n");
for (i = 0; i < l; i++)
alpha[i] *= y[i];
}
private static void solve_nu_svc(Problem prob, Parameter param,
double[] alpha, Solver.SolutionInfo si)
{
int i;
int l = prob.Count;
double nu = param.Nu;
short[] y = new short[l];
for (i = 0; i < l; i++)
if (prob.Y[i] > 0)
y[i] = +1;
else
y[i] = -1;
double sum_pos = nu * l / 2;
double sum_neg = nu * l / 2;
for (i = 0; i < l; i++)
if (y[i] == +1)
{
alpha[i] = Math.Min(1.0, sum_pos);
sum_pos -= alpha[i];
}
else
{
alpha[i] = Math.Min(1.0, sum_neg);
sum_neg -= alpha[i];
}
double[] zeros = new double[l];
for (i = 0; i < l; i++)
zeros[i] = 0;
Solver_NU s = new Solver_NU();
s.Solve(l, new SVC_Q(prob, param, y), zeros, y,
alpha, 1.0, 1.0, param.EPS, si, param.Shrinking);
double r = si.r;
Debug.Write("C = " + 1 / r + "\n");
for (i = 0; i < l; i++)
alpha[i] *= y[i] / r;
si.rho /= r;
si.obj /= (r * r);
si.upper_bound_p = 1 / r;
si.upper_bound_n = 1 / r;
}
private static void solve_one_class(Problem prob, Parameter param,
double[] alpha, Solver.SolutionInfo si)
{
int l = prob.Count;
double[] zeros = new double[l];
short[] ones = new short[l];
int i;
int n = (int)(param.Nu * prob.Count); // # of alpha's at upper bound
for (i = 0; i < n; i++)
alpha[i] = 1;
if (n < prob.Count)
alpha[n] = param.Nu * prob.Count - n;
for (i = n + 1; i < l; i++)
alpha[i] = 0;
for (i = 0; i < l; i++)
{
zeros[i] = 0;
ones[i] = 1;
}
Solver s = new Solver();
s.Solve(l, new ONE_CLASS_Q(prob, param), zeros, ones,
alpha, 1.0, 1.0, param.EPS, si, param.Shrinking);
}
private static void solve_epsilon_svr(Problem prob, Parameter param,
double[] alpha, Solver.SolutionInfo si)
{
int l = prob.Count;
double[] alpha2 = new double[2 * l];
double[] linear_term = new double[2 * l];
short[] y = new short[2 * l];
int i;
for (i = 0; i < l; i++)
{
alpha2[i] = 0;
linear_term[i] = param.P - prob.Y[i];
y[i] = 1;
alpha2[i + l] = 0;
linear_term[i + l] = param.P + prob.Y[i];
y[i + l] = -1;
}
Solver s = new Solver();
s.Solve(2 * l, new SVR_Q(prob, param), linear_term, y,
alpha2, param.C, param.C, param.EPS, si, param.Shrinking);
double sum_alpha = 0;
for (i = 0; i < l; i++)
{
alpha[i] = alpha2[i] - alpha2[i + l];
sum_alpha += Math.Abs(alpha[i]);
}
Debug.Write("nu = " + sum_alpha / (param.C * l) + "\n");
}
private static void solve_nu_svr(Problem prob, Parameter param,
double[] alpha, Solver.SolutionInfo si)
{
int l = prob.Count;
double C = param.C;
double[] alpha2 = new double[2 * l];
double[] linear_term = new double[2 * l];
short[] y = new short[2 * l];
int i;
double sum = C * param.Nu * l / 2;
for (i = 0; i < l; i++)
{
alpha2[i] = alpha2[i + l] = Math.Min(sum, C);
sum -= alpha2[i];
linear_term[i] = -prob.Y[i];
y[i] = 1;
linear_term[i + l] = prob.Y[i];
y[i + l] = -1;
}
Solver_NU s = new Solver_NU();
s.Solve(2 * l, new SVR_Q(prob, param), linear_term, y, alpha2, C, C, param.EPS, si, param.Shrinking);
Debug.Write("epsilon = " + (-si.r) + "\n");
for (i = 0; i < l; i++)
alpha[i] = alpha2[i] - alpha2[i + l];
}
//
// decision_function
//
private class decision_function
{
public double[] alpha;
public double rho;
};
static decision_function svm_train_one(
Problem prob, Parameter param,
double Cp, double Cn)
{
double[] alpha = new double[prob.Count];
Solver.SolutionInfo si = new Solver.SolutionInfo();
switch (param.SvmType)
{
case SvmType.C_SVC:
solve_c_svc(prob, param, alpha, si, Cp, Cn);
break;
case SvmType.NU_SVC:
solve_nu_svc(prob, param, alpha, si);
break;
case SvmType.ONE_CLASS:
solve_one_class(prob, param, alpha, si);
break;
case SvmType.EPSILON_SVR:
solve_epsilon_svr(prob, param, alpha, si);
break;
case SvmType.NU_SVR:
solve_nu_svr(prob, param, alpha, si);
break;
}
Debug.Write("obj = " + si.obj + ", rho = " + si.rho + "\n");
// output SVs
int nSV = 0;
int nBSV = 0;
for (int i = 0; i < prob.Count; i++)
{
if (Math.Abs(alpha[i]) > 0)
{
++nSV;
if (prob.Y[i] > 0)
{
if (Math.Abs(alpha[i]) >= si.upper_bound_p)
++nBSV;
}
else
{
if (Math.Abs(alpha[i]) >= si.upper_bound_n)
++nBSV;
}
}
}
Debug.Write("nSV = " + nSV + ", nBSV = " + nBSV + "\n");
decision_function f = new decision_function();
f.alpha = alpha;
f.rho = si.rho;
return f;
}
// Platt's binary SVM Probablistic Output: an improvement from Lin et al.
private static void sigmoid_train(int l, double[] dec_values, double[] labels,
double[] probAB)
{
double A, B;
double prior1 = 0, prior0 = 0;
int i;
for (i = 0; i < l; i++)
if (labels[i] > 0) prior1 += 1;
else prior0 += 1;
int max_iter = 100; // Maximal number of iterations
double min_step = 1e-10; // Minimal step taken in line search
double sigma = 1e-3; // For numerically strict PD of Hessian
double eps = 1e-5;
double hiTarget = (prior1 + 1.0) / (prior1 + 2.0);
double loTarget = 1 / (prior0 + 2.0);
double[] t = new double[l];
double fApB, p, q, h11, h22, h21, g1, g2, det, dA, dB, gd, stepsize;
double newA, newB, newf, d1, d2;
int iter;
// Initial Point and Initial Fun Value
A = 0.0; B = Math.Log((prior0 + 1.0) / (prior1 + 1.0));
double fval = 0.0;
for (i = 0; i < l; i++)
{
if (labels[i] > 0) t[i] = hiTarget;
else t[i] = loTarget;
fApB = dec_values[i] * A + B;
if (fApB >= 0)
fval += t[i] * fApB + Math.Log(1 + Math.Exp(-fApB));
else
fval += (t[i] - 1) * fApB + Math.Log(1 + Math.Exp(fApB));
}
for (iter = 0; iter < max_iter; iter++)
{
// Update Gradient and Hessian (use H' = H + sigma I)
h11 = sigma; // numerically ensures strict PD
h22 = sigma;
h21 = 0.0; g1 = 0.0; g2 = 0.0;
for (i = 0; i < l; i++)
{
fApB = dec_values[i] * A + B;
if (fApB >= 0)
{
p = Math.Exp(-fApB) / (1.0 + Math.Exp(-fApB));
q = 1.0 / (1.0 + Math.Exp(-fApB));
}
else
{
p = 1.0 / (1.0 + Math.Exp(fApB));
q = Math.Exp(fApB) / (1.0 + Math.Exp(fApB));
}
d2 = p * q;
h11 += dec_values[i] * dec_values[i] * d2;
h22 += d2;
h21 += dec_values[i] * d2;
d1 = t[i] - p;
g1 += dec_values[i] * d1;
g2 += d1;
}
// Stopping Criteria
if (Math.Abs(g1) < eps && Math.Abs(g2) < eps)
break;
// Finding Newton direction: -inv(H') * g
det = h11 * h22 - h21 * h21;
dA = -(h22 * g1 - h21 * g2) / det;
dB = -(-h21 * g1 + h11 * g2) / det;
gd = g1 * dA + g2 * dB;
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -