gensvm_optimize.c

Go to the documentation of this file.

#include "gensvm_optimize.h"

#ifndef GENSVM_PRINT_ITER
  #define GENSVM_PRINT_ITER 100
#endif

void gensvm_optimize(struct GenModel *model, struct GenData *data)
{
    long it = 0;
    double L, Lbar;

    long n = model->n;
    long m = model->m;
    long K = model->K;

    // initialize the workspace
    struct GenWork *work = gensvm_init_work(model);

    // print some info on the dataset and model configuration
    note("Starting main loop.\n");
    note("Dataset:\n");
    note("\tn = %i\n", n);
    note("\tm = %i\n", m);
    note("\tK = %i\n", K);
    note("Parameters:\n");
    note("\tkappa = %f\n", model->kappa);
    note("\tp = %f\n", model->p);
    note("\tlambda = %15.16f\n", model->lambda);
    note("\tepsilon = %g\n", model->epsilon);
    note("\n");

    // compute necessary simplex vectors
    gensvm_simplex(model);
    gensvm_simplex_diff(model);

    // get initial loss
    L = gensvm_get_loss(model, data, work);
    Lbar = L + 2.0*model->epsilon*L;

    // run main loop
    while ((it < model->max_iter) && (Lbar - L)/L > model->epsilon)
    {
        // ensures V contains newest V and Vbar contains V from
        // previous
        gensvm_get_update(model, data, work);
        if (it > 50)
            gensvm_step_doubling(model);

        Lbar = L;
        L = gensvm_get_loss(model, data, work);

        if (it % GENSVM_PRINT_ITER == 0)
            note("iter = %li, L = %15.16f, Lbar = %15.16f, "
                 "reldiff = %15.16f\n", it, L, Lbar, (Lbar - L)/L);
        it++;
    }

    // status == 0 means training was successful
    model->status = 0;

    // print warnings if necessary
    if (L > Lbar) {
        err("[GenSVM Warning]: Negative step occurred in "
                "majorization.\n");
        model->status = 1;
    }

    if (it >= model->max_iter) {
        err("[GenSVM Warning]: maximum number of iterations "
                "reached.\n");
        model->status = 2;
    }

    // print final iteration count and loss
    note("Optimization finished, iter = %li, loss = %15.16f, "
            "rel. diff. = %15.16f\n", it-1, L,
            (Lbar - L)/L);

    // compute and print the number of SVs in the model
    note("Number of support vectors: %li\n", gensvm_num_sv(model));

    // store the training error in the model
    model->training_error = (Lbar - L)/L;

    // store the iteration count in the model
    model->elapsed_iter = it - 1;

    // free the workspace
    gensvm_free_work(work);
}

double gensvm_get_loss(struct GenModel *model, struct GenData *data, 
        struct GenWork *work)
{
    long i, j;
    long n = model->n;
    long K = model->K;
    long m = model->m;

    double value, rowvalue, loss = 0.0;

    gensvm_calculate_errors(model, data, work->ZV);
    gensvm_calculate_huber(model);

    for (i=0; i<n; i++) {
        rowvalue = 0;
        value = 0;
        for (j=0; j<K; j++) {
            if (j == (data->y[i]-1))
                continue;
            value = matrix_get(model->H, K, i, j);
            value = pow(value, model->p);
            rowvalue += value;
        }
        rowvalue = pow(rowvalue, 1.0/(model->p));
        rowvalue *= model->rho[i];
        loss += rowvalue;
    }
    loss /= ((double) n);

    value = 0;
    for (i=1; i<m+1; i++) {
        for (j=0; j<K-1; j++) {
            value += pow(matrix_get(model->V, K-1, i, j), 2.0);
        }
    }
    loss += model->lambda * value;

    return loss;
}

void gensvm_step_doubling(struct GenModel *model)
{
    long i, j;
    double value;

    long m = model->m;
    long K = model->K;

    for (i=0; i<m+1; i++) {
        for (j=0; j<K-1; j++) {
            matrix_mul(model->V, K-1, i, j, 2.0);
            value = - matrix_get(model->Vbar, K-1, i, j);
            matrix_add(model->V, K-1, i, j, value);
        }
    }
}

void gensvm_calculate_huber(struct GenModel *model)
{
    long i, j;
    double q, value;

    for (i=0; i<model->n; i++) {
        for (j=0; j<model->K; j++) {
            q = matrix_get(model->Q, model->K, i, j);
            value = 0.0;
            if (q <= -model->kappa) {
                value = 1.0 - q - (model->kappa+1.0)/2.0;
            } else if (q <= 1.0) {
                value = 1.0/(2.0*model->kappa+2.0)*pow(1.0 - q,
                            2.0);
            }
            matrix_set(model->H, model->K, i, j, value);
        }
    }
}

void gensvm_calculate_errors(struct GenModel *model, struct GenData *data,
        double *ZV)
{
    long i, j;
    double q, *uu_row = NULL;

    long n = model->n;
    long K = model->K;

    gensvm_calculate_ZV(model, data, ZV);

    for (i=0; i<n; i++) {
        for (j=0; j<K; j++) {
            if (j == (data->y[i]-1))
                continue;
            uu_row = &model->UU[((data->y[i]-1)*K+j)*(K-1)];
            q = cblas_ddot(K-1, &ZV[i*(K-1)], 1, uu_row, 1);
            matrix_set(model->Q, K, i, j, q);
        }
    }
}