32 #ifndef SHARK_ALGORITHMS_QP_BOXCONSTRAINEDPROBLEMS_H 33 #define SHARK_ALGORITHMS_QP_BOXCONSTRAINEDPROBLEMS_H 36 #include <shark/Algorithms/QP/Impl/AnalyticProblems.h> 44 template<
class Problem>
45 double operator()(Problem& problem, std::size_t& i, std::size_t& j){
48 double largestGradient = 0;
49 double secondLargestGradient = 0;
51 for (std::size_t a = 0; a < problem.active(); a++){
52 double g = problem.gradient(a);
53 if (!problem.isUpperBound(a) && g > secondLargestGradient){
54 secondLargestGradient = g;
57 if (!problem.isLowerBound(a) && -g > secondLargestGradient){
58 secondLargestGradient = -g;
61 if(secondLargestGradient > largestGradient){
62 std::swap(secondLargestGradient,largestGradient);
66 if(secondLargestGradient == 0)
68 return largestGradient;
78 template<
class Problem>
79 double operator()(Problem& problem, std::size_t& i, std::size_t& j){
81 double value = criterium(problem, i,j);
91 template<
class Problem>
92 double operator()(Problem& problem, std::size_t& i, std::size_t& j){
95 double maxGrad = firstOrder(problem,i,j);
99 double gi = problem.gradient(i);
100 typename Problem::QpFloatType* q = problem.quadratic().row(i, 0, problem.active());
101 double Qii = problem.diagonal(i);
104 double maxGain = 0.0;
105 for (std::size_t a=0; a<problem.active(); a++)
107 if (a == i)
continue;
108 double ga = problem.gradient(a);
110 (!problem.isLowerBound(a) && ga < 0.0)
111 || (!problem.isUpperBound(a) && ga > 0.0)
114 double Qaa = problem.diagonal(a);
115 double gain = detail::maximumGainQuadratic2D(Qii,Qaa,Qia,gi,ga);
136 template<
class SVMProblem>
146 , m_gradient(problem.linear)
147 , m_active (problem.dimensions())
148 , m_alphaStatus(problem.dimensions(),
AlphaFree){
150 for (std::size_t i=0; i != dimensions(); i++){
153 QpFloatType* q = quadratic().row(i, 0, dimensions());
154 for (std::size_t a=0; a < dimensions(); a++)
155 m_gradient(a) -= q[a] * v;
157 updateAlphaStatus(i);
161 return m_problem.dimensions();
181 return isUpperBound(i) && isLowerBound(i);
186 return m_problem.quadratic;
190 return m_problem.linear(i);
194 return m_problem.alpha(i);
198 return m_problem.diagonal(i);
202 return m_gradient(i);
206 return m_problem.permutation[i];
210 RealVector alpha(dimensions());
211 for (std::size_t i=0; i<dimensions(); i++)
212 alpha(m_problem.permutation[i]) = m_problem.alpha(i);
222 QpFloatType* q = quadratic().row(i, 0, active());
226 double mu = -alpha(i);
227 detail::solveQuadraticEdge(m_problem.alpha(i),gradient(i),diagonal(i),boxMin(i),boxMax(i));
231 for (std::size_t a = 0; a < active(); a++)
232 m_gradient(a) -= mu * q[a];
234 updateAlphaStatus(i);
238 double Li = boxMin(i);
239 double Ui = boxMax(i);
240 double Lj = boxMin(j);
241 double Uj = boxMax(j);
244 QpFloatType* qi = quadratic().row(i, 0, active());
245 QpFloatType* qj = quadratic().row(j, 0, active());
249 double mui = -alpha(i);
250 double muj = -alpha(j);
251 detail::solveQuadratic2DBox(m_problem.alpha(i), m_problem.alpha(j),
252 m_gradient(i), m_gradient(j),
253 diagonal(i), qi[j], diagonal(j),
260 for (std::size_t a = 0; a < active(); a++)
261 m_gradient(a) -= mui * qi[a] + muj * qj[a];
263 updateAlphaStatus(i);
264 updateAlphaStatus(j);
269 return 0.5*
inner_prod(m_gradient+m_problem.linear,m_problem.alpha);
279 double alphai = alpha(i);
280 m_problem.alpha(i) = 0;
282 QpFloatType* qi = quadratic().row(i, 0, active());
283 for (std::size_t a = 0; a < active(); a++)
284 m_gradient(a) += alphai * qi[a];
290 updateAlphaStatus(i);
300 m_problem.flipCoordinates(i, j);
301 std::swap( m_gradient[i], m_gradient[j]);
302 std::swap( m_alphaStatus[i], m_alphaStatus[j]);
307 m_gradient(i) -= linear(i);
308 m_gradient(i) += newValue;
309 m_problem.linear(i) = newValue;
313 double maxViolation = 0.0;
314 for(std::size_t i = 0; i != dimensions(); ++i){
315 if(isDeactivated(i))
continue;
316 if(!isUpperBound(i)){
317 maxViolation =
std::max(maxViolation, gradient(i));
319 if(!isLowerBound(i)){
320 maxViolation =
std::max(maxViolation, -gradient(i));
339 if(m_problem.alpha(i) == boxMax(i))
341 if(m_problem.alpha(i) == boxMin(i))
348 template<
class Problem>
359 , m_isUnshrinked(false)
361 , m_gradientEdge(problem.linear){}
363 using base_type::alpha;
364 using base_type::gradient;
365 using base_type::linear;
366 using base_type::active;
367 using base_type::dimensions;
368 using base_type::quadratic;
369 using base_type::isLowerBound;
370 using base_type::isUpperBound;
371 using base_type::boxMin;
372 using base_type::boxMax;
375 double aiOld = alpha(i);
376 double ajOld = alpha(j);
378 base_type::updateSMO(i,j);
379 double ai = alpha(i);
380 double aj = alpha(j);
383 updateGradientEdge(i,aiOld,ai);
384 updateGradientEdge(j,ajOld,aj);
388 if(!m_shrink)
return false;
392 getMaxKKTViolations(largestUp,smallestDown,active());
395 if (!m_isUnshrinked && (largestUp - smallestDown < 10.0 * epsilon))
399 getMaxKKTViolations(largestUp,smallestDown,dimensions());
402 smallestDown =
std::min(smallestDown,0.0);
403 largestUp =
std::max(largestUp,0.0);
404 for (std::size_t a = this->active(); a > 0; --a){
405 if(testShrinkVariable(a-1,largestUp,smallestDown)){
406 flipCoordinates(a-1,active()-1);
415 if (active() == dimensions())
return;
416 m_isUnshrinked =
true;
425 for (std::size_t a = active(); a < dimensions(); a++)
426 this->m_gradient(a) = m_gradientEdge(a);
428 for (std::size_t i = 0; i < active(); i++)
431 if (isUpperBound(i) || isLowerBound(i))
continue;
433 QpFloatType* q = quadratic().row(i, 0, dimensions());
434 for (std::size_t a = active(); a < dimensions(); a++)
435 this->m_gradient(a) -= alpha(i) * q[a] ;
438 this->m_active = dimensions();
442 m_shrink = shrinking;
449 base_type::scaleBoxConstraints(factor,variableScalingFactor);
450 if(factor != variableScalingFactor){
451 for(std::size_t i = 0; i != dimensions(); ++i){
452 m_gradientEdge(i) = linear(i);
456 for(std::size_t i = 0; i != dimensions(); ++i){
457 m_gradientEdge(i) -= linear(i);
458 m_gradientEdge(i) *= factor;
459 m_gradientEdge(i) += linear(i);
466 double alphai = alpha(i);
467 base_type::deactivateVariable(i);
468 updateGradientEdge(i,alphai,0.0);
473 m_gradientEdge(i) -= linear(i);
474 m_gradientEdge(i) += newValue;
475 base_type::setLinear(i,newValue);
480 base_type::flipCoordinates(i,j);
481 std::swap( m_gradientEdge[i], m_gradientEdge[j]);
490 if(!m_shrink || oldAlpha==newAlpha)
return;
491 bool isInsideOld = oldAlpha > boxMin(i) && oldAlpha < boxMax(i);
492 bool isInsideNew = newAlpha > boxMin(i) && newAlpha < boxMax(i);
495 if( (oldAlpha == 0 || isInsideOld) && isInsideNew )
506 QpFloatType* q = quadratic().row(i, 0, dimensions());
507 for(std::size_t a = 0; a != dimensions(); ++a){
508 m_gradientEdge(a) -= diff*q[a];
513 bool testShrinkVariable(std::size_t a,
double largestUp,
double smallestDown)
const{
515 ( isLowerBound(a) && gradient(a) < smallestDown)
516 || ( isUpperBound(a) && gradient(a) >largestUp)
525 void getMaxKKTViolations(
double& largestUp,
double& smallestDown, std::size_t maxIndex){
527 smallestDown = 1e100;
528 for (std::size_t a = 0; a < maxIndex; a++){
529 if (!isLowerBound(a))
530 smallestDown =
std::min(smallestDown,gradient(a));
531 if (!isUpperBound(a))
532 largestUp =
std::max(largestUp,gradient(a));
542 RealVector m_gradientEdge;