30 #ifndef SHARK_ALGORITHMS_QP_SVMPROBLEMS_H 31 #define SHARK_ALGORITHMS_QP_SVMPROBLEMS_H 47 template<
class Problem>
48 double operator()(Problem& problem, std::size_t& i, std::size_t& j)
50 double largestUp = -1e100;
51 double smallestDown = 1e100;
53 for (std::size_t a=0; a < problem.active(); a++)
55 double ga = problem.gradient(a);
56 if (!problem.isUpperBound(a))
64 if (!problem.isLowerBound(a))
66 if (ga < smallestDown)
75 return largestUp - smallestDown;
90 template<
class Problem>
91 double operator()(Problem& problem, std::size_t& i, std::size_t& j)
96 double smallestDown = 1e100;
97 double largestUp = -1e100;
99 for (std::size_t a=0; a < problem.active(); a++)
101 double ga = problem.gradient(a);
103 if (!problem.isUpperBound(a))
112 if (largestUp == -1e100)
return 0.0;
115 typename Problem::QpFloatType* q = problem.quadratic().row(i, 0, problem.active());
117 for (std::size_t a = 0; a < problem.active(); a++){
118 double ga = problem.gradient(a);
119 if (!problem.isLowerBound(a))
121 smallestDown=
std::min(smallestDown,ga);
123 double gain = detail::maximumGainQuadratic2DOnLine(
137 if (best == 0.0)
return 0.0;
140 return largestUp - smallestDown;
156 template<
class Problem>
157 double operator()(Problem& problem, std::size_t& i, std::size_t& j)
159 if (smallProblem || useLibSVM || isInCorner(problem))
162 if(!smallProblem &&
sqr(problem.active()) < problem.quadratic().getMaxCacheSize())
165 double value = libSVMSelection(problem,i, j);
171 MGStep besti = selectMGVariable(problem,last_i);
172 if(besti.violation == 0.0)
174 MGStep bestj = selectMGVariable(problem,last_j);
176 if(bestj.gain > besti.gain){
183 if (problem.gradient(i) < problem.gradient(j))
187 return besti.violation;
192 smallProblem =
false;
196 template<
class Problem>
197 bool isInCorner(Problem
const& problem)
const{
198 double Li = problem.boxMin(last_i);
199 double Ui = problem.boxMax(last_i);
200 double Lj = problem.boxMin(last_j);
201 double Uj = problem.boxMax(last_j);
202 double eps_i = 1e-8 * (Ui - Li);
203 double eps_j = 1e-8 * (Uj - Lj);
205 if ((problem.alpha(last_i) <= Li + eps_i || problem.alpha(last_i) >= Ui - eps_i)
206 && ((problem.alpha(last_j) <= Lj + eps_j || problem.alpha(last_j) >= Uj - eps_j)))
217 template<
class Problem>
218 MGStep selectMGVariable(Problem& problem,std::size_t i)
const{
221 std::size_t bestIndex = 0;
224 double largestUp = -1e100;
225 double smallestDown = 1e100;
228 typename Problem::QpFloatType* q = problem.quadratic().row(i, 0, problem.active());
229 double ab = problem.alpha(i);
230 double db = problem.diagonal(i);
231 double Lb = problem.boxMin(i);
232 double Ub = problem.boxMax(i);
233 double gb = problem.gradient(i);
234 for (std::size_t a = 0; a < problem.active(); a++)
236 double ga = problem.gradient(a);
238 if (!problem.isUpperBound(a))
240 if (!problem.isLowerBound(a))
241 smallestDown =
std::min(smallestDown,ga);
243 if (a == i)
continue;
245 double denominator = (problem.diagonal(a) + db - 2.0 * q[a]);
246 double mu_max = (ga - gb) / denominator;
253 double aa = problem.alpha(a);
254 double La = problem.boxMin(a);
255 double Ua = problem.boxMax(a);
256 double mu_star = mu_max;
257 if (aa + mu_star < La) mu_star = La - aa;
258 else if (mu_star + aa > Ua) mu_star = Ua - aa;
259 if (ab - mu_star < Lb) mu_star = ab - Lb;
260 else if (ab - mu_star > Ub) mu_star = ab - Ub;
262 double gain = mu_star * (2.0 * mu_max - mu_star) * denominator;
272 step.violation= largestUp-smallestDown;
273 step.index = bestIndex;
287 template<
class Problem>
296 , m_gradient(problem.linear)
297 , m_active(problem.dimensions())
298 , m_alphaStatus(problem.dimensions(),
AlphaFree){
300 for (std::size_t i=0; i != dimensions(); i++){
303 QpFloatType* q = quadratic().row(i, 0, dimensions());
304 for (std::size_t a=0; a < dimensions(); a++)
305 m_gradient(a) -= q[a] * v;
307 updateAlphaStatus(i);
311 return m_problem.dimensions();
333 return m_problem.quadratic;
337 return m_problem.linear(i);
341 return m_problem.alpha(i);
345 return m_problem.diagonal(i);
349 return m_gradient(i);
353 return m_problem.permutation[i];
357 RealVector alpha(dimensions());
358 for (std::size_t i=0; i<dimensions(); i++)
359 alpha(m_problem.permutation[i]) = m_problem.alpha(i);
368 QpFloatType* qi = quadratic().row(i, 0, active());
371 double numerator = gradient(i) - gradient(j);
372 double denominator = diagonal(i) + diagonal(j) - 2.0 * qi[j];
373 denominator =
std::max(denominator,1.e-12);
374 double mu = numerator/denominator;
382 return 0.5*
inner_prod(m_gradient+m_problem.linear,m_problem.alpha);
398 for (std::size_t j=0; j<dimensions(); j++){
399 if(alpha(i) == 0.0)
break;
402 applyStep(i,j, -alpha(i));
411 updateAlphaStatus(i);
421 m_problem.flipCoordinates(i, j);
422 std::swap( m_gradient[i], m_gradient[j]);
423 std::swap( m_alphaStatus[i], m_alphaStatus[j]);
428 m_problem.scaleBoxConstraints(factor,variableScalingFactor);
429 for(std::size_t i = 0; i != this->dimensions(); ++i){
432 m_gradient(i) -= linear(i);
433 m_gradient(i) *= variableScalingFactor;
434 m_gradient(i) += linear(i);
435 updateAlphaStatus(i);
441 m_gradient(i) -= linear(i);
442 m_gradient(i) += newValue;
443 m_problem.linear(i) = newValue;
447 double largestUp = -1e100;
448 double smallestDown = 1e100;
449 for (std::size_t a = 0; a < active(); a++){
450 if (!isLowerBound(a))
451 smallestDown =
std::min(smallestDown,gradient(a));
452 if (!isUpperBound(a))
453 largestUp =
std::max(largestUp,gradient(a));
455 return largestUp - smallestDown;
477 virtual void applyStep(std::size_t i, std::size_t j,
double step){
481 double Ui = boxMax(i);
482 double Lj = boxMin(j);
483 double aiOld = m_problem.alpha(i);
484 double ajOld = m_problem.alpha(j);
485 double& ai = m_problem.alpha(i);
486 double& aj = m_problem.alpha(j);
487 if (step >=
std::min(Ui - ai, aj - Lj))
489 if (Ui - ai > aj - Lj)
495 else if (Ui - ai < aj - Lj)
514 if(ai == aiOld && aj == ajOld)
return;
517 QpFloatType* qi = quadratic().row(i, 0, active());
518 QpFloatType* qj = quadratic().row(j, 0, active());
519 for (std::size_t a = 0; a < active(); a++)
520 m_gradient(a) -= step * qi[a] - step * qj[a];
523 updateAlphaStatus(i);
524 updateAlphaStatus(j);
530 if(m_problem.alpha(i) == boxMax(i))
532 if(m_problem.alpha(i) == boxMin(i))
537 template<
class Problem>
548 , m_isUnshrinked(false)
550 , m_gradientEdge(problem.linear){}
552 using base_type::alpha;
553 using base_type::gradient;
554 using base_type::linear;
555 using base_type::active;
556 using base_type::dimensions;
557 using base_type::quadratic;
558 using base_type::isLowerBound;
559 using base_type::isUpperBound;
560 using base_type::boxMin;
561 using base_type::boxMax;
564 if(!m_shrink)
return false;
568 getMaxKKTViolations(largestUp,smallestDown,active());
571 if (!m_isUnshrinked && (largestUp - smallestDown < 10.0 * epsilon))
575 getMaxKKTViolations(largestUp,smallestDown,dimensions());
578 for (std::size_t a = this->active(); a > 0; --a){
579 if(testShrinkVariable(a-1,largestUp,smallestDown))
580 this->shrinkVariable(a-1);
587 if (active() == dimensions())
return;
588 m_isUnshrinked =
true;
597 for (std::size_t a = active(); a < dimensions(); a++)
598 this->m_gradient(a) = m_gradientEdge(a);
600 for (std::size_t i = 0; i < active(); i++)
603 if (isUpperBound(i) || isLowerBound(i))
continue;
605 QpFloatType* q = quadratic().row(i, 0, dimensions());
606 for (std::size_t a = active(); a < dimensions(); a++)
607 this->m_gradient(a) -= alpha(i) * q[a] ;
610 this->m_active = dimensions();
614 m_shrink = shrinking;
621 base_type::scaleBoxConstraints(factor,variableScalingFactor);
622 if(factor != variableScalingFactor){
623 for(std::size_t i = 0; i != dimensions(); ++i){
624 m_gradientEdge(i) = linear(i);
628 for(std::size_t i = 0; i != dimensions(); ++i){
629 m_gradientEdge(i) -= linear(i);
630 m_gradientEdge(i) *= factor;
631 m_gradientEdge(i) += linear(i);
638 m_gradientEdge(i) -= linear(i);
639 m_gradientEdge(i) += newValue;
640 base_type::setLinear(i,newValue);
651 virtual void applyStep(std::size_t i, std::size_t j,
double step){
654 double aiOld = alpha(i);
655 double ajOld = alpha(j);
657 base_type::applyStep(i,j,step);
658 double ai = alpha(i);
659 double aj = alpha(j);
660 if(!m_shrink || ai == aiOld)
return;
663 updateGradientEdge(i,aiOld,ai);
664 updateGradientEdge(j,ajOld,aj);
668 void updateGradientEdge(std::size_t i,
double oldAlpha,
double newAlpha){
670 bool isInsideOld = oldAlpha > boxMin(i) && oldAlpha < boxMax(i);
671 bool isInsideNew = newAlpha > boxMin(i) && newAlpha < boxMax(i);
674 if( (oldAlpha == 0 || isInsideOld) && isInsideNew )
685 QpFloatType* q = quadratic().row(i, 0, dimensions());
686 for(std::size_t a = 0; a != dimensions(); ++a){
687 m_gradientEdge(a) -= diff*q[a];
691 void shrinkVariable(std::size_t i){
693 base_type::flipCoordinates(i,active()-1);
694 std::swap( m_gradientEdge[i], m_gradientEdge[active()-1]);
699 bool testShrinkVariable(std::size_t a,
double largestUp,
double smallestDown)
const{
701 ( isLowerBound(a) && gradient(a) < smallestDown)
702 || ( isUpperBound(a) && gradient(a) >largestUp)
711 void getMaxKKTViolations(
double& largestUp,
double& smallestDown, std::size_t maxIndex){
713 smallestDown = 1e100;
714 for (std::size_t a = 0; a < maxIndex; a++){
715 if (!isLowerBound(a))
716 smallestDown =
std::min(smallestDown,gradient(a));
717 if (!isUpperBound(a))
718 largestUp =
std::max(largestUp,gradient(a));
728 RealVector m_gradientEdge;