Solver Class Reference

Inheritance diagram for Solver:

Collaboration diagram for Solver:

Classes
struct	SolutionInfo

Public Member Functions
	Solver ()
virtual	~Solver ()
void	Solve (int l, const QMatrix &Q, const double *p_, const schar *y_, double *alpha_, double Cp, double Cn, double eps, SolutionInfo *si, int shrinking)

Protected Types
enum	{ LOWER_BOUND, UPPER_BOUND, FREE }

Protected Member Functions
double	get_C (int i)
void	update_alpha_status (int i)
bool	is_upper_bound (int i)
bool	is_lower_bound (int i)
bool	is_free (int i)
void	swap_index (int i, int j)
void	reconstruct_gradient ()
virtual int	select_working_set (int &i, int &j)
virtual double	calculate_rho ()
virtual void	do_shrinking ()

Protected Attributes
int	active_size
schar *	y
double *	G
char *	alpha_status
double *	alpha
const QMatrix *	Q
const double *	QD
double	eps
double	Cp
double	Cn
double *	p
int *	active_set
double *	G_bar
int	l
bool	unshrink

Detailed Description

Definition at line 399 of file svm.cpp.

Member Enumeration Documentation

anonymous enum

anonymous enum

protected

Enumerator
LOWER_BOUND
UPPER_BOUND
FREE

Definition at line 419 of file svm.cpp.

{ LOWER_BOUND, UPPER_BOUND, FREE };

Constructor & Destructor Documentation

Solver()

Solver::Solver ( )

inline

Definition at line 401 of file svm.cpp.

{};

~Solver()

virtual Solver::~Solver ( )

inlinevirtual

Definition at line 402 of file svm.cpp.

{};

Member Function Documentation

calculate_rho()

double Solver::calculate_rho ( )

protectedvirtual

Definition at line 972 of file svm.cpp.

{
    double r;
    int nr_free = 0;
    double ub = INF, lb = -INF, sum_free = 0;
    for(int i=0;i<active_size;i++)
    {
        double yG = y[i]*G[i];

        if(is_upper_bound(i))
        {
            if(y[i]==-1)
                ub = min(ub,yG);
            else
                lb = max(lb,yG);
        }
        else if(is_lower_bound(i))
        {
            if(y[i]==+1)
                ub = min(ub,yG);
            else
                lb = max(lb,yG);
        }
        else
        {
            ++nr_free;
            sum_free += yG;
        }
    }

    if(nr_free>0)
        r = sum_free/nr_free;
    else
        r = (ub+lb)/2;

    return r;
}

do_shrinking()

void Solver::do_shrinking ( )

protectedvirtual

Definition at line 911 of file svm.cpp.

{
    int i;
    double Gmax1 = -INF;        // max { -y_i * grad(f)_i | i in I_up(\alpha) }
    double Gmax2 = -INF;        // max { y_i * grad(f)_i | i in I_low(\alpha) }

    // find maximal violating pair first
    for(i=0;i<active_size;i++)
    {
        if(y[i]==+1)
        {
            if(!is_upper_bound(i))
            {
                if(-G[i] >= Gmax1)
                    Gmax1 = -G[i];
            }
            if(!is_lower_bound(i))
            {
                if(G[i] >= Gmax2)
                    Gmax2 = G[i];
            }
        }
        else
        {
            if(!is_upper_bound(i))
            {
                if(-G[i] >= Gmax2)
                    Gmax2 = -G[i];
            }
            if(!is_lower_bound(i))
            {
                if(G[i] >= Gmax1)
                    Gmax1 = G[i];
            }
        }
    }

    if(unshrink == false && Gmax1 + Gmax2 <= eps*10)
    {
        unshrink = true;
        reconstruct_gradient();
        active_size = l;
        info("*");
    }

    for(i=0;i<active_size;i++)
        if (be_shrunk(i, Gmax1, Gmax2))
        {
            active_size--;
            while (active_size > i)
            {
                if (!be_shrunk(active_size, Gmax1, Gmax2))
                {
                    swap_index(i,active_size);
                    break;
                }
                active_size--;
            }
        }
}

get_C()

double Solver::get_C ( int  i )

inlineprotected

Definition at line 432 of file svm.cpp.

    {
        return (y[i] > 0)? Cp : Cn;
    }

is_free()

bool Solver::is_free ( int  i )

inlineprotected

Definition at line 446 of file svm.cpp.

{ return alpha_status[i] == FREE; }

is_lower_bound()

bool Solver::is_lower_bound ( int  i )

inlineprotected

Definition at line 445 of file svm.cpp.

{ return alpha_status[i] == LOWER_BOUND; }

is_upper_bound()

bool Solver::is_upper_bound ( int  i )

inlineprotected

Definition at line 444 of file svm.cpp.

{ return alpha_status[i] == UPPER_BOUND; }

reconstruct_gradient()

void Solver::reconstruct_gradient ( )

protected

Definition at line 468 of file svm.cpp.

{
    // reconstruct inactive elements of G from G_bar and free variables

    if(active_size == l) return;

    int i,j;
    int nr_free = 0;

    for(j=active_size;j<l;j++)
        G[j] = G_bar[j] + p[j];

    for(j=0;j<active_size;j++)
        if(is_free(j))
            nr_free++;

    if(2*nr_free < active_size)
        info("\nWARNING: using -h 0 may be faster\n");

    if (nr_free*l > 2*active_size*(l-active_size))
    {
        for(i=active_size;i<l;i++)
        {
            const Qfloat *Q_i = Q->get_Q(i,active_size);
            for(j=0;j<active_size;j++)
                if(is_free(j))
                    G[i] += alpha[j] * Q_i[j];
        }
    }
    else
    {
        for(i=0;i<active_size;i++)
            if(is_free(i))
            {
                const Qfloat *Q_i = Q->get_Q(i,l);
                double alpha_i = alpha[i];
                for(j=active_size;j<l;j++)
                    G[j] += alpha_i * Q_i[j];
            }
    }
}

select_working_set()

int Solver::select_working_set ( int &  i, int &  j  )

protectedvirtual

Definition at line 792 of file svm.cpp.

{
    // return i,j such that
    // i: maximizes -y_i * grad(f)_i, i in I_up(\alpha)
    // j: minimizes the decrease of obj value
    //    (if quadratic coefficeint <= 0, replace it with tau)
    //    -y_j*grad(f)_j < -y_i*grad(f)_i, j in I_low(\alpha)

    double Gmax = -INF;
    double Gmax2 = -INF;
    int Gmax_idx = -1;
    int Gmin_idx = -1;
    double obj_diff_min = INF;

    for(int t=0;t<active_size;t++)
        if(y[t]==+1)
        {
            if(!is_upper_bound(t))
                if(-G[t] >= Gmax)
                {
                    Gmax = -G[t];
                    Gmax_idx = t;
                }
        }
        else
        {
            if(!is_lower_bound(t))
                if(G[t] >= Gmax)
                {
                    Gmax = G[t];
                    Gmax_idx = t;
                }
        }

    int i = Gmax_idx;
    const Qfloat *Q_i = NULL;
    if(i != -1) // NULL Q_i not accessed: Gmax=-INF if i=-1
        Q_i = Q->get_Q(i,active_size);

    for(int j=0;j<active_size;j++)
    {
        if(y[j]==+1)
        {
            if (!is_lower_bound(j))
            {
                double grad_diff=Gmax+G[j];
                if (G[j] >= Gmax2)
                    Gmax2 = G[j];
                if (grad_diff > 0)
                {
                    double obj_diff;
                    double quad_coef = QD[i]+QD[j]-2.0*y[i]*Q_i[j];
                    if (quad_coef > 0)
                        obj_diff = -(grad_diff*grad_diff)/quad_coef;
                    else
                        obj_diff = -(grad_diff*grad_diff)/TAU;

                    if (obj_diff <= obj_diff_min)
                    {
                        Gmin_idx=j;
                        obj_diff_min = obj_diff;
                    }
                }
            }
        }
        else
        {
            if (!is_upper_bound(j))
            {
                double grad_diff= Gmax-G[j];
                if (-G[j] >= Gmax2)
                    Gmax2 = -G[j];
                if (grad_diff > 0)
                {
                    double obj_diff;
                    double quad_coef = QD[i]+QD[j]+2.0*y[i]*Q_i[j];
                    if (quad_coef > 0)
                        obj_diff = -(grad_diff*grad_diff)/quad_coef;
                    else
                        obj_diff = -(grad_diff*grad_diff)/TAU;

                    if (obj_diff <= obj_diff_min)
                    {
                        Gmin_idx=j;
                        obj_diff_min = obj_diff;
                    }
                }
            }
        }
    }

    if(Gmax+Gmax2 < eps || Gmin_idx == -1)
        return 1;

    out_i = Gmax_idx;
    out_j = Gmin_idx;
    return 0;
}

Solve()

void Solver::Solve	(	int	l,
		const QMatrix &	Q,
		const double *	p_,
		const schar *	y_,
		double *	alpha_,
		double	Cp,
		double	Cn,
		double	eps,
		SolutionInfo *	si,
		int	shrinking
	)

Definition at line 510 of file svm.cpp.

{
    this->l = l;
    this->Q = &Q;
    QD=Q.get_QD();
    clone(p, p_,l);
    clone(y, y_,l);
    clone(alpha,alpha_,l);
    this->Cp = Cp;
    this->Cn = Cn;
    this->eps = eps;
    unshrink = false;

    // initialize alpha_status
    {
        alpha_status = new char[l];
        for(int i=0;i<l;i++)
            update_alpha_status(i);
    }

    // initialize active set (for shrinking)
    {
        active_set = new int[l];
        for(int i=0;i<l;i++)
            active_set[i] = i;
        active_size = l;
    }

    // initialize gradient
    {
        G = new double[l];
        G_bar = new double[l];
        int i;
        for(i=0;i<l;i++)
        {
            G[i] = p[i];
            G_bar[i] = 0;
        }
        for(i=0;i<l;i++)
            if(!is_lower_bound(i))
            {
                const Qfloat *Q_i = Q.get_Q(i,l);
                double alpha_i = alpha[i];
                int j;
                for(j=0;j<l;j++)
                    G[j] += alpha_i*Q_i[j];
                if(is_upper_bound(i))
                    for(j=0;j<l;j++)
                        G_bar[j] += get_C(i) * Q_i[j];
            }
    }

    // optimization step

    int iter = 0;
    int max_iter = max(10000000, l>INT_MAX/100 ? INT_MAX : 100*l);
    int counter = min(l,1000)+1;

    while(iter < max_iter)
    {
        // show progress and do shrinking

        if(--counter == 0)
        {
            counter = min(l,1000);
            if(shrinking) do_shrinking();
            info(".");
        }

        int i,j;
        if(select_working_set(i,j)!=0)
        {
            // reconstruct the whole gradient
            reconstruct_gradient();
            // reset active set size and check
            active_size = l;
            info("*");
            if(select_working_set(i,j)!=0)
                break;
            else
                counter = 1;    // do shrinking next iteration
        }

        ++iter;

        // update alpha[i] and alpha[j], handle bounds carefully

        const Qfloat *Q_i = Q.get_Q(i,active_size);
        const Qfloat *Q_j = Q.get_Q(j,active_size);

        double C_i = get_C(i);
        double C_j = get_C(j);

        double old_alpha_i = alpha[i];
        double old_alpha_j = alpha[j];

        if(y[i]!=y[j])
        {
            double quad_coef = QD[i]+QD[j]+2*Q_i[j];
            if (quad_coef <= 0)
                quad_coef = TAU;
            double delta = (-G[i]-G[j])/quad_coef;
            double diff = alpha[i] - alpha[j];
            alpha[i] += delta;
            alpha[j] += delta;

            if(diff > 0)
            {
                if(alpha[j] < 0)
                {
                    alpha[j] = 0;
                    alpha[i] = diff;
                }
            }
            else
            {
                if(alpha[i] < 0)
                {
                    alpha[i] = 0;
                    alpha[j] = -diff;
                }
            }
            if(diff > C_i - C_j)
            {
                if(alpha[i] > C_i)
                {
                    alpha[i] = C_i;
                    alpha[j] = C_i - diff;
                }
            }
            else
            {
                if(alpha[j] > C_j)
                {
                    alpha[j] = C_j;
                    alpha[i] = C_j + diff;
                }
            }
        }
        else
        {
            double quad_coef = QD[i]+QD[j]-2*Q_i[j];
            if (quad_coef <= 0)
                quad_coef = TAU;
            double delta = (G[i]-G[j])/quad_coef;
            double sum = alpha[i] + alpha[j];
            alpha[i] -= delta;
            alpha[j] += delta;

            if(sum > C_i)
            {
                if(alpha[i] > C_i)
                {
                    alpha[i] = C_i;
                    alpha[j] = sum - C_i;
                }
            }
            else
            {
                if(alpha[j] < 0)
                {
                    alpha[j] = 0;
                    alpha[i] = sum;
                }
            }
            if(sum > C_j)
            {
                if(alpha[j] > C_j)
                {
                    alpha[j] = C_j;
                    alpha[i] = sum - C_j;
                }
            }
            else
            {
                if(alpha[i] < 0)
                {
                    alpha[i] = 0;
                    alpha[j] = sum;
                }
            }
        }

        // update G

        double delta_alpha_i = alpha[i] - old_alpha_i;
        double delta_alpha_j = alpha[j] - old_alpha_j;

        for(int k=0;k<active_size;k++)
        {
            G[k] += Q_i[k]*delta_alpha_i + Q_j[k]*delta_alpha_j;
        }

        // update alpha_status and G_bar

        {
            bool ui = is_upper_bound(i);
            bool uj = is_upper_bound(j);
            update_alpha_status(i);
            update_alpha_status(j);
            int k;
            if(ui != is_upper_bound(i))
            {
                Q_i = Q.get_Q(i,l);
                if(ui)
                    for(k=0;k<l;k++)
                        G_bar[k] -= C_i * Q_i[k];
                else
                    for(k=0;k<l;k++)
                        G_bar[k] += C_i * Q_i[k];
            }

            if(uj != is_upper_bound(j))
            {
                Q_j = Q.get_Q(j,l);
                if(uj)
                    for(k=0;k<l;k++)
                        G_bar[k] -= C_j * Q_j[k];
                else
                    for(k=0;k<l;k++)
                        G_bar[k] += C_j * Q_j[k];
            }
        }
    }

    if(iter >= max_iter)
    {
        if(active_size < l)
        {
            // reconstruct the whole gradient to calculate objective value
            reconstruct_gradient();
            active_size = l;
            info("*");
        }
        fprintf(stderr,"\nWARNING: reaching max number of iterations\n");
    }

    // calculate rho

    si->rho = calculate_rho();

    // calculate objective value
    {
        double v = 0;
        int i;
        for(i=0;i<l;i++)
            v += alpha[i] * (G[i] + p[i]);

        si->obj = v/2;
    }

    // put back the solution
    {
        for(int i=0;i<l;i++)
            alpha_[active_set[i]] = alpha[i];
    }

    // juggle everything back
    /*{
        for(int i=0;i<l;i++)
            while(active_set[i] != i)
                swap_index(i,active_set[i]);
                // or Q.swap_index(i,active_set[i]);
    }*/

    si->upper_bound_p = Cp;
    si->upper_bound_n = Cn;

    info("\noptimization finished, #iter = %d\n",iter);

    delete[] p;
    delete[] y;
    delete[] alpha;
    delete[] alpha_status;
    delete[] active_set;
    delete[] G;
    delete[] G_bar;
}

swap_index()

void Solver::swap_index ( int  i, int  j  )

protected

Definition at line 456 of file svm.cpp.

{
    Q->swap_index(i,j);
    swap(y[i],y[j]);
    swap(G[i],G[j]);
    swap(alpha_status[i],alpha_status[j]);
    swap(alpha[i],alpha[j]);
    swap(p[i],p[j]);
    swap(active_set[i],active_set[j]);
    swap(G_bar[i],G_bar[j]);
}

update_alpha_status()

void Solver::update_alpha_status ( int  i )

inlineprotected

Definition at line 436 of file svm.cpp.

    {
        if(alpha[i] >= get_C(i))
            alpha_status[i] = UPPER_BOUND;
        else if(alpha[i] <= 0)
            alpha_status[i] = LOWER_BOUND;
        else alpha_status[i] = FREE;
    }

Member Data Documentation

active_set

int* Solver::active_set

protected

Definition at line 427 of file svm.cpp.

active_size

int Solver::active_size

protected

Definition at line 416 of file svm.cpp.

alpha

double* Solver::alpha

protected

Definition at line 421 of file svm.cpp.

alpha_status

char* Solver::alpha_status

protected

Definition at line 420 of file svm.cpp.

Cn

double Solver::Cn

protected

Definition at line 425 of file svm.cpp.

Cp

double Solver::Cp

protected

Definition at line 425 of file svm.cpp.

eps

double Solver::eps

protected

Definition at line 424 of file svm.cpp.

G

double* Solver::G

protected

Definition at line 418 of file svm.cpp.

G_bar

double* Solver::G_bar

protected

Definition at line 428 of file svm.cpp.

l

int Solver::l

protected

Definition at line 429 of file svm.cpp.

p

double* Solver::p

protected

Definition at line 426 of file svm.cpp.

Q

const QMatrix* Solver::Q

protected

Definition at line 422 of file svm.cpp.

QD

const double* Solver::QD

protected

Definition at line 423 of file svm.cpp.

unshrink

bool Solver::unshrink

protected

Definition at line 430 of file svm.cpp.

y

schar* Solver::y

protected

Definition at line 417 of file svm.cpp.

---

The documentation for this class was generated from the following file:

• libsvm/svm.cpp