37 #ifndef SHARK_ALGORITHMS_QP_QPMCDECOMP_H 38 #define SHARK_ALGORITHMS_QP_QPMCDECOMP_H 49 #define ITERATIONS_BETWEEN_SHRINKING 1000 52 #define GRADIENT_UPDATE(r, from, to, mu, q) \ 54 std::size_t a, b, p; \ 55 for (a=(from); a<(to); a++) \ 58 tExample& ex = example[a]; \ 59 typename QpSparseArray<QpFloatType>::Row const& row = M.row(classes * (r) + ex.y); \ 60 QpFloatType def = row.defaultvalue; \ 63 for (b=0; b<row.size; b++) \ 65 p = row.entry[b].index; \ 66 gradient(ex.var[p]) -= (mu) * row.entry[b].value * k; \ 71 for (b=0; b<row.size; b++) \ 73 p = row.entry[b].index; \ 74 gradient(ex.var[p]) -= (mu) * (row.entry[b].value - def) * k; \ 76 double upd = (mu) * def * k; \ 77 for (b=0; b<ex.active; b++) gradient(ex.avar[b]) -= upd; \ 82 #define GAIN_SELECTION_BOX(qif) \ 84 if (f >= activeVar) continue; \ 85 double gf = gradient(f); \ 86 double af = alpha(f); \ 87 if ((af > 0.0 && gf < 0.0) || (af < C && gf > 0.0)) \ 89 double df = variable[f].diagonal; \ 90 double diag_q = di * df; \ 91 double det_q = diag_q - qif * qif; \ 92 if (det_q < 1e-12 * diag_q) \ 94 if (di == 0.0 && df == 0.0) \ 95 { if (f != i) { j = f; return ret; } } \ 98 if (di * gf - df * gi != 0.0) \ 99 { if (f != i) { j = f; return ret; } } \ 102 double g2 = gf*gf + gi*gi; \ 103 double gain = (g2*g2) / (gf*gf*df + 2.0*gf*gi*qif + gi*gi*di); \ 104 if (gain > bestgain) { if (f != i) { bestgain = gain; j = f; } } \ 110 double gain = (gf*gf*di - 2.0*gf*gi*qif + gi*gi*df) / det_q; \ 111 if (gain > bestgain) { if (f != i) { bestgain = gain; j = f; } } \ 115 #define GAIN_SELECTION_TRIANGLE(qif) \ 117 if (f >= activeVar) continue; \ 118 double gf = gradient(f); \ 119 double af = alpha(f); \ 120 if ((af > 0.0 && gf < 0.0) || (varsum < C && gf > 0.0)) \ 122 double df = variable[f].diagonal; \ 123 double diag_q = di * df; \ 124 double det_q = diag_q - qif * qif; \ 125 if (det_q < 1e-12 * diag_q) \ 127 if (di == 0.0 && df == 0.0) \ 128 { if (f != i) { j = f; return ret; } } \ 131 if (di * gf - df * gi != 0.0) \ 132 { if (f != i) { j = f; return ret; } } \ 135 double g2 = gf*gf + gi*gi; \ 136 double gain = (g2*g2) / (gf*gf*df + 2.0*gf*gi*qif + gi*gi*di); \ 137 if (gain > bestgain) { if (f != i) { bestgain = gain; j = f; } } \ 143 double gain = (gf*gf*di - 2.0*gf*gi*qif + gi*gi*df) / det_q; \ 144 if (gain > bestgain) { if (f != i) { bestgain = gain; j = f; } } \ 186 template <
class Matrix>
208 RealMatrix
const& _gamma,
209 UIntVector
const& _rho,
212 bool sumToZeroConstraint)
227 for (p=0; p<cardP; p++) if (rho[p] >=
cardR)
cardR =
rho[p] +1;
237 "[QpMcDecomp::QpMcDecomp] dimension conflict" 243 throw SHARKEXCEPTION(
"[QpMcDecomp::QpMcDecomp] currently this solver supports only bijective or trivial rho");
260 RealVector& solutionAlpha,
263 RealVector* solutionBias = NULL)
268 std::size_t v,
w, i, e;
274 alpha = solutionAlpha;
292 unsigned int y = target.
element(i);
300 for (p=0; p<
cardP; p++, v++)
305 double Q =
M(
classes * (y * cardP + p) + y, p) * k;
309 example[i].varsum += solutionAlpha(v);
310 double lin =
gamma(y, p);
322 e = (std::size_t)(-1);
323 QpFloatType* q = NULL;
327 for (p=0; p<
cardP; p++, v++)
329 double av =
alpha(v);
338 unsigned int r = y*cardP+p;
353 unsigned long long iter = 0;
358 if (acc < dualAccuracy)
419 unsigned int r =
cardP * y + p;
437 double& varsum =
example[i].varsum;
438 double max_mu = C - varsum;
439 double max_alpha = max_mu +
alpha(v);
443 alpha(v) = max_alpha;
454 if (mu >= C -
alpha(v))
469 unsigned int yv =
example[iv].y;
473 unsigned int yw =
example[iw].y;
478 unsigned int rv =
cardP*yv+pv;
479 unsigned int rw =
cardP*yw+pw;
484 double Qvw =
M(
classes * rv + yw, pw) * qv[iw];
502 double& varsum1 =
example[iv].varsum;
503 double& varsum2 =
example[iw].varsum;
504 double U1 = C - (varsum1 -
alpha(v));
505 double U2 = C - (varsum2 -
alpha(w));
513 if (
alpha(v) == U1) varsum1 =
C;
514 else varsum1 += mu_v;
515 if (
alpha(w) == U2) varsum2 =
C;
516 else varsum2 += mu_w;
534 if (checkCounter == 0)
544 if (current_time - start_time > stop.
maxSeconds)
546 if (prop != NULL) prop->type =
QpTimeout;
562 double objective = 0.0;
566 solutionAlpha(w) =
alpha(v);
577 prop->value = objective;
578 prop->iterations = iter;
579 prop->seconds = finish_time - start_time;
607 RealVector& solutionAlpha,
610 RealVector* solutionBias = NULL)
615 std::size_t v,
w, i, e;
621 alpha = solutionAlpha;
639 unsigned int y = target.
element(i);
646 for (p=0; p<
cardP; p++, v++)
655 example[i].varsum += solutionAlpha(v);
656 double lin =
gamma(y, p);
668 QpFloatType* q = NULL;
672 for (p=0; p<
cardP; p++, v++)
674 double av =
alpha(v);
683 unsigned int r = y*cardP+p;
698 unsigned long long iter = 0;
703 if (acc < dualAccuracy)
759 unsigned int yv =
example[iv].y;
768 unsigned int rv =
cardP*yv+pv;
782 throw SHARKEXCEPTION(
"[QpMcDecomp::solveSMO] SMO is implemented only for box constraints");
786 if (v != w)
throw SHARKEXCEPTION(
"[QpMcDecomp::solveSMO] internal error");
804 double a =
alpha(v) + mu_v;
828 if (checkCounter == 0)
838 if (current_time - start_time > stop.
maxSeconds)
840 if (prop != NULL) prop->type =
QpTimeout;
855 double objective = 0.0;
859 solutionAlpha(w) =
alpha(v);
869 prop->value = objective;
870 prop->iterations = iter;
871 prop->seconds = finish_time - start_time;
955 SHARK_CHECK(
bias != NULL,
"[QpMcDecomp::solveForBias] internal error");
959 RealVector stepsize(
classes,epsilon);
974 unsigned int largest_p =
cardP;
975 double largest_value = 0.0;
976 for (p=0; p<
cardP; p++)
978 std::size_t v = ex.
var[p];
981 if (g > largest_value)
987 if (largest_p < cardP)
1020 if (g > 0.0) step(c) = -stepsize(c);
1021 else if (g < 0.0) step(c) = stepsize(c);
1023 double gg = prev(c) * grad(c);
1024 if (gg > 0.0) stepsize(c) *= 1.2;
1025 else stepsize(c) *= 0.5;
1045 for (b=0; b<row.
size; b++)
1054 double maxstep = 0.0;
1055 for (c=0; c<classes; c++) if (stepsize(c) > maxstep) maxstep = stepsize(c);
1056 if (maxstep < 0.01 * epsilon)
break;
1149 if (! boost::math::isnormal(mu))
1164 double a = alpha + mu;
1187 double gi,
double gj,
1188 double Qii,
double Qij,
double Qjj,
1189 double Ui,
double Uj,
1190 double& mui,
double& muj)
1193 double detQ = Qii * Qjj - Qij * Qij;
1194 mui = (Qjj * gi - Qij * gj) / detQ;
1195 muj = (Qii * gj - Qij * gi) / detQ;
1196 double opti = alphai + mui;
1197 double optj = alphaj + muj;
1198 if (boost::math::isnormal(opti) && boost::math::isnormal(optj) && opti > 0.0 && optj > 0.0 && opti < Ui && optj < Uj)
1206 struct tEdgeSolution
1213 tEdgeSolution solution[4];
1214 tEdgeSolution* sol = solution;
1215 tEdgeSolution* best = solution;
1216 double gain, bestgain = 0.0;
1220 sol->alphaj = alphaj;
1222 solveEdge(sol->alphaj, gj + Qij * alphai, Qjj, Uj, sol->muj);
1223 gain = sol->mui * (gi - 0.5 * (Qii*sol->mui + Qij*sol->muj))
1224 + sol->muj * (gj - 0.5 * (Qij*sol->mui + Qjj*sol->muj));
1225 if (gain > bestgain) { bestgain = gain; best = sol; }
1230 sol->alphai = alphai;
1233 solveEdge(sol->alphai, gi + Qij * alphaj, Qii, Ui, sol->mui);
1234 gain = sol->mui * (gi - 0.5 * (Qii*sol->mui + Qij*sol->muj))
1235 + sol->muj * (gj - 0.5 * (Qij*sol->mui + Qjj*sol->muj));
1236 if (gain > bestgain) { bestgain = gain; best = sol; }
1242 sol->alphaj = alphaj;
1243 sol->mui = Ui - alphai;
1244 solveEdge(sol->alphaj, gj - Qij * sol->mui, Qjj, Uj, sol->muj);
1245 gain = sol->mui * (gi - 0.5 * (Qii*sol->mui + Qij*sol->muj))
1246 + sol->muj * (gj - 0.5 * (Qij*sol->mui + Qjj*sol->muj));
1247 if (gain > bestgain) { bestgain = gain; best = sol; }
1252 sol->alphai = alphai;
1254 sol->muj = Uj - alphaj;
1255 solveEdge(sol->alphai, gi - Qij * sol->muj, Qii, Ui, sol->mui);
1256 gain = sol->mui * (gi - 0.5 * (Qii*sol->mui + Qij*sol->muj))
1257 + sol->muj * (gj - 0.5 * (Qij*sol->mui + Qjj*sol->muj));
1258 if (gain > bestgain) { bestgain = gain; best = sol; }
1261 alphai = best->alphai;
1262 alphaj = best->alphaj;
1273 double gi,
double gj,
1274 double Qii,
double Qij,
double Qjj,
1275 double& mui,
double& muj)
1278 double V =
C - alphasum;
1279 double U = V + alphai + alphaj;
1280 double detQ = Qii * Qjj - Qij * Qij;
1281 mui = (Qjj * gi - Qij * gj) / detQ;
1282 muj = (Qii * gj - Qij * gi) / detQ;
1283 double opti = alphai + mui;
1284 double optj = alphaj + muj;
1285 if (boost::math::isnormal(opti) && boost::math::isnormal(optj) && opti > 0.0 && optj > 0.0 && opti + optj < U)
1289 alphasum += mui + muj;
1290 if (alphasum >
C) alphasum =
C;
1295 struct tEdgeSolution
1303 tEdgeSolution solution[3];
1304 tEdgeSolution* sol = solution;
1305 tEdgeSolution* best = NULL;
1306 double gain, bestgain = -1e100;
1310 sol->alphaj = alphaj;
1312 solveEdge(sol->alphaj, gj + Qij * alphai, Qjj, V + alphaj, sol->muj);
1313 sol->alphasum = alphasum + sol->mui + sol->muj;
1314 gain = sol->mui * (gi - 0.5 * (Qii*sol->mui + Qij*sol->muj))
1315 + sol->muj * (gj - 0.5 * (Qij*sol->mui + Qjj*sol->muj));
1316 if (gain > bestgain) { bestgain = gain; best = sol; }
1321 sol->alphai = alphai;
1324 solveEdge(sol->alphai, gi + Qij * alphaj, Qii, V + alphai, sol->mui);
1325 sol->alphasum = alphasum + sol->mui + sol->muj;
1326 gain = sol->mui * (gi - 0.5 * (Qii*sol->mui + Qij*sol->muj))
1327 + sol->muj * (gj - 0.5 * (Qij*sol->mui + Qjj*sol->muj));
1328 if (gain > bestgain) { bestgain = gain; best = sol; }
1333 double a = 0.0, mu = 0.0;
1334 double ggi = gi - (U - alphai) * Qii + alphaj * Qij;
1335 double ggj = gj - (U - alphai) * Qij + alphaj * Qjj;
1336 solveEdge(a, ggj - ggi, Qii + Qjj - 2.0 * Qij, U, mu);
1337 sol->alphai = U - a;
1339 sol->mui = U - a - alphai;
1340 sol->muj = a - alphaj;
1342 gain = sol->mui * (gi - 0.5 * (Qii*sol->mui + Qij*sol->muj))
1343 + sol->muj * (gj - 0.5 * (Qij*sol->mui + Qjj*sol->muj));
1344 if (gain > bestgain) { bestgain = gain; best = sol; }
1348 alphai = best->alphai;
1349 alphaj = best->alphaj;
1350 alphasum = best->alphasum;
1355 if (alphai + alphaj < 1e-12 *
C) alphai = alphaj = alphasum = 0.0;
1356 if (alphasum > (1.0 - 1e-12) *
C)
1359 if (alphai > (1.0 - 1e-12) *
C) { alphai =
C; alphaj = 0.0; }
1360 else if (alphaj > (1.0 - 1e-12) *
C) { alphai = 0.0; alphaj =
C; }
1365 double kktViolationActive(unsigned int example)
1370 double kktViolationAll(unsigned int example)
1383 double a =
alpha(v);
1387 if (g > ret) ret = g;
1391 if (-g > ret) ret = -g;
1399 std::size_t i, p, pc, v;
1404 bool cangrow = (ex.
varsum <
C);
1406 double down = 1e100;
1407 for (p=0; p<pc; p++)
1411 double a =
alpha(v);
1419 if (g < down) down = g;
1422 if (up - down > ret) ret = up - down;
1423 if (up > ret) ret = up;
1424 if (-down > ret) ret = -down;
1452 unsigned int b, bc = ex.
active;
1455 unsigned int b2, a2;
1456 for (b=0; b<bc; b++)
1458 std::size_t a = ex.
avar[b];
1460 double aa =
alpha(a);
1470 for (b2=0; b2<bc; b2++)
1489 double down = 1e100;
1491 for (b=0; b<bc; b++)
1493 std::size_t v = ex.
avar[b];
1495 double a =
alpha(v);
1497 if (g > up) { i_up = v; up = g; }
1498 if (a > 0.0 && g < down) { i_down = v; down = g; }
1500 if (up - down > ret) { two =
true; ret = up - down; i = i_up; j = i_down; }
1501 if (up > ret) { two =
false; ret = up; i = i_up; }
1502 if (-down > ret) { two =
false; ret = -down; i = i_down; }
1527 if (two || ret == 0.0)
return ret;
1530 std::size_t b, f, pf;
1532 std::size_t ii = vari.
i;
1533 unsigned int pi = vari.
p;
1534 unsigned int yi =
example[ii].y;
1539 double bestgain = 0.0;
1540 double gain_i = gi * gi / di;
1545 double varsum = exa.
varsum;
1546 unsigned int ya = exa.
y;
1547 QpFloatType kiia = k[a];
1552 for (pf=0, b=0; b<row.
size; b++)
1558 double af =
alpha(f);
1560 if ((af > 0.0 && gf < 0.0) || (varsum < C && gf > 0.0))
1563 double gain = gain_i + gf * gf / df;
1564 if (gain > bestgain && f != i) { bestgain = gain; j = f; }
1572 for (; pf<
cardP; pf++)
1576 double af =
alpha(f);
1578 if ((af > 0.0 && gf < 0.0) || (varsum < C && gf > 0.0))
1581 double gain = gain_i + gf * gf / df;
1582 if (gain > bestgain && f != i) { bestgain = gain; j = f; }
1588 for (pf=0, b=0; b<row.
size; b++)
1599 for (; pf<
cardP; pf++)
1617 double aa =
alpha(a);
1619 if (ga > ret && aa <
C)
1624 else if (-ga > ret && aa > 0.0)
1630 if (ret == 0.0)
return ret;
1633 std::size_t b, f, pf;
1635 std::size_t ii = vari.
i;
1636 unsigned int pi = vari.
p;
1637 unsigned int yi =
example[ii].y;
1642 double bestgain = 0.0;
1643 double gain_i = gi * gi / di;
1647 unsigned int ya = exa.
y;
1648 QpFloatType kiia = k[a];
1653 for (pf=0, b=0; b<row.
size; b++)
1658 if (f >= activeVar)
continue;
1659 double af =
alpha(f);
1661 if ((af > 0.0 && gf < 0.0) || (af < C && gf > 0.0))
1664 double gain = gain_i + gf * gf / df;
1665 if (gain > bestgain && f != i) { bestgain = gain; j = f; }
1673 for (; pf<
cardP; pf++)
1676 if (f >= activeVar)
continue;
1677 double af =
alpha(f);
1679 if ((af > 0.0 && gf < 0.0) || (af < C && gf > 0.0))
1682 double gain = gain_i + gf * gf / df;
1683 if (gain > bestgain && f != i) { bestgain = gain; j = f; }
1689 for (pf=0, b=0; b<row.
size; b++)
1700 for (; pf<
cardP; pf++)
1725 throw SHARKEXCEPTION(
"[QpMcDecomp::selectWorkingSetSMO] SMO is implemented only for box constraints");
1734 double bestgain = 0.0;
1737 double aa =
alpha(a);
1739 if (ga > 0.0 && aa <
C)
1741 double gain = ga * ga /
variable[a].diagonal;
1742 if (gain > bestgain)
1747 if (ga > ret) ret = ga;
1749 else if (ga < 0.0 && aa > 0.0)
1751 double gain = ga * ga /
variable[a].diagonal;
1752 if (gain > bestgain)
1757 if (-ga > ret) ret = -ga;
1773 double largest = 0.0;
1785 if (largest < 10.0 * epsilon)
1801 if ((v == 0.0 && g <= 0.0) || (v ==
C && g >= 0.0))
1818 for (a =
activeEx - 1; a >= 0; a--)
1851 if (mu == 0.0)
continue;
1855 unsigned int yv =
example[iv].y;
1856 unsigned int r =
cardP * yv + pv;
1857 std::vector<QpFloatType> q(
examples);
1860 std::size_t a, b, f;
1869 for (b=0; b<row.
size; b++)
1877 for (b=0; b<row.
size; b++)
1882 double upd = (mu) * def * (k);
1900 if (! complete)
shrink(epsilon);
1910 unsigned int iv =
variable[v].index;
1914 std::size_t ih = exv->
active - 1;
1915 std::size_t
h = exv->
avar[ih];
1924 unsigned int ij =
variable[j].index;
1953 std::size_t* pe =
example[e].var;
1954 std::size_t* pj =
example[j].var;
1955 for (v = 0; v <
cardP; v++)
2066 #undef GRADIENT_UPDATE 2067 #undef GAIN_SELECTION_BOX 2068 #undef GAIN_SELECTION_TRIANGLE