37 #ifndef SHARK_ALGORITHMS_QP_QPMCSIMPLEXDECOMP_H 38 #define SHARK_ALGORITHMS_QP_QPMCSIMPLEXDECOMP_H 42 #include <shark/Algorithms/QP/Impl/AnalyticProblems.h> 49 template <
class Matrix>
59 template<
class Problem>
60 double operator()(Problem& problem, std::size_t& i, std::size_t& j){
62 return problem.selectWorkingSet(i,j);
78 RealMatrix
const& linearMat,
100 && linearMat.size1() == kernel.size(),
101 "[QpMcDecomp::QpMcDecomp] dimension conflict" 109 unsigned int y = target.
element(i);
118 for (std::size_t p=0; p<
m_cardP; p++, v++)
123 double Q =
m_M(
m_classes * (y * m_cardP + p) + y, p) * k;
148 return solutionMatrix;
157 return solutionGradientMatrix;
203 unsigned int r =
m_cardP * y + p;
234 unsigned int rv =
m_cardP*yv+pv;
235 unsigned int rw =
m_cardP*yw+pw;
274 for(std::size_t p = 0; p !=
m_cardP; ++p){
279 for(std::size_t p = 0; p !=
m_cardP; ++p){
306 std::pair<std::pair<double,std::size_t>,std::pair<double,std::size_t> > pair =
getSimplexMVP(ex);
307 double up = pair.first.first;
308 double down = pair.second.first;
312 if(down > 0 && ex.
varsum ==
m_C && up-down > 0){
314 for(
int p = pc-1; p >= 0; --p){
318 if(a == 0 && g-down < 0){
321 else if (a ==
m_C && up-g < 0){
325 for(
int q = ex.
active; q >= 0; --q){
334 else if(ex.
varsum == 0 && up < 0){
336 for(
int p = pc-1; p >= 0; --p){
362 if (mu == 0.0)
continue;
367 unsigned int r =
m_cardP * yv + pv;
377 for (std::size_t b=0; b<row.
size; b++)
385 double upd = mu * def * k;
388 std::size_t f = ex.
avar[b];
417 double maxGradient = 0;
418 double maxSimplexGradient = 0;
419 std::size_t maxSimplexExample = 0;
428 std::pair<std::pair<double,std::size_t>,std::pair<double,std::size_t> > pair =
getSimplexMVP(ex);
429 double up = pair.first.first;
430 double down = pair.second.first;
432 if(!canGrow && up - down > maxSimplexGradient){
433 maxSimplexGradient = up-down;
434 maxSimplexExample = e;
436 if (canGrow && up > maxGradient) {
438 i = pair.first.second;
440 if (-down > maxGradient) {
442 i = pair.second.second;
448 std::pair<std::pair<std::size_t,std::size_t> ,
double > boxPair =
maxGainBox(i);
449 double bestGain = boxPair.second;
450 std::pair<std::size_t, std::size_t> best = boxPair.first;
454 if(simplexPairi.second > bestGain){
455 best = simplexPairi.first;
456 bestGain = simplexPairi.second;
459 if(maxSimplexGradient > 0){
460 std::pair<std::pair<std::size_t,std::size_t> ,
double > simplexPairBest =
maxGainSimplex(maxSimplexExample);
461 if(simplexPairBest.second > bestGain){
462 best = simplexPairBest.first;
463 bestGain = simplexPairBest.second;
469 return std::max(maxGradient,maxSimplexGradient);
478 std::pair<std::pair<double,std::size_t>,std::pair<double,std::size_t> > pair =
getSimplexMVP(ex);
479 double up = pair.first.first;
480 double down = pair.second.first;
534 std::pair<std::pair<std::size_t,std::size_t>,
double>
maxGainBox(std::size_t i)
const{
542 return std::make_pair(std::make_pair(i,i),0.0);
547 std::size_t bestj = i;
548 double bestGain = gi * gi / Qii;
556 unsigned int ya = exa.
y;
562 for (std::size_t p=0, b=0; p <
m_cardP; p++)
564 std::size_t j = exa.
var[p];
568 double Qij = def * k[a];
579 double gain = detail::maximumGainQuadratic2D(Qii, Qjj, Qij, gi,gj);
580 if( bestGain < gain){
586 return std::make_pair(std::make_pair(i,bestj),bestGain);
594 std::pair<std::pair<std::size_t,std::size_t>,
double>
maxGainSimplex(std::size_t e)
const{
596 std::size_t pc = ex.
active;
597 unsigned int y = ex.
y;
601 double bestGain = -1e100;
602 std::size_t besti = 0;
603 std::size_t bestj = 0;
616 for(std::size_t p1 = 0; p1 != pc; ++p1){
617 std::size_t i = ex.
avar[p1];
623 if((gi < 0 &&
m_alpha(i) > 0.0) || (gi > 0 && canGrow)){
624 double gain = gi * gi / Qii;
642 for(std::size_t p2 = 0, b=0; p2 !=
m_cardP; ++p2){
643 std::size_t j = ex.
var[p2];
649 double Qij = def * Qee;
663 if(!canGrow && gi > 0 && gj > 0){
667 if(aj > 0 && gi-gj > 0){
668 gainUp = detail::maximumGainQuadratic2DOnLine(Qii, Qjj, Qij, gi,gj);
671 if (ai > 0 &&gj-gi > 0){
672 gainDown = detail::maximumGainQuadratic2DOnLine(Qjj, Qii, Qij, gj,gi);
678 else if(!(gi <= 0 && ai == 0) && !(gj<= 0 && aj == 0)){
679 gain = detail::maximumGainQuadratic2D(Qii, Qjj, Qij, gi,gj);
692 return std::make_pair(std::make_pair(besti,bestj),bestGain);
697 std::pair<double,std::size_t>,
698 std::pair<double,std::size_t>
700 unsigned int pc = ex.
active;
703 std::size_t maxUp = ex.
avar[0];
704 std::size_t minDown = ex.
avar[0];
705 for (std::size_t p = 0; p != pc; p++)
707 std::size_t v = ex.
avar[p];
715 if (a > 0.0 && g < down){
720 return std::make_pair(std::make_pair(up,maxUp),std::make_pair(down,minDown));
724 double& varsum =
m_examples[exampleId].varsum;
726 if(varsum > 1.e-12 &&
m_C-varsum > 1.e-12*
m_C)
731 for(std::size_t p = 0; p !=
m_cardP; ++p){
732 std::size_t varIndex =
m_examples[exampleId].var[p];
738 if(
m_C-varsum < 1.e-14*
m_C)
750 for (std::size_t b=0; b<row.
size; b++){
755 double upd = mu* def * k;
756 for (std::size_t b=0; b<ex.
active; b++)
770 std::size_t ih = exv->
active - 1;
771 std::size_t
h = exv->
avar[ih];
816 for (std::size_t v = 0; v <
m_cardP; v++)
890 template<
class Matrix>
900 bool sumToZero =
false 902 std::size_t classes = bias.size();
903 std::size_t numExamples = m_problem->getNumExamples();
904 std::size_t
cardP = m_problem->cardP();
905 RealVector stepsize(classes, 0.01);
906 RealVector prev(classes,0);
907 RealVector step(classes);
916 RealMatrix dualGradient = m_problem->solutionGradient();
918 RealVector grad(classes,0);
919 for (std::size_t i=0; i<numExamples; i++)
921 std::size_t largestP =
cardP;
922 double largest_value = 0.0;
923 for (std::size_t p=0; p<
cardP; p++)
925 if (dualGradient(i,p) > largest_value)
927 largest_value = dualGradient(i,p);
931 if (largestP < cardP)
933 unsigned int y = m_problem->label(i);
935 for (std::size_t b=0; b != row.
size; b++)
943 grad -=
sum(grad) / classes;
947 for (std::size_t c=0; c<classes; c++)
951 step(c) = -stepsize(c);
953 step(c) = stepsize(c);
955 double gg = prev(c) * grad(c);
966 step -=
sum(step) / classes;
971 performBiasUpdate(step,nu);
978 void performBiasUpdate(
981 std::size_t numExamples = m_problem->getNumExamples();
982 std::size_t
cardP = m_problem->cardP();
983 RealMatrix deltaLinear(numExamples,cardP,0.0);
984 for (std::size_t i=0; i<numExamples; i++){
985 for (std::size_t p=0; p<
cardP; p++){
986 unsigned int y = m_problem->label(i);
989 for (std::size_t b=0; b<row.
size; b++)
995 m_problem->addDeltaLinear(deltaLinear);