gensvm_kernel.c

Go to the documentation of this file.

#include "gensvm_kernel.h"
#include "gensvm_print.h"

void gensvm_kernel_copy_kernelparam_to_data(struct GenModel *model, 
        struct GenData *data)
{
    data->kerneltype = model->kerneltype;
    data->gamma = model->gamma;
    data->coef = model->coef;
    data->degree = model->degree;
}

void gensvm_kernel_preprocess(struct GenModel *model, struct GenData *data)
{
    if (model->kerneltype == K_LINEAR) {
        data->r = data->m;
        return;
    }

    long r, n = data->n;
    double *P = NULL,
           *Sigma = NULL,
           *K = NULL;

    // build the kernel matrix
    K = Calloc(double, n*n);
    gensvm_kernel_compute(model, data, K);

    // generate the eigen decomposition
    r = gensvm_kernel_eigendecomp(K, n, model->kernel_eigen_cutoff, &P, 
            &Sigma);

    // build M and set to data (leave RAW intact)
    gensvm_kernel_trainfactor(data, P, Sigma, r);

    // Set Sigma to data->Sigma (need it again for prediction)
    if (data->Sigma != NULL) {
        free(data->Sigma);
        data->Sigma = NULL;
    }
    data->Sigma = Sigma;

    // write kernel params to data
    gensvm_kernel_copy_kernelparam_to_data(model, data);

    free(K);
    free(P);
}

void gensvm_kernel_postprocess(struct GenModel *model,
        struct GenData *traindata, struct GenData *testdata)
{
    if (model->kerneltype == K_LINEAR) {
        testdata->r = testdata->m;
        return;
    }

    // build the cross kernel matrix between train and test
    double *K2 = gensvm_kernel_cross(model, traindata, testdata);

    // generate the data matrix N = K2 * M * Sigma^{-2}
    gensvm_kernel_testfactor(testdata, traindata, K2);

    free(K2);
}

void gensvm_kernel_compute(struct GenModel *model, struct GenData *data,
        double *K)
{
    long i, j;
    long n = data->n;
    double value;
    double *x1 = NULL,
           *x2 = NULL;

    for (i=0; i<n; i++) {
        for (j=i; j<n; j++) {
            x1 = &data->RAW[i*(data->m+1)+1];
            x2 = &data->RAW[j*(data->m+1)+1];
            if (model->kerneltype == K_POLY)
                value = gensvm_kernel_dot_poly(x1, x2, data->m,
                        model->gamma, model->coef, 
                        model->degree);
            else if (model->kerneltype == K_RBF)
                value = gensvm_kernel_dot_rbf(x1, x2, data->m,
                        model-> gamma);
            else if (model->kerneltype == K_SIGMOID)
                value = gensvm_kernel_dot_sigmoid(x1, x2, 
                        data->m, model->gamma, 
                        model->coef);
            else {
                // LCOV_EXCL_START
                err("[GenSVM Error]: Unknown kernel type in "
                        "gensvm_make_kernel\n");
                exit(EXIT_FAILURE);
                // LCOV_EXCL_STOP
            }
            matrix_set(K, n, i, j, value);
            matrix_set(K, n, j, i, value);
        }
    }
}

long gensvm_kernel_eigendecomp(double *K, long n, double cutoff, double **P_ret,
        double **Sigma_ret)
{
    int M, status, LWORK, *IWORK = NULL,
        *IFAIL = NULL;
    long i, j, num_eigen, cutoff_idx;
    double max_eigen, abstol, *WORK = NULL,
           *Sigma = NULL,
           *P = NULL;

    double *tempSigma = Malloc(double, n);
    double *tempP = Malloc(double, n*n);

    IWORK = Malloc(int, 5*n);
    IFAIL = Malloc(int, n);

    // highest precision eigenvalues, may reduce for speed
    abstol = 2.0*dlamch('S');

    // first perform a workspace query to determine optimal size of the
    // WORK array.
    WORK = Malloc(double, 1);
    status = dsyevx('V', 'A', 'U', n, K, n, 0, 0, 0, 0, abstol, &M,
            tempSigma, tempP, n, WORK, -1, IWORK, IFAIL);
    LWORK = WORK[0];

    // allocate the requested memory for the eigendecomposition
    WORK = (double *)realloc(WORK, LWORK*sizeof(double));
    status = dsyevx('V', 'A', 'U', n, K, n, 0, 0, 0, 0, abstol, &M,
            tempSigma, tempP, n, WORK, LWORK, IWORK, IFAIL);

    if (status != 0) {
        // LCOV_EXCL_START
        err("[GenSVM Error]: Nonzero exit status from dsyevx.\n");
        exit(EXIT_FAILURE);
        // LCOV_EXCL_STOP
    }

    // Select the desired number of eigenvalues, depending on their size.
    // dsyevx sorts eigenvalues in ascending order.
    max_eigen = tempSigma[n-1];
    cutoff_idx = 0;

    for (i=0; i<n; i++) {
        if (tempSigma[i]/max_eigen > cutoff) {
            cutoff_idx = i;
            break;
        }
    }

    num_eigen = n - cutoff_idx;

    Sigma = Calloc(double, num_eigen);
    for (i=0; i<num_eigen; i++) {
        Sigma[i] = tempSigma[n-1 - i];
    }

    // revert P to row-major order and copy only the the columns
    // corresponding to the selected eigenvalues
    P = Calloc(double, n*num_eigen);
    for (j=n-1; j>n-1-num_eigen; j--) {
        for (i=0; i<n; i++) {
            P[i*num_eigen + (n-1)-j] = tempP[i + j*n];
        }
    }

    free(tempSigma);
    free(tempP);
    free(IWORK);
    free(IFAIL);
    free(WORK);

    *Sigma_ret = Sigma;
    *P_ret = P;

    return num_eigen;
}

double *gensvm_kernel_cross(struct GenModel *model, struct GenData *data_train,
        struct GenData *data_test)
{
    long i, j;
    long n_train = data_train->n;
    long n_test = data_test->n;
    long m = data_test->m;
    double value, *x1 = NULL,
           *x2 = NULL,
           *K2 = Calloc(double, n_test * n_train);

    for (i=0; i<n_test; i++) {
        for (j=0; j<n_train; j++) {
            x1 = &data_test->RAW[i*(m+1)+1];
            x2 = &data_train->RAW[j*(m+1)+1];
            if (model->kerneltype == K_POLY)
                value = gensvm_kernel_dot_poly(x1, x2, m, 
                        model->gamma, model->coef, 
                        model->degree);
            else if (model->kerneltype == K_RBF)
                value = gensvm_kernel_dot_rbf(x1, x2, m, 
                        model->gamma);
            else if (model->kerneltype == K_SIGMOID)
                value = gensvm_kernel_dot_sigmoid(x1, x2, m,
                        model->gamma, model->coef);
            else {
                // LCOV_EXCL_START
                err("[GenSVM Error]: Unknown kernel type in "
                        "gensvm_make_crosskernel\n");
                exit(EXIT_FAILURE);
                // LCOV_EXCL_STOP
            }
            matrix_set(K2, n_train, i, j, value);
        }
    }
    return K2;
}

void gensvm_kernel_trainfactor(struct GenData *data, double *P, double *Sigma,
        long r)
{
    long i, j, n = data->n;
    double value;

    // allocate Z
    data->Z = Calloc(double, n*(r+1));

    // Write data->Z = [1 M] = [1 P*Sigma]
    for (i=0; i<n; i++) {
        for (j=0; j<r; j++) {
            value = matrix_get(P, r, i, j);
            value *= matrix_get(Sigma, 1, j, 0);
            matrix_set(data->Z, r+1, i, j+1, value);
        }
        matrix_set(data->Z, r+1, i, 0, 1.0);
    }

    // Set data->r to r so data knows the width of Z
    data->r = r;
}

void gensvm_kernel_testfactor(struct GenData *testdata,
        struct GenData *traindata, double *K2)
{
    long n1, n2, r, i, j;
    double value, *N = NULL,
           *M = NULL;

    n1 = traindata->n;
    n2 = testdata->n;
    r = traindata->r;

    N = Calloc(double, n2*r);
    M = Calloc(double, n1*r);

    // copy M from traindata->Z because we need it in dgemm without column
    // of 1's.
    for (i=0; i<n1; i++) {
        for (j=0; j<r; j++) {
            value = matrix_get(traindata->Z, r+1, i, j+1);
            matrix_set(M, r, i, j, value);
        }
    }

    // Multiply K2 with M and store in N
    cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, n2, r, n1, 1.0,
            K2, n1, M, r, 0.0, N, r);

    // Multiply N with Sigma^{-2}
    for (j=0; j<r; j++) {
        value = pow(matrix_get(traindata->Sigma, 1, j, 0), -2.0);
        for (i=0; i<n2; i++)
            matrix_mul(N, r, i, j, value);
    }

    // write N to Z with a column of ones
    testdata->Z = Calloc(double, n2*(r+1));
    for (i=0; i<n2; i++) {
        for (j=0; j<r; j++) {
            matrix_set(testdata->Z, r+1, i, j+1,
                    matrix_get(N, r, i, j));
        }
        matrix_set(testdata->Z, r+1, i, 0, 1.0);
    }
    // Set r to testdata
    testdata->r = r;

    free(M);
    free(N);
}

double gensvm_kernel_dot_rbf(double *x1, double *x2, long n, double gamma)
{
    long i;
    double value = 0.0;

    for (i=0; i<n; i++)
        value += (x1[i] - x2[i]) * (x1[i] - x2[i]);
    value *= -gamma;
    return exp(value);
}

double gensvm_kernel_dot_poly(double *x1, double *x2, long n, double gamma,
        double coef, double degree)
{
    double value = cblas_ddot(n, x1, 1, x2, 1);
    value *= gamma;
    value += coef;
    return pow(value, degree);
}

double gensvm_kernel_dot_sigmoid(double *x1, double *x2, long n, double gamma,
        double coef)
{
    double value = cblas_ddot(n, x1, 1, x2, 1);
    value *= gamma;
    value += coef;
    return tanh(value);
}

int dsyevx(char JOBZ, char RANGE, char UPLO, int N, double *A, int LDA,
        double VL, double VU, int IL, int IU, double ABSTOL, int *M,
        double *W, double *Z, int LDZ, double *WORK, int LWORK,
        int *IWORK, int *IFAIL)
{
    extern void dsyevx_(char *JOBZ, char *RANGE, char *UPLO, int *Np,
            double *A, int *LDAp, double *VLp, double *VUp,
            int *ILp, int *IUp, double *ABSTOLp, int *M,
            double *W, double *Z, int *LDZp, double *WORK,
            int *LWORKp, int *IWORK, int *IFAIL, int *INFOp);
    int INFO;
    dsyevx_(&JOBZ, &RANGE, &UPLO, &N, A, &LDA, &VL, &VU, &IL, &IU, &ABSTOL,
            M, W, Z, &LDZ, WORK, &LWORK, IWORK, IFAIL, &INFO);
    return INFO;
}

double dlamch(char CMACH)
{
    extern double dlamch_(char *CMACH);
    return dlamch_(&CMACH);
}