How does a Relevance Vector Machine (RVM) work?

Question

Relevance Vector Machines (RVMs) are really interesting models when contrasted with the highly geometrical (and popular) SVMs.

In the light of a question like How does a Support Vector Machine (SVM) work?, and how RVMs are substantially different to SVMs, e.g. What is the difference between Informative (IVM) and Relevance (RVM) vector machines, I think this is a good question to be made.

What are the assumptions made (if any) in RVMs, and what is the general optimization problem?

Answer 1

$\newcommand{\wv}{\mathbf{w}} \newcommand{\alv}{\boldsymbol{\alpha}} \newcommand{\thv}{\boldsymbol{\theta}} \newcommand{\muv}{\boldsymbol{\mu}} \newcommand{\ev}{\mathbf{e}} \newcommand{\fv}{\mathbf{f}} \newcommand{\Xv}{\mathbf{X}} \newcommand{\xv}{\mathbf{x}} \newcommand{\yv}{\mathbf{y}} \newcommand{\vv}{\mathbf{v}} $ The RVM method combines four techniques:

• dual model
• Bayesian approach
• sparsity promoting prior
• kernel trick

The application of this scheme to regression is called Relevance Vector Regression (RVR), and the application to classification is called Relevance Vector Classification (RVC). Since RVC uses logistic regression (or softmax), which is in essence regression, too, the procedure is in principle the same in both cases, which is why I will only describe RVR.

Dual model
Note, that the word "dual" is quite overloaded and is used for many different notions. The usage here follows that in Prince. Consider linear regression: $$ \yv = \Xv\vv + \ev,\tag{1}\label{ols} $$ with $\yv\in\mathbb{R}^n$ the dependent variables in our $n$ observations, $\Xv\in\mathbb{R}^{n\times d}$ the design matrix, i.e. each row contains the $d$ features of a single observation, $\vv\in \mathbb{R}^d$ the parameters that have to be learned, and $\ev\sim N(0,\sigma^2\mathbf 1)$ is Gaussian noise.
Furthermore let's define $S := span(\{\xv\}_{i=1}^n)$, the span of all the observed inputs, i.e. of all the rows of $\Xv$. Then, the dual model uses the fact that $\vv$ can always be chosen as a linear combination of the observations $\mathbf{x}_i$, i.e. $v\in S$. That is because the response $y$ for some input $\xv$ will be computed as the scalar product $y = \langle\vv, \xv\rangle$, and adding a component $\mathbf f$ to $\vv$ with $\mathbf f$ being perpendicular to S, i.e. with $\langle \mathbf f, \xv \rangle = 0$ for all rows $\xv_i$ from $\Xv$, would not change anything, since: $$ \begin{align} \langle\vv+\mathbf f, \xv\rangle &= \langle\vv, \xv\rangle + \langle\mathbf f, \xv\rangle\\ &= \langle\vv, \xv\rangle + 0. \end{align} $$

Therefore, the extra component $\fv$ is not justified by anything in our data $\Xv$ and is thus left out. I.e. it suffices to consider $\vv$ such that there is a (not necessarily unique) $\mathbf{w}$ with: $$ \mathbf{v} = \mathbf{X}^T\mathbf{w}, \tag{2}\label{dual} $$ which turns $\eqref{ols}$ into $$ \yv = \Xv\Xv^T\wv + \ev.\tag{3}\label{dualOls} $$ If there are more dimensions than observations ($d>n$), we have reduced the dimensionality, while, if $d<n$, we have increased. it.
Now, switching from $\mathbf{v}$ to $\mathbf{w}$ means switching to the dual model.

Note, that now each observation $\mathbf{x}_i$ has its own parameter $w_i$. In fact, $w_i$ can be understood as a weight that describes the "relevance" of $\mathbf{x}_i$ for $\mathbf{v}$ in $\eqref{dual}$. Conversely, the smaller the fitted parameter $w_i$ is, the less the influence of $\mathbf{x}_i$ on $\mathbf{v}$. And if $w_i\approx 0$, we see that this observation is of no relevance to the fitted model.

Bayesian approach
Using the Bayesian approach means that we don't just fit the parameters $\wv$ via optimization (that’s the frequentist approach), but rather consider them as random variables of their own. Learning then consists in obtaining their distributions, or at least some approximations thereof. For this approach, we need to choose a prior distribution $p(\wv)$, which contains our prior knowledge of the parameters, independent of the observations. This prior is then combined with the observations to obtain the posterior $p(\wv|\mathbf{X}, \mathbf{y})$. Equipped with the posterior, we can e.g. predict responses $y^\ast$ for new input $\xv^\ast$: $$ p(y^\ast|\Xv, \yv) = \int p(y^\ast|\xv^\ast, \wv)\; p(\wv|\Xv,\yv)\; d\wv.\tag{4}\label{predict} $$

Sparsity promoting prior
This is the most defining part of RVMs. Sparsity is obtained by learning a special prior that assigns more prior probability density to parameters $\mathbf{w}$ that are sparse, i.e. many parameters $w_i$ are zero. This learning procedure is a form of Automatic Relevance Determination (ARD) (MacKay, Neal), and can be understood as follows (see e.g. Tipping, Bishop, or Prince for more details):
The standard choice for a prior of parameters $\mathbf w$ is an isotropic Gaussian $\mathbf w \sim N(0, \alpha^{-1}\mathbf 1)$, where the scalar $\alpha\in\mathbb R$ is the precision. But this would give equal emphasis to all directions, while our goal is to distinguish between coordinates so a subset of them can become zero, giving sparsity. Hence, the major idea for RVMs is to choose as prior $p(\mathbf w)$ a diagonal Gaussian with a different precision $\alpha_i$ for each coordinate $w_i$, i.e. $$ p(\mathbf w) = N(0, \; diag(\alpha_1^{-1}, \ldots, \alpha_n^{-1})). $$ The following procedure is known as maximization of the marginal likelihood (or evidence), which, intuitively, works as follows: Instead of going full Bayesian, i.e. treating all parameters $(\wv, \alv, \sigma^2)$ as random variables, we adopt a "partial" Bayesian approach where only the $\wv$ are random variables while $(\alv, \sigma^2)$ are chosen via optimization. Namely, we pick the hyperparameter pair $(\alv_m, \sigma^2_m)$ that gives the most probability to the observed data $(\Xv, \yv)$ when averaged over all values of $\wv$: $$ \begin{align} (\alv_m, \sigma^2_m) &:= argmax_{(\alv, \sigma^2)} p(\yv| \Xv, \alv, \sigma^2)\\ &= argmax_{(\alv, \sigma^2)} \int p(\yv|\Xv, \wv, \sigma^2)\; p(\wv|\alv) \; d\wv, \end{align} $$ and then use as posterior $p(\wv|\Xv, \yv, \alv_m, \sigma^2_m)$ for predictions (see $\eqref{predict}$). This optimization can be done in various ways and is usually an iterative approximation. Those methods are not completely straightforward and for the gory details, I again refer you to e.g. Tipping, Bishop, or Prince. I will, however, consider the following two aspects:

1. How good an approximation is it?
2. Why does this lead to sparse $\wv$?

With respect to the first question, Tipping points out that setting $(\alv, \sigma^2) = (\alv_m, \sigma_m^2)$ actually does not lead to a good approximation of the posterior. However, we only need it to work well in prediction integrals (see $\eqref{predict}$) and there:

Tipping: "All the evidence from the experiments presented in this paper suggests that this predictive approximation is very effective in general."

For the second question, there exists some nice geometric explanation. Let's consider the space $\mathbb{R}^n$ where the response vectors $\yv\in\mathbb{R}^n$ live. Here, one single point $\yv$ represents all $n$ observed responses $y_i$. The marginal likelihood $p(\yv| \Xv, \alv, \sigma^2)$, which we want to maximize, is a distribution in this space. So our task is to find the distribution that gives the response vector $\yv$ the most probability. It is not difficult to show (see e.g. Tipping, equation (34)), that this marginal likelihood is a family of Gaussians parameterized by $\alv$ and $\sigma^2$ like this: $$ \begin{align} p(\yv| \Xv, \alv, \sigma^2) &= N(\yv | 0, C)\\ C &= \sigma^2\mathbf 1 + \sum_{i=0}^n \alpha_i^{-1}\Xv^c_i {\Xv^c_i}^T, \end{align} $$ where $\Xv^c_i$ is the $i$-th column vector of the design matrix $\Xv$, i.e. the $i$-th components of all the observed $\xv$. The term $\Xv^c_i {\Xv^c_i}^T$ can be understood as the singular variance of a pdf that has positive density only on the line in the direction of $\Xv^c_i$.

Now, first, consider the case where all $\alpha_i^{-1}$ are nearly zero. Then the only relevant part of $C$ is $\sigma^2\mathbf 1$. This means the contour lines of $C$ are circular. That would be the blue circle in the figure below.

Next, imagine that one of the $\alpha_i^{-1}$ is nonzero, say for the index $k$, but the belonging $\Xv^c_k$ is not pointing much in the direction of $\yv$. The sum $C = \sigma^2\mathbf 1 + \alpha_k^{-1}\Xv^c_k {\Xv^c_k}^T$ is then an ellipse that is elongated along $\Xv^c_k$ which makes it leaner and thus it has less density in the direction of $\yv$. This would be the green ellipse in the figure. That means we would want $\alpha_k^{-1}$ to be zero. Thus, for all the columns $\Xv^c_i$ of the design matrix $\Xv$ that are badly aligned with $\yv$, it is not unreasonable to expect that the optimization will drive their precisions $\alpha_i$ to infinity, so the belonging $w_i$ become zero, resulting in sparsity.

Kernel trick
Fortunately, the methods described here only depend on the scalar product of the $\xv$. Thus we can use the well-known kernel trick, meaning we replace the scalar product with a kernel, which is tantamount to using much more complex basis functions.

---

Complexity
One has to distinguish between the training and application complexity of RVM. The application of the trained RVM, i,e. the computation of the prediction $\eqref{predict}$ of the response $\yv^\ast$ of a new input $\xv^\ast$, involves the scalar product of the posterior mean $\muv_\wv$ of $\wv$ with $\xv^\ast$. In the dual model, $\wv$ is $n$-dimensional, so this scalar product is of order $n$, but since it is sparse, the complexity reduces to the order of the number of relevance vectors.

One of the main weaknesses, however, of RVMs, is that the training, in the vanilla version, is cubic in the number of observations. However, there are methods that can reduce this complexity.

Comparison RVM SVM
Superficially, RVMs and SVMs are quite similar. Both use the dual model and both reduce the input set to only a few important observations. However, internally they are quite different. A short comparison might include the following points:

• RVM is Bayesian, SVM is not.
• RVM provides probabilities, SVM doesn’t (in its original variant)
• RVM is sparser and application is faster:

Bishop: “For a wide range of regression and classification tasks, the RVM is found to give models that are typically an order of magnitude more compact than the corresponding support vector machine, resulting in a significant improvement in the speed of processing on test data. Remarkably, this greater sparsity is achieved with little or no reduction in generalization error compared with the corresponding SVM.”

• SVM is a simpler model with a convex objective function, i.e. a global optimum is guaranteed.
• RVM struggles with training complexity:

Bishop: “The principal disadvantage of the relevance vector machine is the relatively long training times compared with the SVM.”

• SVMs achieve sparsity via the maximum margin (classification) or the epsilon-tube (regression) approach, which is geometrically intuitive. RVM, on the other hand, achieves sparsity via special priors and uses a nontrivial approximate optimization of partial posteriors, which is arguably more complex.

The improved sparsity can be observed in the following plot which I reproduced from here using the sklearn-rvm package:

While SVR (left) is using 27 support vectors, RVR (right) uses only ten.