510{
511 this->l = l;
512 this->Q = &Q;
513 QD=Q.get_QD();
514 clone(p, p_,l);
515 clone(y, y_,l);
516 clone(alpha,alpha_,l);
517 this->Cp = Cp;
518 this->Cn = Cn;
519 this->eps = eps;
520 unshrink = false;
521
522
523 {
524 alpha_status = new char[l];
525 for(int i=0;i<l;i++)
526 update_alpha_status(i);
527 }
528
529
530 {
531 active_set = new int[l];
532 for(int i=0;i<l;i++)
533 active_set[i] = i;
534 active_size = l;
535 }
536
537
538 {
539 G = new double[l];
540 G_bar = new double[l];
541 int i;
542 for(i=0;i<l;i++)
543 {
544 G[i] = p[i];
545 G_bar[i] = 0;
546 }
547 for(i=0;i<l;i++)
548 if(!is_lower_bound(i))
549 {
550 const Qfloat *Q_i = Q.get_Q(i,l);
551 double alpha_i = alpha[i];
552 int j;
553 for(j=0;j<l;j++)
554 G[j] += alpha_i*Q_i[j];
555 if(is_upper_bound(i))
556 for(j=0;j<l;j++)
557 G_bar[j] += get_C(i) * Q_i[j];
558 }
559 }
560
561
562
563 int iter = 0;
564 int max_iter = max(10000000, l>INT_MAX/100 ? INT_MAX : 100*l);
565 int counter = min(l,1000)+1;
566
567 while(iter < max_iter)
568 {
569
570
571 if(--counter == 0)
572 {
573 counter = min(l,1000);
574 if(shrinking) do_shrinking();
575 if(verbose)
576 info(".");
577 }
578
579 int i,j;
580 if(select_working_set(i,j)!=0)
581 {
582
583 reconstruct_gradient();
584
585 active_size = l;
586 if(verbose)
587 info("*");
588 if(select_working_set(i,j)!=0)
589 break;
590 else
591 counter = 1;
592 }
593
594 ++iter;
595
596
597
598 const Qfloat *Q_i = Q.get_Q(i,active_size);
599 const Qfloat *Q_j = Q.get_Q(j,active_size);
600
601 double C_i = get_C(i);
602 double C_j = get_C(j);
603
604 double old_alpha_i = alpha[i];
605 double old_alpha_j = alpha[j];
606
607 if(y[i]!=y[j])
608 {
609 double quad_coef = QD[i]+QD[j]+2*Q_i[j];
610 if (quad_coef <= 0)
611 quad_coef = TAU;
612 double delta = (-G[i]-G[j])/quad_coef;
613 double diff = alpha[i] - alpha[j];
614 alpha[i] += delta;
615 alpha[j] += delta;
616
617 if(diff > 0)
618 {
619 if(alpha[j] < 0)
620 {
621 alpha[j] = 0;
622 alpha[i] = diff;
623 }
624 }
625 else
626 {
627 if(alpha[i] < 0)
628 {
629 alpha[i] = 0;
630 alpha[j] = -diff;
631 }
632 }
633 if(diff > C_i - C_j)
634 {
635 if(alpha[i] > C_i)
636 {
637 alpha[i] = C_i;
638 alpha[j] = C_i - diff;
639 }
640 }
641 else
642 {
643 if(alpha[j] > C_j)
644 {
645 alpha[j] = C_j;
646 alpha[i] = C_j + diff;
647 }
648 }
649 }
650 else
651 {
652 double quad_coef = QD[i]+QD[j]-2*Q_i[j];
653 if (quad_coef <= 0)
654 quad_coef = TAU;
655 double delta = (G[i]-G[j])/quad_coef;
656 double sum = alpha[i] + alpha[j];
657 alpha[i] -= delta;
658 alpha[j] += delta;
659
660 if(sum > C_i)
661 {
662 if(alpha[i] > C_i)
663 {
664 alpha[i] = C_i;
665 alpha[j] = sum - C_i;
666 }
667 }
668 else
669 {
670 if(alpha[j] < 0)
671 {
672 alpha[j] = 0;
673 alpha[i] = sum;
674 }
675 }
676 if(sum > C_j)
677 {
678 if(alpha[j] > C_j)
679 {
680 alpha[j] = C_j;
681 alpha[i] = sum - C_j;
682 }
683 }
684 else
685 {
686 if(alpha[i] < 0)
687 {
688 alpha[i] = 0;
689 alpha[j] = sum;
690 }
691 }
692 }
693
694
695
696 double delta_alpha_i = alpha[i] - old_alpha_i;
697 double delta_alpha_j = alpha[j] - old_alpha_j;
698
699 for(int k=0;k<active_size;k++)
700 {
701 G[k] += Q_i[k]*delta_alpha_i + Q_j[k]*delta_alpha_j;
702 }
703
704
705
706 {
707 bool ui = is_upper_bound(i);
708 bool uj = is_upper_bound(j);
709 update_alpha_status(i);
710 update_alpha_status(j);
711 int k;
712 if(ui != is_upper_bound(i))
713 {
714 Q_i = Q.get_Q(i,l);
715 if(ui)
716 for(k=0;k<l;k++)
717 G_bar[k] -= C_i * Q_i[k];
718 else
719 for(k=0;k<l;k++)
720 G_bar[k] += C_i * Q_i[k];
721 }
722
723 if(uj != is_upper_bound(j))
724 {
725 Q_j = Q.get_Q(j,l);
726 if(uj)
727 for(k=0;k<l;k++)
728 G_bar[k] -= C_j * Q_j[k];
729 else
730 for(k=0;k<l;k++)
731 G_bar[k] += C_j * Q_j[k];
732 }
733 }
734 }
735
736 if(iter >= max_iter)
737 {
738 if(active_size < l)
739 {
740
741 reconstruct_gradient();
742 active_size = l;
743 if(verbose)
744 info("*");
745 }
746 info("\nWARNING: reaching max number of iterations");
747 }
748
749
750
751 si->rho = calculate_rho();
752
753
754 {
755 double v = 0;
756 int i;
757 for(i=0;i<l;i++)
758 v += alpha[i] * (G[i] + p[i]);
759
760 si->obj = v/2;
761 }
762
763
764 {
765 for(int i=0;i<l;i++)
766 alpha_[active_set[i]] = alpha[i];
767 }
768
769
770
771
772
773
774
775
776
777 si->upper_bound_p = Cp;
778 si->upper_bound_n = Cn;
779
780 if(verbose)
781 info("\noptimization finished, #iter = %d\n",iter);
782
783 delete[] p;
784 delete[] y;
785 delete[] alpha;
786 delete[] alpha_status;
787 delete[] active_set;
788 delete[] G;
789 delete[] G_bar;
790}