/** @file svm.c ** @brief Support Vector Machines (SVM) - Implementation ** @author Milan Sulc ** @author Daniele Perrone ** @author Andrea Vedaldi **/ /* 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). */ /** @file svm.h ** @see @ref svm. **/ /** @page svm Support Vector Machines (SVM) @author Milan Sulc @author Daniele Perrone @author Andrea Vedaldi @tableofcontents *Support Vector Machines* (SVMs) are one of the most popular types of discriminate classifiers. VLFeat implements two solvers, SGD and SDCA, capable of learning linear SVMs on a large scale. These linear solvers can be combined with explicit feature maps to learn non-linear models as well. The solver supports a few variants of the standard SVM formulation, including using loss functions other than the hinge loss. @ref svm-starting demonstrates how to use VLFeat to learn an SVM. Information on SVMs and the corresponding optimization algorithms as implemented by VLFeat are given in: - @subpage svm-fundamentals - Linear SVMs and their learning. - @subpage svm-advanced - Loss functions, dual objective, and other details. - @subpage svm-sgd - The SGD algorithm. - @subpage svm-sdca - The SDCA algorithm. @section svm-starting Getting started This section demonstrates how to learn an SVM by using VLFeat. SVM learning is implemented by the ::VlSvm object type. Let's start by a complete example: @code #include #include int main() { vl_size const numData = 4 ; vl_size const dimension = 2 ; double x [dimension * numData] = { 0.0, -0.5, 0.6, -0.3, 0.0, 0.5 0.6, 0.0} ; double y [numData] = {1, 1, -1, 1} ; double lambda = 0.01; double * const model ; double bias ; VlSvm * svm = vl_svm_new(VlSvmSolverSgd, x, dimension, numData, y, lambda) ; vl_svm_train(svm) ; model = vl_svm_get_model(svm) ; bias = vl_svm_get_bias(svm) ; printf("model w = [ %f , %f ] , bias b = %f \n", model[0], model[1], bias); vl_svm_delete(svm) ; return 0; } @endcode This code learns a binary linear SVM using the SGD algorithm on four two-dimensional points using 0.01 as regularization parameter. ::VlSvmSolverSdca can be specified in place of ::VlSvmSolverSdca in orer to use the SDCA algorithm instead. Convergence and other diagnostic information can be obtained after training by using the ::vl_svm_get_statistics function. Algorithms regularly check for convergence (usally after each pass over the data). The ::vl_svm_set_diagnostic_function can be used to specify a callback to be invoked when diagnostic is run. This can be used, for example, to dump information on the screen as the algorithm progresses. Convergence is reached after a maximum number of iterations (::vl_svm_set_max_num_iterations) or after a given criterion falls below a threshold (::vl_svm_set_epsilon). The meaning of these may depend on the specific algorithm (see @ref svm for further details). ::VlSvm is a quite powerful object. Algorithms only need to perform inner product and accumulation operation on the data (see @ref svm-advanced). This is used to abstract from the data type and support almost anything by speciying just two functions (::vl_svm_set_data_functions). A simple interface to this advanced functionality is provided by the ::VlSvmDataset object. This supports natively @c float and @c double data types, as well as applying on the fly the homogeneous kernel map (@ref homkermap). This is exemplified in @ref svmdataset-starting. */ /** @page svm-fundamentals SVM fundamentals @tableofcontents This page introduces the SVM formulation used in VLFeat. See @ref svm for more information on VLFeat SVM support. Let $ \bx \in \real^d $ be a vector representing, for example, an image, an audio track, or a fragment of text. Our goal is to design a *classifier*, i.e. a function that associates to each vector $\bx$ a positive or negative label based on a desired criterion, for example the fact that the image contains or not a cat, that the audio track contains or not English speech, or that the text is or not a scientific paper. The vector $\bx$ is classified by looking at the sign of a *linear scoring function* $\langle \bx, \bw \rangle$. The goal of learning is to estimate the parameter $\bw \in \real^d$ in such a way that the score is positive if the vector $\bx$ belongs to the positive class and negative otherwise. In fact, in the standard SVM formulation the the goal is to have a score of *at least 1* in the first case, and of *at most -1* in the second one, imposing a *margin*. The parameter $\bw$ is estimated or *learned* by fitting the scoring function to a training set of $n$ example pairs $(\bx_i,y_i), i=1,\dots,n$. Here $y_i \in \{-1,1\}$ are the *ground truth labels* of the corresponding example vectors. The fit quality is measured by a *loss function* which, in standard SVMs, is the *hinge loss*: \[ \ell_i(\langle \bw,\bx\rangle) = \max\{0, 1 - y_i \langle \bw,\bx\rangle\}. \] Note that the hinge loss is zero only if the score $\langle \bw,\bx\rangle$ is at least 1 or at most -1, depending on the label $y_i$. Fitting the training data is usually insufficient. In order for the scoring function *generalize to future data* as well, it is usually preferable to trade off the fitting accuracy with the *regularity* of the learned scoring function $\langle \bx, \bw \rangle$. Regularity in the standard formulation is measured by the norm of the parameter vector $\|\bw\|^2$ (see @ref svm-advanced). Averaging the loss on all training samples and adding to it the regularizer weighed by a parameter $\lambda$ yields the *regularized loss objective* @f{equation}{ \boxed{\displaystyle E(\bw) = \frac{\lambda}{2} \left\| \bw \right\|^2 + \frac{1}{n} \sum_{i=1}^n \max\{0, 1 - y_i \langle \bw,\bx\rangle\}. \label{e:svm-primal-hinge} } @f} Note that this objective function is *convex*, so that there exists a single global optimum. The scoring function $\langle \bx, \bw \rangle$ considered so far has been linear and unbiased. @ref svm-bias discusses how a bias term can be added to the SVM and @ref svm-feature-maps shows how non-linear SVMs can be reduced to the linear case by computing suitable feature maps. @ref svm-learning shows how VLFeat can be used to learn an SVM by minimizing $E(\bw)$. @section svm-learning Learning Learning an SVM amounts to finding the minimizer $\bw^*$ of the cost function $E(\bw)$. While there are dozens of methods that can be used to do so, VLFeat implements two large scale methods, designed to work with linear SVMs (see @ref svm-feature-maps to go beyond linear): - @ref svm-sgd - @ref svm-sdca Using these solvers is exemplified in @ref svm-starting. @section svm-bias Adding a bias It is common to add to the SVM scoring function a *bias term* $b$, and to consider the score $\langle \bx,\bw \rangle + b$. In practice the bias term can be crucial to fit the training data optimally, as there is no reason why the inner products $\langle \bx,\bw \rangle$ should be naturally centered at zero. Some SVM learning algorithms can estimate both $\bw$ and $b$ directly. However, other algorithms such as SGD and SDCA cannot. In this case, a simple workaround is to add a constant component $B > 0$ (we call this constant the *bias multiplier*) to the data, i.e. consider the extended data vectors: \[ \bar \bx = \begin{bmatrix} \bx \\ B \end{bmatrix}, \quad \bar \bw = \begin{bmatrix} \bw \\ w_b \end{bmatrix}. \] In this manner the scoring function incorporates implicitly a bias $b = B w_b$: \[ \langle \bar\bx, \bar\bw \rangle = \langle \bx, \bw \rangle + B w_b. \] The disadvantage of this reduction is that the term $w_b^2$ becomes part of the SVM regularizer, which shrinks the bias $b$ towards zero. This effect can be alleviated by making $B$ sufficiently large, because in this case $\|\bw\|^2 \gg w_b^2$ and the shrinking effect is negligible. Unfortunately, making $B$ too large makes the problem numerically unbalanced, so a reasonable trade-off between shrinkage and stability is generally sought. Typically, a good trade-off is obtained by normalizing the data to have unitary Euclidean norm and then choosing $B \in [1, 10]$. Specific implementations of SGD and SDCA may provide explicit support to learn the bias in this manner, but it is important to understand the implications on speed and accuracy of the learning if this is done. @section svm-feature-maps Non-linear SVMs and feature maps So far only linear scoring function $\langle \bx,\bw \rangle$ have been considered. Implicitly, however, this assumes that the objects to be classified (e.g. images) have been encoded as vectors $\bx$ in a way that makes linear classification possible. This encoding step can be made explicit by introducing the *feature map* $\Phi(\bx) \in \real^d$. Including the feature map yields a scoring function *non-linear* in $\bx$: \[ \bx\in\mathcal{X} \quad\longrightarrow\quad \langle \Phi(\bx), \bw \rangle. \] The nature of the input space $\mathcal{X}$ can be arbitrary and might not have a vector space structure at all. The representation or encoding captures a notion of *similarity* between objects: if two vectors $\Phi(\bx_1)$ and $\Phi(\bx_2)$ are similar, then their scores will also be similar. Note that choosing a feature map amounts to incorporating this information in the model *prior* to learning. The relation of feature maps to similarity functions is formalized by the notion of a *kernel*, a positive definite function $K(\bx,\bx')$ measuring the similarity of a pair of objects. A feature map defines a kernel by \[ K(\bx,\bx') = \langle \Phi(\bx),\Phi(\bx') \rangle. \] Viceversa, any kernel function can be represented by a feature map in this manner, establishing an equivalence. So far, all solvers in VLFeat assume that the feature map $\Psi(\bx)$ can be explicitly computed. Although classically kernels were introduced to generalize solvers to non-linear SVMs for which a feature map *cannot* be computed (e.g. for a Gaussian kernel the feature map is infinite dimensional), in practice using explicit feature representations allow to use much faster solvers, so it makes sense to *reverse* this process. */ /** @page svm-advanced Advanced SVM topics @tableofcontents This page discusses advanced SVM topics. For an introduction to SVMs, please refer to @ref svm and @ref svm-fundamentals. @section svm-loss-functions Loss functions The SVM formulation given in @ref svm-fundamentals uses the hinge loss, which is only one of a variety of loss functions that are often used for SVMs. More in general, one can consider the objective @f{equation}{ E(\bw) = \frac{\lambda}{2} \left\| \bw \right\|^2 + \frac{1}{n} \sum_{i=1}^n \ell_i(\langle \bw,\bx\rangle). \label{e:svm-primal} @f} where the loss $\ell_i(z)$ is a convex function of the scalar variable $z$. Losses differ by: (i) their purpose (some are suitable for classification, other for regression), (ii) their smoothness (which usually affects how quickly the SVM objective function can be minimized), and (iii) their statistical interpretation (for example the logistic loss can be used to learn logistic models). Concrete examples are the:
Name Loss $\ell_i(z)$ Description
Hinge $\max\{0, 1-y_i z\}$ The standard SVM loss function.
Square hinge $\max\{0, 1-y_i z\}^2$ The standard SVM loss function, but squared. This version is smoother and may yield numerically easier problems.
Square or l2 $(y_i - z)^2$ This loss yields the ridge regression model (l2 regularised least square).
Linear or l1 $|y_i - z|$ Another loss suitable for regression, usually more robust but harder to optimize than the squared one.
Insensitive l1 $\max\{0, |y_i - z| - \epsilon\}$. This is a variant of the previous loss, proposed in the original Support Vector Regression formulation. Differently from the previous two losses, the insensitivity may yield to a sparse selection of support vectors.
Logistic $\log(1 + e^{-y_i z})$ This corresponds to regularized logisitc regression. The loss can be seen as a negative log-likelihood: $\ell_i(z) = -\log P[y_i | z] = - \log \sigma(y_iz/2)$, where $\sigma(z) = e^z/(1 + e^z)$ is the sigmoid function, mapping a score $z$ to a probability. The $1/2$ factor in the sigmoid is due to the fact that labels are in $\{-1,1\}$ rather than $\{0,1\}$ as more common for the standard sigmoid model.
@section svm-data-abstraction Data abstraction: working with compressed data VLFeat learning algorithms (SGD and SDCA) access the data by means of only two operations: - *inner product*: computing the inner product between the model and a data vector, i.e. $\langle \bw, \bx \rangle$. - *accumulation*: summing a data vector to the model, i.e. $\bw \leftarrow \bw + \beta \bx$. VLFeat learning algorithms are *parameterized* in these two operations. As a consequence, the data can be stored in any format suitable to the user (e.g. dense matrices, sparse matrices, block-sparse matrices, disk caches, and so on) provided that these two operations can be implemented efficiently. Differently from the data, however, the model vector $\bw$ is represented simply as a dense array of doubles. This choice is adequate in almost any case. A particularly useful aspect of this design choice is that the training data can be store in *compressed format* (for example by using product quantization (PQ)). Furthermore, higher-dimensional encodings such as the homogeneous kernel map (@ref homkermap) and the intersection kernel map can be *computed on the fly*. Such techniques are very important when dealing with GBs of data. @section svm-dual-problem Dual problem In optimization, the *dual objective* $D(\balpha)$ of the SVM objective $E(\bw)$ is of great interest. To obtain the dual objective, one starts by approximating each loss term from below by a family of planes: \[ \ell_i(z) = \sup_{u} (u z - \ell_i^*(u) ), \qquad \ell_i^*(u) = \sup_{z} (z u - \ell_i(z) ) \] where $\ell_i^*(u)$ is the *dual conjugate* of the loss and gives the intercept of each approximating plane as a function of the slope. When the loss function is convex, the approximation is in fact exact. Examples include:
Name Loss $\ell_i(z)$ Conjugate loss $\ell_i^*(u)$
Hinge $\max\{0, 1-y_i z\}$ \[ \ell_i^*(u) = \begin{cases} y_i u, & -1 \leq y_i u \leq 0, \\ +\infty, & \text{otherwise} \end{cases} \]
Square hinge $\max\{0, 1-y_i z\}^2$ \[\ell_i^*(u) = \begin{cases} y_i u + \frac{u^2}{4}, & y_i u \leq 0, \\ +\infty, & \text{otherwise} \\ \end{cases}\]
Linear or l1 $|y_i - z|$ \[\ell_i^*(u) = \begin{cases} y_i u, & -1 \leq y_i u \leq 1, \\ +\infty, & \text{otherwise} \\ \end{cases}\]
Square or l2 $(y_i - z)^2$ \[\ell_i^*(u)=y_iu + \frac{u^2}{4}\]
Insensitive l1 $\max\{0, |y_i - z| - \epsilon\}$.
Logistic $\log(1 + e^{-y_i z})$ \[\ell_i^*(u) = \begin{cases} (1+u) \log(1+u) - u \log(-u), & -1 \leq y_i u \leq 0, \\ +\infty, & \text{otherwise} \\ \end{cases}\]
Since each plane $- z \alpha_i - \ell^*_i(-\alpha_i) \leq \ell_i(z)$ bounds the loss from below, by substituting in $E(\bw)$ one can write a lower bound for the SVM objective \[ F(\bw,\balpha) = \frac{\lambda}{2} \|\bw\|^2 - \frac{1}{n}\sum_{i=1}^n (\bw^\top \bx_i\alpha_i + \ell_i^*(-\alpha_i)) \leq E(\bw). \] for each setting of the *dual variables* $\alpha_i$. The dual objective function $D(\balpha)$ is obtained by minimizing the lower bound $F(\bw,\balpha)$ w.r.t. to $\bw$: \[ D(\balpha) = \inf_{\bw} F(\bw,\balpha) \leq E(\bw). \] The minimizer and the dual objective are now easy to find: \[ \boxed{\displaystyle \bw(\balpha) = \frac{1}{\lambda n} \sum_{i=1}^n \bx_i \alpha_i = \frac{1}{\lambda n} X\balpha, \quad D(\balpha) = - \frac{1}{2\lambda n^2} \balpha^\top X^\top X \balpha + \frac{1}{n} \sum_{i=1}^n - \ell_i^*(-\alpha_i) } \] where $X = [\bx_1, \dots, \bx_n]$ is the data matrix. Since the dual is uniformly smaller than the primal, one has the *duality gap* bound: \[ D(\balpha) \leq P(\bw^*) \leq P(\bw(\balpha)) \] This bound can be used to evaluate how far off $\bw(\balpha)$ is from the primal minimizer $\bw^*$. In fact, due to convexity, this bound can be shown to be zero when $\balpha^*$ is the dual maximizer (strong duality): \[ D(\balpha^*) = P(\bw^*) = P(\bw(\balpha^*)), \quad \bw^* = \bw(\balpha^*). \] @section svm-C Parametrization in C Often a slightly different form of the SVM objective is considered, where a parameter $C$ is used to scale the loss instead of the regularizer: \[ E_C(\bw) = \frac{1}{2} \|\bw\|^2 + C \sum_{i=1}^n \ell_i(\langle \bx_i, \bw\rangle) \] This and the objective function $E(\bw)$ in $\lambda$ are equivalent (proportional) if \[ \lambda = \frac{1}{nC}, \qquad C = \frac{1}{n\lambda}. \] up to an overall scaling factor to the problem. **/ /** @page svm-sdca Stochastic Dual Coordinate Ascent @tableofcontents This page describes the *Stochastic Dual Coordinate Ascent* (SDCA) linear SVM solver. Please see @ref svm for an overview of VLFeat SVM support. SDCA maximizes the dual SVM objective (see @ref svm-dual-problem for a derivation of this expression): \[ D(\balpha) = - \frac{1}{2\lambda n^2} \balpha^\top X^\top X \balpha + \frac{1}{n} \sum_{i=1}^n - \ell_i^*(-\alpha_i) \] where $X$ is the data matrix. Recall that the primal parameter corresponding to a given setting of the dual variables is: \[ \bw(\balpha) = \frac{1}{\lambda n} \sum_{i=1}^n \bx_i \alpha_i = \frac{1}{\lambda n} X\balpha \] In its most basic form, the *SDCA algorithm* can be summarized as follows: - Let $\balpha_0 = 0$. - Until the duality gap $P(\bw(\balpha_t)) - D(\balpha_t) < \epsilon$ - Pick a dual variable $q$ uniformly at random in $1, \dots, n$. - Maximize the dual with respect to this variable: $\Delta\alpha_q = \max_{\Delta\alpha_q} D(\balpha_t + \Delta\alpha_q \be_q )$ - Update $\balpha_{t+1} = \balpha_{t} + \be_q \Delta\alpha_q$. In VLFeat, we partially use the nomenclature from @cite{shwartz13a-dual} and @cite{hsieh08a-dual}. @section svm-sdca-dual-max Dual coordinate maximization The updated dual objective can be expanded as: \[ D(\balpha_t + \be_q \Delta\alpha_q) = \text{const.} - \frac{1}{2\lambda n^2} \bx_q^\top \bx_q (\Delta\alpha_q)^2 - \frac{1}{n} \bx_q^\top \frac{X\alpha_t}{\lambda n} \Delta\alpha_q - \frac{1}{n} \ell^*_q(- \alpha_q - \Delta\alpha_q) \] This can also be written as @f{align*} D(\balpha_t + \be_q \Delta\alpha_q) &\propto - \frac{A}{2} (\Delta\alpha_q)^2 - B \Delta\alpha_q - \ell^*_q(- \alpha_q - \Delta\alpha_q), \\ A &= \frac{1}{\lambda n} \bx_q^\top \bx_q = \frac{1}{\lambda n} \| \bx_q \|^2, \\ B &= \bx_q^\top \frac{X\balpha_t}{\lambda n} = \bx_q^\top \bw_t. @f} Maximizing this quantity in the scalar variable $\Delta\balpha$ is usually not difficult. It is convenient to store and incrementally update the model $\bw_t$ after the optimal step $\Delta\balpha$ has been determined: \[ \bw_t = \frac{X \balpha_t}{\lambda n}, \quad \bw_{t+1} = \bw_t + \frac{1}{\lambda n }\bx_q \be_q \Delta\alpha_q. \] For example, consider the hinge loss as given in @ref svm-advanced : \[ \ell_q^*(u) = \begin{cases} y_q u, & -1 \leq y_q u \leq 0, \\ +\infty, & \text{otherwise}. \end{cases} \] The maximizer $\Delta\alpha_q$ of the update objective must be in the range where the conjugate loss is not infinite. Ignoring such bounds, the update can be obtained by setting the derivative of the objective to zero, obtaining \[ \tilde {\Delta \alpha_q}= \frac{y_q - B}{A}. \] Note that $B$ is simply current score associated by the SVM to the sample $\bx_q$. Incorporating the constraint $-1 \leq - y_q (\alpha_q + \Delta \alpha_q) \leq 0$, i.e. $0 \leq y_q (\alpha_q + \Delta \alpha_q) \leq 1$, one obtains the update \[ \Delta\alpha_q = y_q \max\{0, \min\{1, \tilde {\Delta\alpha_q } + y_q \alpha_q\}\} - \alpha_q. \] @section svm-sdca-details Implementation details Rather than visiting points completely at random, VLFeat SDCA follows the best practice of visiting all the points at every epoch (pass through the data), changing the order of the visit randomly by picking every time a new random permutation. **/ /** @page svm-sgd Stochastic Gradient Descent @tableofcontents This page describes the *Stochastic Gradient Descent* (SGD) linear SVM solver. SGD minimizes directly the primal SVM objective (see @ref svm): \[ E(\bw) = \frac{\lambda}{2} \left\| \bw \right\|^2 + \frac{1}{n} \sum_{i=1}^n \ell_i(\langle \bw,\bx\rangle) \] Firts, rewrite the objective as the average \[ E(\bw) = \frac{1}{n} \sum_{i=1}^n E_i(\bw), \quad E_i(\bw) = \frac{\lambda}{2} \left\| \bw \right\|^2 + \ell_i(\langle \bw,\bx\rangle). \] Then SGD performs gradient steps by considering at each iteration one term $E_i(\bw)$ selected at random from this average. In its most basic form, the algorithm is: - Start with $\bw_0 = 0$. - For $t=1,2,\dots T$: - Sample one index $i$ in $1,\dots,n$ uniformly at random. - Compute a subgradient $\bg_t$ of $E_i(\bw)$ at $\bw_t$. - Compute the learning rate $\eta_t$. - Update $\bw_{t+1} = \bw_t - \eta_t \bg_t$. Provided that the learning rate $\eta_t$ is chosen correctly, this simple algorithm is guaranteed to converge to the minimizer $\bw^*$ of $E$. @section svm-sgd-convergence Convergence and speed The goal of the SGD algorithm is to bring the *primal suboptimality* below a threshold $\epsilon_P$: \[ E(\bw_t) - E(\bw^*) \leq \epsilon_P. \] If the learning rate $\eta_t$ is selected appropriately, SGD can be shown to converge properly. For example, @cite{shalev-shwartz07pegasos} show that, since $E(\bw)$ is $\lambda$-strongly convex, then using the learning rate \[ \boxed{\eta_t = \frac{1}{\lambda t}} \] guarantees that the algorithm reaches primal-suboptimality $\epsilon_P$ in \[ \tilde O\left( \frac{1}{\lambda \epsilon_P} \right). \] iterations. This particular SGD variant is sometimes known as PEGASOS @cite{shalev-shwartz07pegasos} and is the version implemented in VLFeat. The *convergence speed* is not sufficient to tell the *learning speed*, i.e. how quickly an algorithm can learn an SVM that performs optimally on the test set. The following two observations can be used to link convergence speed to learning speed: - The regularizer strength is often heuristically selected to be inversely proportional to the number of training samples: $\lambda = \lambda_0 /n$. This reflects the fact that with more training data the prior should count less. - The primal suboptimality $\epsilon_P$ should be about the same as the estimation error of the SVM primal. This estimation error is due to the finite training set size and can be shown to be of the order of $1/\lambda n = 1 / \lambda_0$. Under these two assumptions, PEGASOS can learn a linear SVM in time $\tilde O(n)$, which is *linear in the number of training examples*. This fares much better with $O(n^2)$ or worse of non-linear SVM solvers. @section svm-sgd-bias The bias term Adding a bias $b$ to the SVM scoring function $\langle \bw, \bx \rangle +b$ is done, as explained in @ref svm-bias, by appending a constant feature $B$ (the *bias multiplier*) to the data vectors $\bx$ and a corresponding weight element $w_b$ to the weight vector $\bw$, so that $b = B w_b$ As noted, the bias multiplier should be relatively large in order to avoid shrinking the bias towards zero, but small to make the optimization stable. In particular, setting $B$ to zero learns an unbiased SVM (::vl_svm_set_bias_multiplier). To counter instability caused by a large bias multiplier, the learning rate of the bias is slowed down by multiplying the overall learning rate $\eta_t$ by a bias-specific rate coefficient (::vl_svm_set_bias_learning_rate). As a rule of thumb, if the data vectors $\bx$ are $l^2$ normalized (as they typically should for optimal performance), then a reasonable bias multiplier is in the range 1 to 10 and a reasonable bias learning rate is somewhere in the range of the inverse of that (in this manner the two parts of the extended feature vector $(\bx, B)$ are balanced). @section svm-sgd-starting-iteration Adjusting the learning rate Initially, the learning rate $\eta_t = 1/\lambda t$ is usually too fast: as usually $\lambda \ll 1$, $\eta_1 \gg 1$. But this is clearly excessive (for example, without a loss term, the best learning rate at the first iteration is simply $\eta_1=1$, as this nails the optimum in one step). Thus, the learning rate formula is modified to be $\eta_t = 1 / \lambda (t + t_0)$, where $t_0 \approx 0.5/\lambda$, which is equivalent to start $t_0$ iterations later. In this manner $\eta_1 \approx 1/2$. @subsection svm-sgd-warm-start Warm start Starting from a given model $\bw$ is easy in SGD as the optimization runs in the primal. However, the starting iteration index $t$ should also be advanced for a warm start, as otherwise the initial setting of $\bw$ is rapidly forgot (::vl_svm_set_model, ::vl_svm_set_bias, ::vl_svm_set_iteration_number). @section svm-sgd-details Implementation details @par "Random sampling of points" Rather than visiting points completely at random, VLFeat SDCA follows the best practice of visiting all the points at every epoch (pass through the data), changing the order of the visit randomly by picking every time a new random permutation. @par "Factored representation" At each iteration, the SGD algorithm updates the vector $\bw$ (including the additional bias component $w_b$) as $\bw_{t+1} \leftarrow \bw_t - \eta_t \bw_t - \eta \bg_t$, where $\eta_t$ is the learning rate. If the subgradient of the loss function $\bg_t$ is zero at a given iteration, this amounts to simply shrink $\bw$ towards the origin by multiplying it by the factor $1 - \eta$. Thus such an iteration can be accelerated significantly by representing internally $\bw_t = f_t \bu_t$, where $f_t$ is a scaling factor. Then, the update becomes \[ f_{t+1} \bu_{t+1} = f_{t} \bu_{t} - \eta_t f_{t} \bu_{t} - \eta_t \bg_t = (1-\eta_t) f_{t} \bu_{t} - \eta \bg_t. \] Setting $f_{t+1} = (1-\eta_t) f_{t}$, this gives the update equation for $\bu_t$ \[ \bu_{t+1} = \bu_{t} - \frac{\eta_t}{f_{t+1}} \bg_t. \] but this step can be skipped whenever $\bg_t$ is equal to zero. When the bias component has a different learning rate, this scheme must be adjusted slightly by adding a separated factor for the bias, but it is otherwise identical. **/ /* @section svm-pegasos PEGASOS @subsection svm-pegasos-algorithm Algorithm PEGASOS @cite{shalev-shwartz07pegasos} is a stochastic subgradient optimizer. At the t-th iteration the algorithm: - Samples uniformly at random as subset @f$ A_t @f$ of k of training pairs @f$(x,y)@f$ from the m 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 feature map 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 i-th columns of V). 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 struct VlSvm_ { VlSvmSolverType solver ; /**< SVM solver type. */ vl_size dimension ; /**< Model dimension. */ double * model ; /**< Model ($\bw$ vector). */ double bias ; /**< Bias. */ double biasMultiplier ; /**< Bias feature multiplier. */ /* valid during a run */ double lambda ; /**< Regularizer multiplier. */ void const * data ; vl_size numData ; double const * labels ; /**< Data labels. */ double const * weights ; /**< Data weights. */ VlSvmDataset * ownDataset ; /**< Optional owned dataset. */ VlSvmDiagnosticFunction diagnosticFn ; void * diagnosticFnData ; vl_size diagnosticFrequency ; /**< Frequency of diagnostic. */ VlSvmLossFunction lossFn ; VlSvmLossFunction conjugateLossFn ; VlSvmLossFunction lossDerivativeFn ; VlSvmDcaUpdateFunction dcaUpdateFn ; VlSvmInnerProductFunction innerProductFn ; VlSvmAccumulateFunction accumulateFn ; vl_size iteration ; /**< Current iterations number. */ vl_size maxNumIterations ; /**< Maximum number of iterations. */ double epsilon ; /**< Stopping threshold. */ /* Book keeping */ VlSvmStatistics statistics ; /**< Statistcs. */ double * scores ; /* SGD specific */ double biasLearningRate ; /**< Bias learning rate. */ /* SDCA specific */ double * alpha ; /**< Dual variables. */ } ; /* ---------------------------------------------------------------- */ /** @brief Create a new object with plain data. ** @param type type of SMV solver. ** @param data a pointer to a matrix of data. ** @param dimension dimension of the SVM model. ** @param numData number of training samples. ** @param labels training labels. ** @param lambda regularizer parameter. ** @return the new object. ** ** @a data has one column per sample, in @c double format. ** More advanced inputs can be used with ::vl_svm_new_with_dataset ** and ::vl_svm_new_with_abstract_data. ** ** @sa ::vl_svm_delete **/ 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 ; } /** @brief Create a new object with a dataset. ** @param solver type of SMV solver. ** @param dataset SVM dataset object ** @param labels training samples labels. ** @param lambda regularizer parameter. ** @return the new object. ** @sa ::vl_svm_delete **/ 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 ; } /** @brief Create a new object with abstract data. ** @param solver type of SMV solver. ** @param data pointer to the data. ** @param dimension dimension of the SVM model. ** @param numData num training samples. ** @param labels training samples labels. ** @param lambda regularizer parameter. ** @return the new object. ** ** After calling this function, ::vl_svm_set_data_functions *must* ** be used to setup suitable callbacks for the inner product ** and accumulation operations (@see svm-data-abstraction). ** ** @sa ::vl_svm_delete **/ 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 ; } /** @brief Delete object. ** @param self object. ** @sa ::vl_svm_new **/ 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 */ /* ---------------------------------------------------------------- */ /** @brief Set the convergence threshold ** @param self object ** @param epsilon threshold (non-negative). **/ void vl_svm_set_epsilon (VlSvm *self, double epsilon) { assert(self) ; assert(epsilon >= 0) ; self->epsilon = epsilon ; } /** @brief Get the convergence threshold ** @param self object ** @return epsilon threshold. **/ double vl_svm_get_epsilon (VlSvm const *self) { assert(self) ; return self->epsilon ; } /** @brief Set the bias learning rate ** @param self object ** @param rate bias learning rate (positive). ** ** This parameter applies only to the SGD solver. **/ void vl_svm_set_bias_learning_rate (VlSvm *self, double rate) { assert(self) ; assert(rate > 0) ; self->biasLearningRate = rate ; } /** @brief Get the bias leraning rate. ** @param self object ** @return bias learning rate. **/ double vl_svm_get_bias_learning_rate (VlSvm const *self) { assert(self) ; return self->biasLearningRate ; } /** @brief Set the bias multiplier. ** @param self object ** @param b bias multiplier. ** ** The *bias multiplier* is the value of the constant feature ** appended to the data vectors to implement the bias (@ref svm-bias). **/ void vl_svm_set_bias_multiplier (VlSvm * self, double b) { assert(self) ; assert(b >= 0) ; self->biasMultiplier = b ; } /** @brief Get the bias multiplier. ** @param self object. ** @return bias multiplier. **/ double vl_svm_get_bias_multiplier (VlSvm const * self) { assert(self) ; return self->biasMultiplier ; } /** @brief Set the current iteratio number. ** @param self object. ** @param n iteration number. ** ** If called before training, ** this can be used with SGD for a warm start, as the net ** effect is to slow down the learning rate. **/ void vl_svm_set_iteration_number (VlSvm *self, vl_uindex n) { assert(self) ; self->iteration = n ; } /** @brief Get the current iteration number. ** @param self object. ** @return current iteration number. **/ vl_size vl_svm_get_iteration_number (VlSvm const *self) { assert(self) ; return self->iteration ; } /** @brief Set the maximum number of iterations. ** @param self object. ** @param n maximum number of iterations. **/ void vl_svm_set_max_num_iterations (VlSvm *self, vl_size n) { assert(self) ; self->maxNumIterations = n ; } /** @brief Get the maximum number of iterations. ** @param self object. ** @return maximum number of iterations. **/ vl_size vl_svm_get_max_num_iterations (VlSvm const *self) { assert(self) ; return self->maxNumIterations ; } /** @brief Set the diagnostic frequency. ** @param self object. ** @param f diagnostic frequency (@c >= 1). ** ** A diagnostic round (to test for convergence and to printout ** information) is performed every @a f iterations. **/ void vl_svm_set_diagnostic_frequency (VlSvm *self, vl_size f) { assert(self) ; assert(f > 0) ; self->diagnosticFrequency = f ; } /** @brief Get the diagnostic frequency. ** @param self object. ** @return diagnostic frequency. **/ vl_size vl_svm_get_diagnostic_frequency (VlSvm const *self) { assert(self) ; return self->diagnosticFrequency ; } /** @brief Get the SVM solver type. ** @param self object. ** @return SVM solver type. **/ VlSvmSolverType vl_svm_get_solver (VlSvm const * self) { assert(self) ; return self->solver ; } /** @brief Set the regularizer parameter lambda. ** @param self object. ** @param lambda regularizer parameter. ** ** Note that @a lambda is usually set when calling a ** constructor for ::VlSvm as certain parameters, such ** as the maximum number of iterations, are tuned accordingly. ** This tuning is not performed when @a lambda is changed ** using this function. **/ void vl_svm_set_lambda (VlSvm * self, double lambda) { assert(self) ; assert(lambda >= 0) ; self->lambda = lambda ; } /** @brief Get the regularizer parameter lambda. ** @param self object. ** @return diagnostic frequency. **/ double vl_svm_get_lambda (VlSvm const * self) { assert(self) ; return self->lambda ; } /** @brief Set the data weights. ** @param self object. ** @param weights data weights. ** ** @a weights must be an array of non-negative weights. ** The loss of each data point is multiplied by the corresponding ** weight. ** ** Set @a weights to @c NULL to weight the data uniformly by 1 (default). ** ** Note that the @a weights array is *not* copied and must be valid ** througout the object lifetime (unless it is replaced). **/ void vl_svm_set_weights (VlSvm * self, double const *weights) { assert(self) ; self->weights = weights ; } /** @brief Get the data weights. ** @param self object. ** @return data weights. **/ double const *vl_svm_get_weights (VlSvm const * self) { assert(self) ; return self->weights ; } /* ---------------------------------------------------------------- */ /* Get data */ /* ---------------------------------------------------------------- */ /** @brief Get the model dimenison. ** @param self object. ** @return model dimension. ** ** This is the dimensionality of the weight vector $\bw$. **/ vl_size vl_svm_get_dimension (VlSvm *self) { assert(self) ; return self->dimension ; } /** @brief Get the number of data samples. ** @param self object. ** @return model number of data samples ** ** This is the dimensionality of the weight vector $\bw$. **/ vl_size vl_svm_get_num_data (VlSvm *self) { assert(self) ; return self->numData ; } /** @brief Get the SVM model. ** @param self object. ** @return model. ** ** This is the weight vector $\bw$. **/ double const * vl_svm_get_model (VlSvm const *self) { assert(self) ; return self->model ; } /** @brief Set the SVM model. ** @param self object. ** @param model model. ** ** The function *copies* the content of the vector @a model to the ** internal model buffer. This operation can be used for warm start ** with the SGD algorithm, but has undefined effect with the SDCA algorithm. **/ 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)) ; } /** @brief Set the SVM bias. ** @param self object. ** @param b bias. ** ** The function set the internal representation of the SVM bias to ** be equal to @a b (the bias multiplier ** is applied). The same remark ** that applies to ::vl_svm_set_model applies here too. **/ void vl_svm_set_bias (VlSvm *self, double b) { assert(self); if (self->biasMultiplier) { self->bias = b / self->biasMultiplier ; } } /** @brief Get the value of the bias. ** @param self object. ** @return bias $b$. ** ** The value of the bias returned already include the effect of ** bias mutliplier. **/ double vl_svm_get_bias (VlSvm const *self) { assert(self) ; return self->bias * self->biasMultiplier ; } /** @brief Get the solver statistics. ** @param self object. ** @return statistics. **/ VlSvmStatistics const * vl_svm_get_statistics (VlSvm const *self) { assert(self) ; return &self->statistics ; } /** @brief Get the scores of the data points. ** @param self object. ** @return vector of scores. ** ** After training or during the diagnostic callback, ** this function can be used to retrieve the scores ** of the points, i.e. $\langle \bx_i, \bw \rangle + b$. **/ double const * vl_svm_get_scores (VlSvm const *self) { return self->scores ; } /* ---------------------------------------------------------------- */ /* Callbacks */ /* ---------------------------------------------------------------- */ /** @typedef VlSvmDiagnosticFunction ** @brief SVM diagnostic function pointer. ** @param svm is an instance of ::VlSvm . **/ /** @typedef VlSvmAccumulateFunction ** @brief Pointer to a function that adds to @a model the data point at ** position @a element multiplied by the constant @a multiplier. **/ /** @typedef VlSvmInnerProductFunction ** @brief Pointer to a function that defines the inner product ** between the data point at position @a element and the SVM model **/ /** @brief Set the diagnostic function callback ** @param self object. ** @param f diagnostic function pointer. ** @param data pointer to data used by the diagnostic function. **/ void vl_svm_set_diagnostic_function (VlSvm *self, VlSvmDiagnosticFunction f, void *data) { self->diagnosticFn = f ; self->diagnosticFnData = data ; } /** @brief Set the data functions. ** @param self object. ** @param inner inner product function. ** @param acc accumulate function. ** ** See @ref svm-data-abstraction. **/ void vl_svm_set_data_functions (VlSvm *self, VlSvmInnerProductFunction inner, VlSvmAccumulateFunction acc) { assert(self) ; assert(inner) ; assert(acc) ; self->innerProductFn = inner ; self->accumulateFn = acc ; } /** @brief Set the loss function callback. ** @param self object. ** @param f loss function callback. ** ** Note that setting up a loss requires specifying more than just one ** callback. See @ref svm-loss-functions for details. **/ void vl_svm_set_loss_function (VlSvm *self, VlSvmLossFunction f) { assert(self) ; self->lossFn = f ; } /** @brief Set the loss derivative function callback. ** @copydetails vl_svm_set_loss_function. **/ void vl_svm_set_loss_derivative_function (VlSvm *self, VlSvmLossFunction f) { assert(self) ; self->lossDerivativeFn = f ; } /** @brief Set the conjugate loss function callback. ** @copydetails vl_svm_set_loss_function. **/ void vl_svm_set_conjugate_loss_function (VlSvm *self, VlSvmLossFunction f) { assert(self) ; self->conjugateLossFn = f ; } /** @brief Set the DCA update function callback. ** @copydetails vl_svm_set_loss_function. **/ void vl_svm_set_dca_update_function (VlSvm *self, VlSvmDcaUpdateFunction f) { assert(self) ; self->dcaUpdateFn = f ; } /** @brief Set the loss function to one of the default types. ** @param self object. ** @param loss type of loss function. ** @sa @ref svm-loss-functions. **/ 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 */ /* ---------------------------------------------------------------- */ /** @typedef VlSvmLossFunction ** @brief SVM loss function pointer. ** @param inner inner product between sample and model $\bw^\top \bx$. ** @param label sample label $y$. ** @return value of the loss. ** ** The interface is the same for a loss function, its derivative, ** or the conjugate loss. ** ** @sa @ref svm-fundamentals **/ /** @typedef VlSvmDcaUpdateFunction ** @brief SVM SDCA update function pointer. ** @param alpha current value of the dual variable. ** @param inner inner product $\bw^\top \bx$ of the sample with the SVM model. ** @param norm2 normalization factor $\|\bx\|^2/\lambda n$. ** @param label label $y$ of the sample. ** @return incremental update $\Delta\alpha$ of the dual variable. ** ** @sa @ref svm-sdca **/ /** @brief SVM hinge loss ** @copydetails VlSvmLossFunction */ double vl_svm_hinge_loss (double inner, double label) { return VL_MAX(1 - label * inner, 0.0); } /** @brief SVM hinge loss derivative ** @copydetails VlSvmLossFunction */ double vl_svm_hinge_loss_derivative (double inner, double label) { if (label * inner < 1.0) { return - label ; } else { return 0.0 ; } } /** @brief SVM hinge loss conjugate ** @param u dual variable. ** @param label label value. ** @return conjugate loss. **/ 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 ; } } /** @brief SVM hinge loss DCA update ** @copydetails VlSvmDcaUpdateFunction */ 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 ; } /** @brief SVM square hinge loss ** @copydetails VlSvmLossFunction */ double vl_svm_hinge2_loss (double inner,double label) { double z = VL_MAX(1 - label * inner, 0.0) ; return z*z ; } /** @brief SVM square hinge loss derivative ** @copydetails VlSvmLossFunction */ double vl_svm_hinge2_loss_derivative (double inner, double label) { if (label * inner < 1.0) { return 2 * (inner - label) ; } else { return 0 ; } } /** @brief SVM square hinge loss conjugate ** @copydetails vl_svm_hinge_conjugate_loss */ double vl_svm_hinge2_conjugate_loss (double u, double label) { if (label * u <= 0) { return (label + u/4) * u ; } else { return VL_INFINITY_D ; } } /** @brief SVM square hinge loss DCA update ** @copydetails VlSvmDcaUpdateFunction */ 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 ; } /** @brief SVM l1 loss ** @copydetails VlSvmLossFunction */ double vl_svm_l1_loss (double inner,double label) { return vl_abs_d(label - inner) ; } /** @brief SVM l1 loss derivative ** @copydetails VlSvmLossFunction */ double vl_svm_l1_loss_derivative (double inner, double label) { if (label > inner) { return - 1.0 ; } else { return + 1.0 ; } } /** @brief SVM l1 loss conjugate ** @copydetails vl_svm_hinge_conjugate_loss */ double vl_svm_l1_conjugate_loss (double u, double label) { if (vl_abs_d(u) <= 1) { return label*u ; } else { return VL_INFINITY_D ; } } /** @brief SVM l1 loss DCA update ** @copydetails VlSvmDcaUpdateFunction */ 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 ; } } /** @brief SVM l2 loss ** @copydetails VlSvmLossFunction */ double vl_svm_l2_loss (double inner,double label) { double z = label - inner ; return z*z ; } /** @brief SVM l2 loss derivative ** @copydetails VlSvmLossFunction */ double vl_svm_l2_loss_derivative (double inner, double label) { return - 2 * (label - inner) ; } /** @brief SVM l2 loss conjugate ** @copydetails vl_svm_hinge_conjugate_loss */ double vl_svm_l2_conjugate_loss (double u, double label) { return (label + u/4) * u ; } /** @brief SVM l2 loss DCA update ** @copydetails VlSvmDcaUpdateFunction */ double vl_svm_l2_dca_update (double alpha, double inner, double norm2, double label) { return (label - inner - 0.5*alpha) / (norm2 + 0.5) ; } /** @brief SVM l2 loss ** @copydetails VlSvmLossFunction */ 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) ; } } /** @brief SVM l2 loss derivative ** @copydetails VlSvmLossFunction */ 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) ; } /** @brief SVM l2 loss conjugate ** @copydetails vl_svm_hinge_conjugate_loss */ 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 ; } } /** @brief SVM l2 loss DCA update ** @copydetails VlSvmDcaUpdateFunction */ 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 ; } /* ---------------------------------------------------------------- */ /** @internal @brief Update SVM statistics ** @param self object. **/ 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: = f * 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 */ /* ---------------------------------------------------------------- */ /** @brief Run the SVM solver ** @param self object. ** ** The data on which the SVM operates is passed upon the cration of ** the ::VlSvm object. This function runs a solver to learn a ** corresponding model. See @ref svm-starting. **/ 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) ; } }