svm.c

Go to the documentation of this file.

/*
Copyright (C) 2013 Milan Sulc.
Copyright (C) 2012 Daniele Perrone.
Copyright (C) 2011-13 Andrea Vedaldi.

All rights reserved.

This file is part of the VLFeat library and is made available under
the terms of the BSD license (see the COPYING file).
*/

/*

<!-- ------------------------------------------------------------ --->
@section svm-pegasos PEGASOS
<!-- ------------------------------------------------------------ --->

<!-- ------------------------------------------------------------ --->
@subsection svm-pegasos-algorithm Algorithm
<!-- ------------------------------------------------------------ --->

PEGASOS @cite{shalev-shwartz07pegasos} is a stochastic subgradient
optimizer. At the <em>t</em>-th iteration the algorithm:

- Samples uniformly at random as subset @f$ A_t @f$ of <em>k</em> of
training pairs @f$(x,y)@f$ from the <em>m</em> pairs provided for
training (this subset is called mini batch).
- Computes a subgradient @f$ \nabla_t @f$ of the function @f$ E_t(w) =
\frac{1}{2}\|w\|^2 + \frac{1}{k} \sum_{(x,y) \in A_t} \ell(w;(x,y))
@f$ (this is the SVM objective function restricted to the
minibatch).
- Compute an intermediate weight vector @f$ w_{t+1/2} @f$ by doing a
step @f$ w_{t+1/2} = w_t - \alpha_t \nabla_t @f$ with learning rate
@f$ \alpha_t = 1/(\eta t) @f$ along the subgradient. Note that the
learning rate is inversely proportional to the iteration number.
- Back projects the weight vector @f$ w_{t+1/2} @f$ on the
hypersphere of radius @f$ \sqrt{\lambda} @f$ to obtain the next
model estimate @f$ w_{t+1} @f$:
@f[
w_t = \min\{1, \sqrt{\lambda}/\|w\|\} w_{t+1/2}.
@f]
The hypersphere is guaranteed to contain the optimal weight vector
@f$ w^* @f$.

VLFeat implementation fixes to one the size of the mini batches @f$ k
@f$.


<!-- ------------------------------------------------------------ --->
@subsection svm-pegasos-permutation Permutation
<!-- ------------------------------------------------------------ --->

VLFeat PEGASOS can use a user-defined permutation to decide the order
in which data points are visited (instead of using random
sampling). By specifying a permutation the algorithm is guaranteed to
visit each data point exactly once in each loop. The permutation needs
not to be bijective. This can be used to visit certain data samples
more or less often than others, implicitly reweighting their relative
importance in the SVM objective function. This can be used to balance
the data.

<!-- ------------------------------------------------------------ --->
@subsection svm-pegasos-kernels Non-linear kernels
<!-- ------------------------------------------------------------ --->

PEGASOS can be extended to non-linear kernels, but the algorithm is
not particularly efficient in this setting [1]. When possible, it may
be preferable to work with explicit feature maps.

Let @f$ k(x,y) @f$ be a positive definite kernel. A <em>feature
map</em> is a function @f$ \Psi(x) @f$ such that @f$ k(x,y) = \langle
\Psi(x), \Psi(y) \rangle @f$. Using this representation the non-linear
SVM learning objective function writes:

@f[
\min_{w} \frac{\lambda}{2} \|w\|^2 + \frac{1}{m} \sum_{i=1}^n
\ell(w; (\Psi(x)_i,y_i)).
@f]

Thus the only difference with the linear case is that the feature @f$
\Psi(x) @f$ is used in place of the data @f$ x @f$.

@f$ \Psi(x) @f$ can be learned off-line, for instance by using the
incomplete Cholesky decomposition @f$ V^\top V @f$ of the Gram matrix
@f$ K = [k(x_i,x_j)] @f$ (in this case @f$ \Psi(x_i) @f$ is the
<em>i</em>-th columns of <em>V</em>). Alternatively, for additive
kernels (e.g. intersection, Chi2) the explicit feature map computed by
@ref homkermap.h can be used.

For additive kernels it is also possible to perform the feature
expansion online inside the solver, setting the specific feature map
via ::vl_svmdataset_set_map. This is particular useful to keep the
size of the training data small, when the number of the samples is big
or the memory is limited.
*/

#include "svm.h"
#include "mathop.h"
#include <string.h>

struct VlSvm_ {
  VlSvmSolverType solver ;      
  vl_size dimension ;           
  double * model ;              
  double bias ;                 
  double biasMultiplier ;       
  /* valid during a run */
  double lambda ;               
  void const * data ;
  vl_size numData ;
  double const * labels ;       
  double const * weights ;      
  VlSvmDataset * ownDataset ;   
  VlSvmDiagnosticFunction diagnosticFn ;
  void * diagnosticFnData ;
  vl_size diagnosticFrequency ; 
  VlSvmLossFunction lossFn ;
  VlSvmLossFunction conjugateLossFn ;
  VlSvmLossFunction lossDerivativeFn ;
  VlSvmDcaUpdateFunction dcaUpdateFn ;
  VlSvmInnerProductFunction innerProductFn ;
  VlSvmAccumulateFunction accumulateFn ;

  vl_size iteration ;           
  vl_size maxNumIterations ;    
  double epsilon ;              
  /* Book keeping */
  VlSvmStatistics statistics ;  
  double * scores ;

  /* SGD specific */
  double  biasLearningRate ;    
  /* SDCA specific */
  double * alpha ;              
} ;

/* ---------------------------------------------------------------- */

VlSvm *
vl_svm_new (VlSvmSolverType type,
            double const * data,
            vl_size dimension,
            vl_size numData,
            double const * labels,
            double lambda)
{
  VlSvmDataset * dataset = vl_svmdataset_new(VL_TYPE_DOUBLE, (void*)data, dimension, numData) ;
  VlSvm * self = vl_svm_new_with_dataset (type, dataset, labels, lambda) ;
  self->ownDataset = dataset ;
  return self ;
}

VlSvm *
vl_svm_new_with_dataset (VlSvmSolverType solver,
                         VlSvmDataset * dataset,
                         double const * labels,
                         double lambda)
{
  VlSvm * self = vl_svm_new_with_abstract_data (solver,
                                             dataset,
                                             vl_svmdataset_get_dimension(dataset),
                                             vl_svmdataset_get_num_data(dataset),
                                             labels,
                                             lambda) ;
  vl_svm_set_data_functions (self,
                             vl_svmdataset_get_inner_product_function(dataset),
                             vl_svmdataset_get_accumulate_function(dataset)) ;
  return self ;
}

VlSvm *
vl_svm_new_with_abstract_data (VlSvmSolverType solver,
                               void * data,
                               vl_size dimension,
                               vl_size numData,
                               double const * labels,
                               double lambda)
{
  VlSvm * self = vl_calloc(1,sizeof(VlSvm)) ;

  assert(dimension >= 1) ;
  assert(numData >= 1) ;
  assert(labels) ;

  self->solver = solver ;

  self->dimension = dimension ;
  self->model = 0 ;
  self->bias = 0 ;
  self->biasMultiplier = 1.0 ;

  self->lambda = lambda ;
  self->data = data ;
  self->numData = numData ;
  self->labels = labels ;

  self->diagnosticFrequency = numData ;
  self->diagnosticFn = 0 ;
  self->diagnosticFnData = 0 ;

  self->lossFn = vl_svm_hinge_loss ;
  self->conjugateLossFn = vl_svm_hinge_conjugate_loss ;
  self->lossDerivativeFn = vl_svm_hinge_loss_derivative ;
  self->dcaUpdateFn = vl_svm_hinge_dca_update ;

  self->innerProductFn = 0 ;
  self->accumulateFn = 0 ;

  self->iteration = 0 ;
  self->maxNumIterations = VL_MAX((double)numData, vl_ceil_f(10.0 / lambda)) ;
  self->epsilon = 1e-2 ;

  /* SGD */
  self->biasLearningRate = 0.01 ;

  /* SDCA */
  self->alpha = 0 ;

  /* allocations */
  self->model = vl_calloc(dimension, sizeof(double)) ;
  if (self->model == NULL) goto err_alloc ;

  if (self->solver == VlSvmSolverSdca) {
    self->alpha = vl_calloc(self->numData, sizeof(double)) ;
    if (self->alpha == NULL) goto err_alloc ;
  }

  self->scores = vl_calloc(numData, sizeof(double)) ;
  if (self->scores == NULL) goto err_alloc ;

  return self ;

err_alloc:
  if (self->scores) {
    vl_free (self->scores) ;
    self->scores = 0 ;
  }
  if (self->model) {
    vl_free (self->model) ;
    self->model = 0 ;
  }
  if (self->alpha) {
    vl_free (self->alpha) ;
    self->alpha = 0 ;
  }
  return 0 ;
}

void
vl_svm_delete (VlSvm * self)
{
  if (self->model) {
    vl_free (self->model) ;
    self->model = 0 ;
  }
  if (self->alpha) {
    vl_free (self->alpha) ;
    self->alpha = 0 ;
  }
  if (self->ownDataset) {
    vl_svmdataset_delete(self->ownDataset) ;
    self->ownDataset = 0 ;
  }
  vl_free (self) ;
}

/* ---------------------------------------------------------------- */
/*                                              Setters and getters */
/* ---------------------------------------------------------------- */

void vl_svm_set_epsilon (VlSvm *self, double epsilon)
{
  assert(self) ;
  assert(epsilon >= 0) ;
  self->epsilon = epsilon ;
}

double vl_svm_get_epsilon (VlSvm const *self)
{
  assert(self) ;
  return self->epsilon ;
}

void vl_svm_set_bias_learning_rate (VlSvm *self, double rate)
{
  assert(self) ;
  assert(rate > 0) ;
  self->biasLearningRate = rate ;
}

double vl_svm_get_bias_learning_rate (VlSvm const *self)
{
  assert(self) ;
  return self->biasLearningRate ;
}

void vl_svm_set_bias_multiplier (VlSvm * self, double b)
{
  assert(self) ;
  assert(b >= 0) ;
  self->biasMultiplier = b ;
}

double vl_svm_get_bias_multiplier (VlSvm const * self)
{
  assert(self) ;
  return self->biasMultiplier ;
}

void vl_svm_set_iteration_number (VlSvm *self, vl_uindex n)
{
  assert(self) ;
  self->iteration = n ;
}

vl_size vl_svm_get_iteration_number (VlSvm const *self)
{
  assert(self) ;
  return self->iteration ;
}

void vl_svm_set_max_num_iterations (VlSvm *self, vl_size n)
{
  assert(self) ;
  self->maxNumIterations = n ;
}

vl_size vl_svm_get_max_num_iterations (VlSvm const *self)
{
  assert(self) ;
  return self->maxNumIterations ;
}

void vl_svm_set_diagnostic_frequency (VlSvm *self, vl_size f)
{
  assert(self) ;
  assert(f > 0) ;
  self->diagnosticFrequency = f ;
}

vl_size vl_svm_get_diagnostic_frequency (VlSvm const *self)
{
  assert(self) ;
  return self->diagnosticFrequency ;
}

VlSvmSolverType vl_svm_get_solver (VlSvm const * self)
{
  assert(self) ;
  return self->solver ;
}

void vl_svm_set_lambda (VlSvm * self, double lambda)
{
  assert(self) ;
  assert(lambda >= 0) ;
  self->lambda = lambda ;
}

double vl_svm_get_lambda (VlSvm const * self)
{
  assert(self) ;
  return self->lambda ;
}

void vl_svm_set_weights (VlSvm * self, double const *weights)
{
  assert(self) ;
  self->weights = weights ;
}

double const *vl_svm_get_weights (VlSvm const * self)
{
  assert(self) ;
  return self->weights ;
}

/* ---------------------------------------------------------------- */
/*                                                         Get data */
/* ---------------------------------------------------------------- */

vl_size vl_svm_get_dimension (VlSvm *self)
{
  assert(self) ;
  return self->dimension ;
}

vl_size vl_svm_get_num_data (VlSvm *self)
{
  assert(self) ;
  return self->numData ;
}

double const * vl_svm_get_model (VlSvm const *self)
{
  assert(self) ;
  return self->model ;
}

void vl_svm_set_model (VlSvm *self, double const *model)
{
  assert(self) ;
  assert(model) ;
  memcpy(self->model, model, sizeof(double) * vl_svm_get_dimension(self)) ;
}

void vl_svm_set_bias (VlSvm *self, double b)
{
  assert(self);
  if (self->biasMultiplier) {
    self->bias = b / self->biasMultiplier ;
  }
}

double vl_svm_get_bias (VlSvm const *self)
{
  assert(self) ;
  return self->bias * self->biasMultiplier ;
}

VlSvmStatistics const * vl_svm_get_statistics (VlSvm const *self)
{
  assert(self) ;
  return &self->statistics ;
}

double const * vl_svm_get_scores (VlSvm const *self)
{
  return self->scores ;
}

/* ---------------------------------------------------------------- */
/*                                                        Callbacks */
/* ---------------------------------------------------------------- */

void
vl_svm_set_diagnostic_function (VlSvm *self, VlSvmDiagnosticFunction f, void *data) {
  self->diagnosticFn = f ;
  self->diagnosticFnData = data ;
}

void vl_svm_set_data_functions (VlSvm *self, VlSvmInnerProductFunction inner, VlSvmAccumulateFunction acc)
{
  assert(self) ;
  assert(inner) ;
  assert(acc) ;
  self->innerProductFn = inner ;
  self->accumulateFn = acc ;
}

void vl_svm_set_loss_function (VlSvm *self, VlSvmLossFunction f)
{
  assert(self) ;
  self->lossFn = f ;
}

void vl_svm_set_loss_derivative_function (VlSvm *self, VlSvmLossFunction f)
{
  assert(self) ;
  self->lossDerivativeFn = f ;
}

void vl_svm_set_conjugate_loss_function (VlSvm *self, VlSvmLossFunction f)
{
  assert(self) ;
  self->conjugateLossFn = f ;
}

void vl_svm_set_dca_update_function (VlSvm *self, VlSvmDcaUpdateFunction f)
{
  assert(self) ;
  self->dcaUpdateFn = f ;
}

void
vl_svm_set_loss (VlSvm *self, VlSvmLossType loss)
{
#define SETLOSS(x,y) \
case VlSvmLoss ## x: \
  vl_svm_set_loss_function(self, vl_svm_ ## y ## _loss) ; \
  vl_svm_set_loss_derivative_function(self, vl_svm_ ## y ## _loss_derivative) ; \
  vl_svm_set_conjugate_loss_function(self, vl_svm_ ## y ## _conjugate_loss) ; \
  vl_svm_set_dca_update_function(self, vl_svm_ ## y ## _dca_update) ; \
  break;

  switch (loss) {
      SETLOSS(Hinge, hinge) ;
      SETLOSS(Hinge2, hinge2) ;
      SETLOSS(L1, l1) ;
      SETLOSS(L2, l2) ;
      SETLOSS(Logistic, logistic) ;
    default:
      assert(0) ;
  }
}

/* ---------------------------------------------------------------- */
/*                                               Pre-defined losses */
/* ---------------------------------------------------------------- */

double
vl_svm_hinge_loss (double inner, double label)
{
  return VL_MAX(1 - label * inner, 0.0);
}

double
vl_svm_hinge_loss_derivative (double inner, double label)
{
  if (label * inner < 1.0) {
    return - label ;
  } else {
    return 0.0 ;
  }
}

double
vl_svm_hinge_conjugate_loss (double u, double label) {
  double z = label * u ;
  if (-1 <= z && z <= 0) {
    return label * u ;
  } else {
    return VL_INFINITY_D ;
  }
}

double
vl_svm_hinge_dca_update (double alpha, double inner, double norm2, double label) {
  double palpha = (label - inner) / norm2 + alpha ;
  return label * VL_MAX(0, VL_MIN(1, label * palpha)) - alpha ;
}

double
vl_svm_hinge2_loss (double inner,double label)
{
  double z = VL_MAX(1 - label * inner, 0.0) ;
  return z*z ;
}

double
vl_svm_hinge2_loss_derivative (double inner, double label)
{
  if (label * inner < 1.0) {
    return 2 * (inner - label) ;
  } else {
    return 0 ;
  }
}

double
vl_svm_hinge2_conjugate_loss (double u, double label) {
  if (label * u <= 0) {
    return (label + u/4) * u ;
  } else {
    return VL_INFINITY_D ;
  }
}

double
vl_svm_hinge2_dca_update (double alpha, double inner, double norm2, double label) {
  double palpha = (label - inner - 0.5*alpha) / (norm2 + 0.5) + alpha ;
  return label * VL_MAX(0, label * palpha) - alpha ;
}

double
vl_svm_l1_loss (double inner,double label)
{
  return vl_abs_d(label - inner) ;
}

double
vl_svm_l1_loss_derivative (double inner, double label)
{
  if (label > inner) {
    return - 1.0 ;
  } else {
    return + 1.0 ;
  }
}

double
vl_svm_l1_conjugate_loss (double u, double label) {
  if (vl_abs_d(u) <= 1) {
    return label*u ;
  } else {
    return VL_INFINITY_D ;
  }
}

double
vl_svm_l1_dca_update (double alpha, double inner, double norm2, double label) {
  if (vl_abs_d(alpha) <= 1) {
    double palpha = (label - inner) / norm2 + alpha ;
    return VL_MAX(-1.0, VL_MIN(1.0, palpha)) - alpha ;
  } else {
    return VL_INFINITY_D ;
  }
}

double
vl_svm_l2_loss (double inner,double label)
{
  double z = label - inner ;
  return z*z ;
}

double
vl_svm_l2_loss_derivative (double inner, double label)
{
  return - 2 * (label - inner) ;
}

double
vl_svm_l2_conjugate_loss (double u, double label) {
  return (label + u/4) * u ;
}

double
vl_svm_l2_dca_update (double alpha, double inner, double norm2, double label) {
  return (label - inner - 0.5*alpha) / (norm2 + 0.5) ;
}

double
vl_svm_logistic_loss (double inner,double label)
{
  double z = label * inner ;
  if (z >= 0) {
    return log(1.0 + exp(-z)) ;
  } else {
    return -z + log(exp(z) + 1.0) ;
  }
}

double
vl_svm_logistic_loss_derivative (double inner, double label)
{
  double z = label * inner ;
  double t = 1 / (1 + exp(-z)) ; /* this is stable for z << 0 too */
  return label * (t - 1) ; /*  = -label exp(-z) / (1 + exp(-z)) */
}

VL_INLINE double xlogx(double x)
{
  if (x <= 1e-10) return 0 ;
  return x*log(x) ;
}

double
vl_svm_logistic_conjugate_loss (double u, double label) {
  double z = label * u ;
  if (-1 <= z && z <= 0) {
    return xlogx(-z) + xlogx(1+z) ;
  } else {
    return VL_INFINITY_D ;
  }
}

double
vl_svm_logistic_dca_update (double alpha, double inner, double norm2, double label) {
  /*
   The goal is to solve the problem

   min_delta A/2 delta^2 + B delta + l*(-alpha - delta|y),  -1 <= - y (alpha+delta) <= 0

   where A = norm2, B = inner, and y = label. To simplify the notation, we set

     f(beta) = beta * log(beta) + (1 - beta) * log(1 - beta)

   where beta = y(alpha + delta) such that

     l*(-alpha - delta |y) = f(beta).

   Hence 0 <= beta <= 1, delta = + y beta - alpha. Substituting

     min_beta A/2 beta^2 + y (B - A alpha) beta + f(beta) + const

   The Newton step is then given by

     beta = beta - (A beta + y(B - A alpha) + df) / (A + ddf).

   However, the function is singluar for beta=0 and beta=1 (infinite
   first and second order derivatives). Since the function is monotonic
   (second derivarive always strictly greater than zero) and smooth,
   we canuse bisection to find the zero crossing of the first derivative.
   Once one is sufficiently close to the optimum, a one or two Newton
   steps are sufficien to land on it with excellent accuracy.
   */

  double  df, ddf, der, dder ;
  vl_index t ;

  /* bisection */
  double beta1 = 0 ;
  double beta2 = 1 ;
  double beta = 0.5 ;

  for (t = 0 ; t < 5 ; ++t) {
    df = log(beta) - log(1-beta) ;
    der = norm2 * beta + label * (inner - norm2*alpha) + df ;
    if (der >= 0) {
      beta2 = beta ;
    } else {
      beta1 = beta ;
    }
    beta = 0.5 * (beta1 + beta2) ;
  }

#if 1
  /* a final Newton step, but not too close to the singularities */
  for (t = 0 ; (t < 2) & (beta > VL_EPSILON_D) & (beta < 1-VL_EPSILON_D) ; ++t) {
    df = log(beta) - log(1-beta) ;
    ddf = 1 / (beta * (1-beta)) ;
    der = norm2 * beta + label * (inner - norm2*alpha) + df ;
    dder = norm2 + ddf ;
    beta -= der / dder ;
    beta = VL_MAX(0, VL_MIN(1, beta)) ;
  }
#endif

  return label * beta - alpha ;
}

/* ---------------------------------------------------------------- */

void _vl_svm_update_statistics (VlSvm *self)
{
  vl_size i, k ;
  double inner, p ;

  memset(&self->statistics, 0, sizeof(VlSvmStatistics)) ;

  self->statistics.regularizer = self->bias * self->bias ;
  for (i = 0; i < self->dimension; i++) {
    self->statistics.regularizer += self->model[i] * self->model[i] ;
  }
  self->statistics.regularizer *= self->lambda * 0.5 ;

  for (k = 0; k < self->numData ; k++) {
    p = (self->weights) ? self->weights[k] : 1.0 ;
    if (p <= 0) continue ;
    inner = self->innerProductFn(self->data, k, self->model) ;
    inner += self->bias * self->biasMultiplier ;
    self->scores[k] = inner ;
    self->statistics.loss += p * self->lossFn(inner, self->labels[k]) ;
    if (self->solver == VlSvmSolverSdca) {

      self->statistics.dualLoss -= p * self->conjugateLossFn(- self->alpha[k] / p, self->labels[k]) ;
    }
  }

  self->statistics.loss /= self->numData ;
  self->statistics.objective = self->statistics.regularizer + self->statistics.loss ;

  if (self->solver == VlSvmSolverSdca) {
    self->statistics.dualLoss /= self->numData ;
    self->statistics.dualObjective = - self->statistics.regularizer + self->statistics.dualLoss ;
    self->statistics.dualityGap = self->statistics.objective - self->statistics.dualObjective ;
  }
}

/* ---------------------------------------------------------------- */
/*                                       Evaluate rather than solve */
/* ---------------------------------------------------------------- */

void _vl_svm_evaluate (VlSvm *self)
{
  double startTime = vl_get_cpu_time () ;

  _vl_svm_update_statistics (self) ;

  self->statistics.elapsedTime = vl_get_cpu_time() - startTime ;
  self->statistics.iteration = 0 ;
  self->statistics.epoch = 0 ;
  self->statistics.status = VlSvmStatusConverged ;

  if (self->diagnosticFn) {
    self->diagnosticFn(self, self->diagnosticFnData) ;
  }
}

/* ---------------------------------------------------------------- */
/*                         Stochastic Dual Coordinate Ascent Solver */
/* ---------------------------------------------------------------- */

void _vl_svm_sdca_train (VlSvm *self)
{
  double * norm2 ;
  vl_index * permutation ;
  vl_uindex i, t  ;
  double inner, delta, multiplier, p ;

  double startTime = vl_get_cpu_time () ;
  VlRand * rand = vl_get_rand() ;

  norm2 = (double*) vl_calloc(self->numData, sizeof(double));
  permutation = vl_calloc(self->numData, sizeof(vl_index)) ;

  {
    double * buffer = vl_calloc(self->dimension, sizeof(double)) ;
    for (i = 0 ; i < (unsigned)self->numData; i++) {
      double n2 ;
      permutation [i] = i ;
      memset(buffer, 0, self->dimension * sizeof(double)) ;
      self->accumulateFn (self->data, i, buffer, 1) ;
      n2 = self->innerProductFn (self->data, i, buffer) ;
      n2 += self->biasMultiplier * self->biasMultiplier ;
      norm2[i] = n2 / (self->lambda * self->numData) ;
    }
    vl_free(buffer) ;
  }

  for (t = 0 ; 1 ; ++t) {

    if (t % self->numData == 0) {
      /* once a new epoch is reached (all data have been visited),
       change permutation */
      vl_rand_permute_indexes(rand, permutation, self->numData) ;
    }

    /* pick a sample and compute update */
    i = permutation[t % self->numData] ;
    p = (self->weights) ? self->weights[i] : 1.0 ;
    if (p > 0) {
      inner = self->innerProductFn(self->data, i, self->model) ;
      inner += self->bias * self->biasMultiplier ;
      delta = p * self->dcaUpdateFn(self->alpha[i] / p, inner, p * norm2[i], self->labels[i]) ;
    } else {
      delta = 0 ;
    }

    /* apply update */
    if (delta != 0) {
      self->alpha[i] += delta ;
      multiplier = delta / (self->numData * self->lambda) ;
      self->accumulateFn(self->data,i,self->model,multiplier) ;
      self->bias += self->biasMultiplier * multiplier;
    }

    /* call diagnostic occasionally */
    if ((t + 1) % self->diagnosticFrequency == 0 || t + 1 == self->maxNumIterations) {
      _vl_svm_update_statistics (self) ;
      self->statistics.elapsedTime = vl_get_cpu_time() - startTime ;
      self->statistics.iteration = t ;
      self->statistics.epoch = t / self->numData ;

      self->statistics.status = VlSvmStatusTraining ;
      if (self->statistics.dualityGap < self->epsilon) {
        self->statistics.status = VlSvmStatusConverged ;
      }
      else if (t + 1 == self->maxNumIterations) {
        self->statistics.status = VlSvmStatusMaxNumIterationsReached ;
      }

      if (self->diagnosticFn) {
        self->diagnosticFn(self, self->diagnosticFnData) ;
      }

      if (self->statistics.status != VlSvmStatusTraining) {
        break ;
      }
    }
  } /* next iteration */

  vl_free (norm2) ;
  vl_free (permutation) ;
}

/* ---------------------------------------------------------------- */
/*                               Stochastic Gradient Descent Solver */
/* ---------------------------------------------------------------- */

void _vl_svm_sgd_train (VlSvm *self)
{
  vl_index * permutation ;
  double * scores ;
  double * previousScores ;
  vl_uindex i, t, k ;
  double inner, gradient, rate, biasRate, p ;
  double factor = 1.0 ;
  double biasFactor = 1.0 ; /* to allow slower bias learning rate */
  vl_index t0 = VL_MAX(2, vl_ceil_d(1.0 / self->lambda)) ;
  //t0=2 ;

  double startTime = vl_get_cpu_time () ;
  VlRand * rand = vl_get_rand() ;

  permutation = vl_calloc(self->numData, sizeof(vl_index)) ;
  scores = vl_calloc(self->numData * 2, sizeof(double)) ;
  previousScores = scores + self->numData ;

  for (i = 0 ; i < (unsigned)self->numData; i++) {
    permutation [i] = i ;
    previousScores [i] = - VL_INFINITY_D ;
  }

  /*
   We store the w vector as the product fw (factor * model).
   We also use a different factor for the bias: biasFactor * biasMultiplier
   to enable a slower learning rate for the bias.

   Given this representation, it is easy to carry the two key operations:

   * Inner product: <fw,x> = f <w,x>

   * Model update: fp wp = fw - rate * lambda * w - rate * g
                         = f(1 - rate * lambda) w - rate * g

     Thus the update equations are:

                   fp = f(1 - rate * lambda), and
                   wp = w + rate / fp * g ;

   * Realization of the scaling factor. Before the statistics function
     is called, or training finishes, the factor (and biasFactor)
     are explicitly applied to the model and the bias.
  */

  for (t = 0 ; 1 ; ++t) {

    if (t % self->numData == 0) {
      /* once a new epoch is reached (all data have been visited),
       change permutation */
      vl_rand_permute_indexes(rand, permutation, self->numData) ;
    }

    /* pick a sample and compute update */
    i = permutation[t % self->numData] ;
    p = (self->weights) ? self->weights[i] : 1.0 ;
    p = VL_MAX(0.0, p) ; /* we assume non-negative weights, so this is just for robustness */
    inner = factor * self->innerProductFn(self->data, i, self->model) ;
    inner += biasFactor * (self->biasMultiplier * self->bias) ;
    gradient = p * self->lossDerivativeFn(inner, self->labels[i]) ;
    previousScores[i] = scores[i] ;
    scores[i] = inner ;

    /* apply update */
    rate = 1.0 /  (self->lambda * (t + t0)) ;
    biasRate = rate * self->biasLearningRate ;
    factor *= (1.0 - self->lambda * rate) ;
    biasFactor *= (1.0 - self->lambda * biasRate) ;

    /* debug: realize the scaling factor all the times */
    /*
    for (k = 0 ; k < self->dimension ; ++k) self->model[k] *= factor ;
    self->bias *= biasFactor;
    factor = 1.0 ;
    biasFactor = 1.0 ;
    */

    if (gradient != 0) {
      self->accumulateFn(self->data, i, self->model, - gradient * rate / factor) ;
      self->bias += self->biasMultiplier * (- gradient * biasRate / biasFactor) ;
    }

    /* call diagnostic occasionally */
    if ((t + 1) % self->diagnosticFrequency == 0 || t + 1 == self->maxNumIterations) {

      /* realize factor before computing statistics or completing training */
      for (k = 0 ; k < self->dimension ; ++k) self->model[k] *= factor ;
      self->bias *= biasFactor;
      factor = 1.0 ;
      biasFactor = 1.0 ;

      _vl_svm_update_statistics (self) ;

      for (k = 0 ; k < self->numData ; ++k) {
        double delta = scores[k] - previousScores[k] ;
        self->statistics.scoresVariation += delta * delta ;
      }
      self->statistics.scoresVariation = sqrt(self->statistics.scoresVariation) / self->numData ;

      self->statistics.elapsedTime = vl_get_cpu_time() - startTime ;
      self->statistics.iteration = t ;
      self->statistics.epoch = t / self->numData ;

      self->statistics.status = VlSvmStatusTraining ;
      if (self->statistics.scoresVariation < self->epsilon) {
        self->statistics.status = VlSvmStatusConverged ;
      }
      else if (t + 1 == self->maxNumIterations) {
        self->statistics.status = VlSvmStatusMaxNumIterationsReached ;
      }

      if (self->diagnosticFn) {
        self->diagnosticFn(self, self->diagnosticFnData) ;
      }

      if (self->statistics.status != VlSvmStatusTraining) {
        break ;
      }
    }
  } /* next iteration */

  vl_free (scores) ;
  vl_free (permutation) ;
}

/* ---------------------------------------------------------------- */
/*                                                       Dispatcher */
/* ---------------------------------------------------------------- */

void vl_svm_train (VlSvm * self)
{
  assert (self) ;
  switch (self->solver) {
    case VlSvmSolverSdca:
      _vl_svm_sdca_train(self) ;
      break ;
    case VlSvmSolverSgd:
      _vl_svm_sgd_train(self) ;
      break ;
    case VlSvmSolverNone:
      _vl_svm_evaluate(self) ;
      break ;
    default:
      assert(0) ;
  }
}