gensvm_update.c

Go to the documentation of this file.

#include "gensvm_update.h"

#ifndef GENSVM_BLOCK_SIZE
  #define GENSVM_BLOCK_SIZE 512
#endif

double gensvm_calculate_omega(struct GenModel *model, struct GenData *data,
        long i)
{
    long j;
    double h, omega = 0.0,
           p = model->p;

    for (j=0; j<model->K; j++) {
        if (j == (data->y[i]-1))
            continue;
        h = matrix_get(model->H, model->K, i, j);
        omega += pow(h, p);
    }
    omega = (1.0/p)*pow(omega, 1.0/p - 1.0);

    return omega;
}

bool gensvm_majorize_is_simple(struct GenModel *model, struct GenData *data,
        long i)
{
    long j;
    double h, value = 0;
    for (j=0; j<model->K; j++) {
        if (j == (data->y[i]-1))
            continue;
        h = matrix_get(model->H, model->K, i, j);
        value += h > 0;
        if (value > 1)
            return false;
    }
    return true;
}

void gensvm_calculate_ab_non_simple(struct GenModel *model, long i, long j,
        double *a, double *b_aq)
{
    double q = matrix_get(model->Q, model->K, i, j);
    double p = model->p;
    double kappa = model->kappa;
    const double a2g2 = 0.25*p*(2.0*p - 1.0)*pow((kappa+1.0)/2.0,p-2.0);

    if (2.0 - model->p < 1e-2) {
        if (q <= - kappa) {
            *b_aq = 0.5 - kappa/2.0 - q;
        } else if ( q <= 1.0) {
            *b_aq = pow(1.0 - q, 3.0)/(2.0*pow(kappa + 1.0, 2.0));
        } else {
            *b_aq = 0;
        }
        *a = 1.5;
    } else {
        if (q <= (p + kappa - 1.0)/(p - 2.0)) {
            *a = 0.25*pow(p, 2.0)*pow(0.5 - kappa/2.0 - q, p - 2.0);
        } else if (q <= 1.0) {
            *a = a2g2;
        } else {
            *a = 0.25*pow(p, 2.0)*pow((p/(p - 2.0))*(0.5 -
                        kappa/2.0 - q), p - 2.0);
            *b_aq = (*a)*(2.0*q + kappa - 1.0)/(p - 2.0) +
                0.5*p*pow(p/(p - 2.0)*(0.5 - kappa/2.0 - q),
                        p - 1.0);
        }
        if (q <= -kappa) {
            *b_aq = 0.5*p*pow(0.5 - kappa/2.0 - q, p - 1.0);
        } else if ( q <= 1.0) {
            *b_aq = p*pow(1.0 - q, 2.0*p - 1.0)/pow(2*kappa+2.0, p);
        }
    }
}

void gensvm_calculate_ab_simple(struct GenModel *model, long i, long j,
        double *a, double *b_aq)
{
    double q = matrix_get(model->Q, model->K, i, j);

    if (q <= - model->kappa) {
        *a = 0.25/(0.5 - model->kappa/2.0 - q);
        *b_aq = 0.5;
    } else if (q <= 1.0) {
        *a = 1.0/(2.0*model->kappa + 2.0);
        *b_aq = (1.0 - q)*(*a);
    } else {
        *a = -0.25/(0.5 - model->kappa/2.0 - q);
        *b_aq = 0;
    }
}

double gensvm_get_alpha_beta(struct GenModel *model, struct GenData *data,
        long i, double *beta)
{
    bool simple;
    long j, K = model->K;
    double omega, a, b_aq = 0.0,
           alpha = 0.0;
    double *uu_row = NULL;
    const double in = 1.0/((double) model->n);

    simple = gensvm_majorize_is_simple(model, data, i);
    omega = simple ? 1.0 : gensvm_calculate_omega(model, data, i);

    Memset(beta, double, K-1);
    for (j=0; j<K; j++) {
        // skip the class y_i = k
        if (j == (data->y[i]-1))
            continue;

        // calculate the a_ijk and (b_ijk - a_ijk q_i^(kj)) values
        if (simple) {
            gensvm_calculate_ab_simple(model, i, j, &a, &b_aq);
        } else {
            gensvm_calculate_ab_non_simple(model, i, j, &a, &b_aq);
        }

        // daxpy on beta and UU
        // daxpy does: y = a*x + y
        // so y = beta, UU_row = x, a = factor
        b_aq *= model->rho[i] * omega * in;
        uu_row = &model->UU[((data->y[i]-1)*K+j)*(K-1)];
        cblas_daxpy(K-1, b_aq, uu_row, 1, beta, 1);

        // increment Avalue
        alpha += a;
    }
    alpha *= omega * model->rho[i] * in;
    return alpha;
}

void gensvm_get_update(struct GenModel *model, struct GenData *data,
        struct GenWork *work)
{
    int status;
    long i, j;

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

    // compute the ZAZ and ZB matrices
    gensvm_get_ZAZ_ZB(model, data, work);

    // Calculate right-hand side of system we want to solve
    // dsymm performs ZB := 1.0 * (ZAZ) * Vbar + 1.0 * ZB
    // the right-hand side is thus stored in ZB after this call
    // Note: LDB and LDC are second dimensions of the matrices due to
    // Row-Major order
    cblas_dsymm(CblasRowMajor, CblasLeft, CblasUpper, m+1, K-1, 1,
            work->ZAZ, m+1, model->V, K-1, 1.0, work->ZB, K-1);

    // Calculate left-hand side of system we want to solve
    // Add lambda to all diagonal elements except the first one. Recall
    // that ZAZ is of size m+1 and is symmetric.
    for (i=m+2; i<=m*(m+2); i+=m+2)
        work->ZAZ[i] += model->lambda;

    // Lapack uses column-major order, so we transform the ZB matrix to
    // correspond to this.
    for (i=0; i<m+1; i++)
        for (j=0; j<K-1; j++)
            work->ZBc[j*(m+1)+i] = work->ZB[i*(K-1)+j];

    // Solve the system using dposv. Note that above the upper triangular
    // part has always been used in row-major order for ZAZ. This
    // corresponds to the lower triangular part in column-major order.
    status = dposv('L', m+1, K-1, work->ZAZ, m+1, work->ZBc, m+1);

    // Use dsysv as fallback, for when the ZAZ matrix is not positive
    // semi-definite for some reason (perhaps due to rounding errors).
    // This step shouldn't be necessary but is included for safety.
    if (status != 0) {
        err("[GenSVM Warning]: Received nonzero status from "
                "dposv: %i\n", status);
        int *IPIV = Malloc(int, m+1);
        double *WORK = Malloc(double, 1);
        status = dsysv('L', m+1, K-1, work->ZAZ, m+1, IPIV, work->ZBc,
                m+1, WORK, -1);

        int LWORK = WORK[0];
        WORK = Realloc(WORK, double, LWORK);
        status = dsysv('L', m+1, K-1, work->ZAZ, m+1, IPIV, work->ZBc,
                m+1, WORK, LWORK);
        if (status != 0)
            err("[GenSVM Warning]: Received nonzero "
                    "status from dsysv: %i\n", status);
        free(WORK);
        WORK = NULL;
        free(IPIV);
        IPIV = NULL;
    }

    // the solution is now stored in ZBc, in column-major order. Here we
    // convert this back to row-major order
    for (i=0; i<m+1; i++)
        for (j=0; j<K-1; j++)
            work->ZB[i*(K-1)+j] = work->ZBc[j*(m+1)+i];

    // copy the old V to Vbar and the new solution to V
    for (i=0; i<m+1; i++) {
        for (j=0; j<K-1; j++) {
            matrix_set(model->Vbar, K-1, i, j,
                    matrix_get(model->V, K-1, i, j));
            matrix_set(model->V, K-1, i, j,
                    matrix_get(work->ZB, K-1, i, j));
        }
    }
}

void gensvm_get_ZAZ_ZB_dense(struct GenModel *model, struct GenData *data,
        struct GenWork *work)
{
    long i;
    double alpha, sqalpha;

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

    // generate Z'*A*Z and Z'*B by rank 1 operations
    for (i=0; i<n; i++) {
        alpha = gensvm_get_alpha_beta(model, data, i, work->beta);

        // calculate row of matrix LZ, which is a scalar
        // multiplication of sqrt(alpha_i) and row z_i' of Z
        // Note that we use the fact that the first column of Z is
        // always 1, by only computing the product for m values and
        // copying the first element over.
        sqalpha = sqrt(alpha);
        work->LZ[i*(m+1)] = sqalpha;
        cblas_daxpy(m, sqalpha, &data->Z[i*(m+1)+1], 1,
                &work->LZ[i*(m+1)+1], 1);

        // rank 1 update of matrix Z'*B
        // Note: LDA is the second dimension of ZB because of
        // Row-Major order
        cblas_dger(CblasRowMajor, m+1, K-1, 1, &data->Z[i*(m+1)], 1,
                work->beta, 1, work->ZB, K-1);
    }

    // calculate Z'*A*Z by symmetric multiplication of LZ with itself
    // (ZAZ = (LZ)' * (LZ)
    cblas_dsyrk(CblasRowMajor, CblasUpper, CblasTrans, m+1, n, 1.0,
            work->LZ, m+1, 0.0, work->ZAZ, m+1);
}

void gensvm_get_ZAZ_ZB_sparse(struct GenModel *model, struct GenData *data,
        struct GenWork *work)
{
    long *Zia = NULL,
         *Zja = NULL;
    long b, i, j, k, K, jj, kk, jj_start, jj_end, blk_start, blk_end,
         rem_size, n_blocks, n_row = data->spZ->n_row,
         n_col = data->spZ->n_col;
    double temp, alpha, z_ij, *vals = NULL;

    K = model->K;
    Zia = data->spZ->ia;
    Zja = data->spZ->ja;
    vals = data->spZ->values;

    // calculate ZAZ using blocks of rows of Z. This helps avoiding 
    // rounding errors, which increases precision, and in turn helps 
    // convergence of the IM algorithm.
    // see also: http://stackoverflow.com/q/40286989/
    n_blocks = floor(n_row / GENSVM_BLOCK_SIZE);
    rem_size = n_row % GENSVM_BLOCK_SIZE;

    for (b=0; b<=n_blocks; b++) {
        blk_start = b * GENSVM_BLOCK_SIZE;
        blk_end = blk_start;
        blk_end += (b == n_blocks) ? rem_size : GENSVM_BLOCK_SIZE;

        Memset(work->tmpZAZ, double, n_col*n_col);
        for (i=blk_start; i<blk_end; i++) {
            alpha = gensvm_get_alpha_beta(model, data, i, 
                    work->beta);
            jj_start = Zia[i];
            jj_end = Zia[i+1];

            for (jj=jj_start; jj<jj_end; jj++) {
                j = Zja[jj];
                z_ij = vals[jj];
                cblas_daxpy(K-1, z_ij, work->beta, 1,
                        &work->ZB[j*(K-1)], 1);
                z_ij *= alpha;
                for (kk=jj; kk<jj_end; kk++) {
                    matrix_add(work->tmpZAZ, n_col, j, 
                            Zja[kk], 
                            z_ij * vals[kk]);
                }
            }
        }

        // copy the intermediate results over to the actual ZAZ matrix
        for (j=0; j<n_col; j++) {
            for (k=j; k<n_col; k++) {
                temp = matrix_get(work->tmpZAZ, n_col, j, k);
                matrix_add(work->ZAZ, n_col, j, k, temp);
            }
        }
    }
}

void gensvm_get_ZAZ_ZB(struct GenModel *model, struct GenData *data,
        struct GenWork *work)
{
    gensvm_reset_work(work);

    if (data->Z == NULL)
        gensvm_get_ZAZ_ZB_sparse(model, data, work);
    else
        gensvm_get_ZAZ_ZB_dense(model, data, work);
}

int dposv(char UPLO, int N, int NRHS, double *A, int LDA, double *B,
        int LDB)
{
    extern void dposv_(char *UPLO, int *Np, int *NRHSp, double *A,
            int *LDAp, double *B, int *LDBp, int *INFOp);
    int INFO;
    dposv_(&UPLO, &N, &NRHS, A, &LDA, B, &LDB, &INFO);
    return INFO;
}

int dsysv(char UPLO, int N, int NRHS, double *A, int LDA, int *IPIV,
        double *B, int LDB, double *WORK, int LWORK)
{
    extern void dsysv_(char *UPLO, int *Np, int *NRHSp, double *A,
            int *LDAp, int *IPIV, double *B, int *LDBp,
            double *WORK, int *LWORK, int *INFOp);
    int INFO;
    dsysv_(&UPLO, &N, &NRHS, A, &LDA, IPIV, B, &LDB, WORK, &LWORK, &INFO);
    return INFO;
}