---
abstract: |
    Bayesian max-margin models have shown great superiority in various machine learning tasks with a likelihood regularization, while the probabilistic Monte Carlo sampling for these models still remains challenging, especially for large-scale settings. In analogy to the data augmentation technique to tackle with non-smoothness of the hinge loss, we present a stochastic subgradient MCMC method which is easy to implement and computationally efficient. We investigate the variants that use adaptive stepsizes and thermostats to improve mixing speeds for Bayesian linear SVM. Furthermore, we design a stochastic subgradient HMC within Gibbs method and a doubly stochastic HMC algorithm for mixture of SVMs, a popular extension of linear classifiers. Experimental results on a wide range of problems demonstrate the effectiveness of our approach.

    Many Bayesian models involve continuous but non-differentiable log-posteriors, including the sparse Bayesian methods with a Laplace prior and the regularized Bayesian methods with max-margin posterior regularization that acts like a likelihood term. In analogy to the popular stochastic subgradient methods for deterministic optimization, we present the stochastic subgradient MCMC for efficient posterior inference in such Bayesian models in order to deal with large-scale applications. . Experimental results on a wide range of problems demonstrate the effectiveness of our approach.
author:
- |
    Wenbo Hu, Jun Zhu, Minjie Xu and Bo Zhang\
    State Key Laboratory of Intelligent Technology and Systems (LITS)\
    Tsinghua National Laboratory for Information Science and Technology (TNList)\
    Department of Computer Science and Technology, Tsinghua University, Beijing 100084, China\
    `{hwb13@mails,dcszj@mail}.tsinghua.edu.cn,`\
    `chokkyvista06@gmail.com,dcszb@mail.tsinghua.edu.cn`
title: 'Fast Sampling for Bayesian Max-Margin Models'
---

Introduction
============

In supervised learning, especially classification problems, max-margin learning has proved one of the most important types of methods due to its theoretically well-justified generalization performance [@vapnik2000nature]. Support Vector Machines (SVM) and extensions stand as a good example of this and have been widely used in various discriminative tasks [@altun2003hidden; @cortes1995support]. On the other hand, Bayesian learning methods are also gaining increasing popularity over recent years due to their flexibility to incorporate more delicate latent structures [@gelman2014bayesian] and robustness in the learnt model distributions.

One possible way to combine the above two seemingly irrelevant learning paradigms in order to benefit from a bit of both is called Maximum Entropy Discrimination (MED) [@jaakkola2000maximum]. And it has been more recently extended to the regularized Bayesian inference (RegBayes) framework [@Zhu:jmlr14], with successful applications in topic modeling [@zhu2009medlda], matrix factorization [@xu2012nonparametric], mixture of SVMs for classification [@zhu2011infinite] and Hidden Markov Models [@zhang2014max] for structured prediction.

However, performing posterior inference in such models generally turns out challenging and it is particularly so with MCMC sampling methods because of the hinge loss terms therein that are both non-smooth and hard for direct manipulation. Unsurprisingly, Metropolis-Hastings MCMC with a random walk proposal [@chib1995understanding] is bound to suffer from a fairly low sample efficiency especially in a high dimensional space; While the non-smoothness of the posterior unfortunately disqualifies direct application of the gradient-based Hamiltonian Monte Carlo methods [@neal2011mcmc]. It might be noteworthy that, with a newly discovered data augmentation technique [@polson2011data], people have developed a simple Gibbs sampling algorithm that is arguably the first workable sampling solution for such models. However, built from very local updates, Gibbs sampling is normally slow in convergence, and to be forced to sample from an augmented space only further deteriorates the efficiency.

When it comes to optimizations, a range of subgradient descent methods and their stochastic variants [@shor1985minimization] have been well developed for optimizing non-smooth objectives, examples of which include stochastic subgradient descent (SSGD) solvers for SVMs (Pegasos) [@shai11pegasos], structured SVMs [@Ratliff:2007], sparse Lasso regression [@bertsekas1999nonlinear], as well as several other extensions [@duchi2012projected; @usai2009lasso]. However, none of them has been systematically investigated for performing efficient Bayesian inference.

In this paper, we investigate incorporating subgradients of the posteriors into Hamilton Monte Carlo and come up with subgradient-based MCMC methods for Bayesian models with non-differentiable log-posteriors, i.e. Bayesian max-margin models. We discuss its effect on detailed balance and explore several variants including stochastic subgradient HMC (SSGHMC) which leverages stochastic subsampling of the data [@welling2011bayesian; @chen2014stochastic] to improve efficiency, stochastic subgradient Langevin Dynamics (SSGLD) which uses only one leapfrog step, and stochastic subgradient Nose-Hoover thermostats [@ding2014bayesian] (SSGNHT). We then apply our sampling methods to the mixture of SVMs model and devise a way to effectively integrate out the latent component assignment variables to further improve efficiency over the traditional HMC-within-Gibbs scheme [@meeds2006alternative]. By annealing the stepsizes, our stochastic subgradient MCMC methods can approximately converge to target posteriors fairly efficiently. We note that there have been several previous attempts in using subgradient information in Hamiltonian Monte Carlo or Langevin Monte Carlo [@welling2011bayesian; @neal2011mcmc], yet our work stands as a first close investigation and we also carry out extensive empirical studies on stochastic subgradient MCMC for Bayesian max-margin learning.

The rest of the paper is organized as follows. Section \[sec:background\] reviews Hamiltonian Monte Carlo (HMC) and its extension with stochastic subsampling. Section \[sec:subgradient\] proposes the subgradient HMC method along with its stochastic version. Section \[sec:sampling\_mixture\_of\_svms\] presents a concrete example of the above idea with two specific SSGHMC methods for mixture of SVMs. Section \[sec:experiment\] presents experimental results on several Bayesian max-margin models including Bayesian linear SVMs and mixture of SVMs. Finally, Section \[sec:conclusion\] concludes.

Bayesian methods are becoming increasingly relevant in big data applications to protect rich models from overfitting, to account for uncertainty, and to adaptively infer the model complexity via nonparametric techniques. However, with the fast growth of data volume and model size, developing scalable inference algorithms has become a key step in order to apply Bayesian models. Fortunately, much progress has been done on scalable inference with both variational and Markov chain Monte Carlo (MCMC) methods. We refer the readers to [@zhu2014big] for a review. In particular, stochastic gradient-based MCMC methods have proven effective in exploring high-dimensional sampling spaces, by conjoining the stochastic gradient-based optimization techniques [@robbins1951stochastic] and the Markov chain theory. Representative examples include the stochastic gradient Langevin dynamics (SGLD) [@welling2011bayesian] and its successors on dealing with manifolds [@PatTeh2013a] and exploring higher-order Hamiltonian dynamics [@chen2014stochastic].

However, one implicit assumption typically made in such gradient-based MCMC methods is that the energy function (i.e., log-posteriors) is differentiable, while little work has been done to systematically analyze the continuous but non-differentiable cases, which are not uncommon. For example, sparse Bayesian methods [@park2008bayesian] with a doubly exponential prior is non-smooth in the log-space. Another example arises from Bayesian support vector machines (SVMs) [@polson2011data] and the recent work on max-margin Bayesian latent variable models [@zhu2013g; @Zhu:jmlr14], whose posterior distributions are regularized by the non-smooth hinge loss. For such non-differentiable models, the posterior inference has been largely hindered by the slow batch algorithms. First, the straightforward application of Metropolis-Hastings MCMC with a random walk proposal [@metropolis1953equation] is likely to have very low sample efficiency in high dimensional spaces. Second, the Gibbs samplers with data augmentation techniques [@park2008bayesian; @polson2011data] are not efficient either in high-dimensional spaces as they often involve inverting large matrices. Moreover, the benefit of introducing extra variables would be counteracted in the view of the extra computation on dealing with the extra sampling variables [@roberts2002langevin].

In contrast, stochastic subgradient descent (SSGD) methods [@boyd2008stochastic] have been extensively analyzed in optimizing non-differentiable objectives with many examples, including the stochastic subgradient solvers (Pegasos) for SVMs [@shai11pegasos], the SSGD method for structured SVMs [@Ratliff:2007], and the SSGD methods for sparse Lasso regressor [@bertsekas1999nonlinear] and its various extensions [@duchi2012projected; @usai2009lasso]. However, none of them has been systematically investigated for efficient Bayesian inference.

In this paper, we conjoin the ideas of stochastic subgradient optimization and Markov chain theory and systematically investigate such stochastic subgradient-based MCMC methods for Bayesian models with non-differentiable log-posteriors. By generalizing the Hamiltonian dynamics to the subgradient case, we are able to analyze its properties, including volume preservation and the detailed balance. We further explore the stochastic techniques to reduce the computational cost on calculating the full subgradient by replacing it with an unbiased stochastic estimate. By annealing the stepsizes, our stochastic subgradient MCMC methods can converge efficiently to the target posteriors. We empirically demonstrate the effectiveness on the tasks of learning Bayesian max-margin and sparse Bayesian methods with a diverse range of datasets.

We note that there were several previous attempts to use subgradient information in the Hamiltonian Monte Carlo or Langevin Monte Carlo [@welling2011bayesian; @neal2011mcmc], but little of them gave a rigorous distinction and a theoretical study on subgradient MCMC. One exception is [@NSHMC2015] which distinguished the HMC method for non-smooth posteriors. Compared with [@NSHMC2015], we would like to present the subgradient HMC with a different theoretical study and stochastic subsampling technique. In summary, the stochastic subgradient HMC method is presented with a theoretical guarantee and applied to scalable max margin and sparse applications.

The rest of the paper is organized as follows. Section \[sec:background\] reviews the stochastic Hamiltonian Monte Carlo. Section \[sec:subgradient\] presents the stochastic subgradient Hamiltonian Monte Carlo. Section \[sec:experiment\] presents the experimental results on both max-margin models and sparse Bayesian models. Finally, Section \[sec:conclusion\] concludes.

Preliminaries {#sec:background}
=============

Hamiltonian Dynamics
--------------------

Hamiltonian dynamics is a physical dynamics that generally describes the mechanics [@arnol1989mathematical]. The dynamics is described by two variables, namely the position variable $q$ and the momentum variable $p$. Then the time evolution of the Hamiltonian is determined by the Hamilton’s equation: $$\begin{aligned}
  \frac{{\rm d}q}{{\rm d}t} = \frac{{\partial}H}{{\partial}p}  ~~~~~~~~  \frac{{\rm d}p}{{\rm d}t} = -\frac{{\partial}H}{{\partial}q},\end{aligned}$$ where the Hamiltonian $H$ is a function of $p$ and $q$. If we define $U(q)$ to be the potential energy and $K(p)=p^{\top}M^{-1}p/2$ to be the the kinetic energy where $M$ is a symmetric positive-definite mass matrix, we may decompose the Hamiltonian function into $ H(q,p)=U(q)+K(p)$. The time evolution of the dynamics is $$\label{eqn:Hamiltonian function}
  \frac{{\rm d}q}{{\rm d}t} = M^{-1}p ~~~~~~~~  \frac{{\rm d}p}{{\rm d}t} = -\frac{\partial U}{\partial q}.$$

Hamiltonian Monte Carlo
-----------------------

Hamiltonian Monte Carlo, also known as the Hybrid Monte Carlo (HMC) [@duane1987hybrid; @neal2011mcmc], is one of the classic MCMC methods that combine physical dynamics simulations with statistical sampling.

In the most general sense of the statistical learning, a posterior distribution is any well-defined distribution that takes into account of the prior belief and data [@ghosh2003bayesian]. Besides the conventional Bayes’ rule, a posterior distribution may be derived from other rules, such as regularized Bayesian inference (RegBayes) [@Zhu:jmlr14] which solves an optimization problem with nontrivial constraints.

Formally, we consider the most general sense of statistical learning, where a posterior distribution can be any well-defined distribution that jointly takes into account the prior belief and data. It can thus be represented as $P(\theta | {\mathcal{D}}) \propto \exp( - U(\theta; {\mathcal{D}}) )$, where ${\mathcal{D}}$ is the observation dataset and $U$ is the potential energy function. For conventional Bayesian models, $U$ is typically written as: $$\begin{aligned}
U(\theta;{\mathcal{D}}) = -\log P_0(\theta) - \sum_{i=1}^N \log P(x_i | \theta),\end{aligned}$$ where $P_0(\theta)$ is the prior and $P({\mathcal{D}}| \theta)=\prod_{i=1}^{N}P(x_{i}|\theta)$ is the likelihood given the common assumption. Note that by convention, we use $\theta$ to denote the variables of interest which corresponds to the position variable $q$ in Hamiltonian dynamics. The HMC sampler then simulates the joint distribution over $(\theta, p)$: $$\begin{aligned}
\label{eqn:posterior}
P(\theta, p ) \propto \exp\left( - U(\theta; {\mathcal{D}}) -  p^\top M^{-1} p/2 \right).\end{aligned}$$ Assuming a differentiable potential energy $U(q)$, we can derive an HMC sampler to infer the posterior distribution via simulating the the dynamics  with some discretization methods such as the Euler or leapfrog method. Specifically, using the conventional leapfrog integrator with a stepsize of $\epsilon$, the HMC method updates the iterations through: $$\begin{aligned}
\label{eqn:HMC}
\left\{ \begin{array}{rl}
p_{t+1/2} &= p_{t}-\frac{\epsilon}{2} \nabla U(\theta_{t}|{\mathcal{D}})   \\
\theta_{t+1} &= \theta_{t}-\epsilon M^{-1}p_{t+1/2}\\
p_{t+1} &= p_{t+1/2}-\frac{\epsilon}{2} \nabla U(\theta_{t+1}|{\mathcal{D}}),
\end{array}\right.\end{aligned}$$ where $p_{0}$ is initialized as $p_{0}\sim{\mathcal{N}}(0,M)$. Having obtained samples of $(\theta, p)$, we may simply discard the augmented momentum variable $p$ and get samples from our target posterior $P(\theta | {\mathcal{D}})$.

In particular, if only one leapfrog step is used and $M$ is set to be the identity matrix $I$, then we obtain Langevin Monte Carlo (LMC) by omitting the unused momentum variable [@neal2011mcmc]. $$\begin{aligned}
\label{eqn:LMC}
\left\{ \begin{array}{rl}
\theta_{t+1} &= \theta_{t}-\frac{\epsilon^{2}}{2}\nabla U(\theta_{t+1}|{\mathcal{D}})+\epsilon p_{t} \\
p_{t+1} &= p_{t}-\frac{\epsilon}{2} \nabla U(\theta_{t}|{\mathcal{D}})-\frac{\epsilon}{2} \nabla U(\theta_{t+1}|{\mathcal{D}}).
\end{array}\right.\end{aligned}$$ To compensate for the discretization error, a Metropolis-Hastings correction step is employed to retain the invariance of the target distribution.

Stochastic Gradient HMC
-----------------------

One challenge of these gradient-based MCMC methods on dealing with massive data is that the evaluation of the full gradient $\nabla_{\theta} U(\theta; {\mathcal{D}})$ is often rather expensive to compute. To address this problem, a noisy gradient estimate $\widetilde{\nabla}_{\theta} U(\theta; {\mathcal{D}})$ can be constructed by subsampling the whole dataset, as is the case with stochastic optimization [@robbins1951stochastic; @zhang2004solving].

This idea was first proposed in [@welling2011bayesian] to develop the stochastic gradient Langevin dynamics (SGLD), and was later extended by [@chen2014stochastic] for HMC with friction and by [@ding2014bayesian] for HMC with thermostats. In these stochastic MCMC methods, the gradient of the log-posterior is estimated as $$\begin{aligned}
\label{eqn:noise gradient estimate}
\widetilde{\nabla} U(\theta;{\mathcal{D}}) &=  -\nabla\log P(\theta)-\frac{N}{\widetilde{N}}\nabla\log P(\widehat{{\mathcal{D}}}|\theta),\end{aligned}$$ where $\widehat{{\mathcal{D}}}$ is the randomly drawn subset with size $\widetilde{N}$. Since $\widetilde{N} \ll N$, computing this noisy gradient estimate turns out much cheaper, hence rendering the overall algorithm a lot more scalable.

We now briefly review the stochastic gradient HMC with thermostats, or stochastic gradient Nose-Hoover thermostat (SGNHT). SGNHT uses the simple Euler integrator and introduces augmented thermostat variables to control momentum fluctuations as well as injected noise. The dynamics is simulated via the following equations: $$\begin{aligned}
\label{eqn:SGNHT}
\left\{ \begin{array}{rl}
p_{t+1} &= p_{t}-\epsilon\xi_{t} p_{t}-\epsilon\widetilde{\nabla} U(\theta_{t}|{\mathcal{D}})+\sqrt{2A}{\mathcal{N}}(0,\epsilon)\\
\theta_{t+1} &= \theta_{t}+\epsilon p_{t+1} \\
\xi_{t+1} &= \xi_{t}+\epsilon(\frac{1}{n}p^{\top}_{t}p_{t}-1),
\end{array}\right.\end{aligned}$$ where $\epsilon$ is the stepsize parameter and $A$ is the diffusion factor parameter. $p_{0}$ is initialized as $\mathcal{N}(0,\mathbf{I})$ and $\xi_{0}$ is initialized as $A$.

Similarly, the stochastic gradient LMC [@welling2011bayesian] gives samples through: $$\begin{aligned}
\label{eqn:SGLD}
\left\{ \begin{array}{rl}
\theta_{t+1} &= \theta_{t}-\frac{\epsilon_t^{2}}{2}\widetilde{\nabla} U(\theta_{t}|{\mathcal{D}})+\epsilon_t \eta_t \\
\eta_{t} &\sim  {\mathcal{N}}(0,I),
\end{array}\right.\end{aligned}$$ where $\epsilon_{t}^2$ is the stepsize at step $t$. In [@welling2011bayesian], the polynomially decaying stepsize $\epsilon_{t}^{2} = a(b+t)^{-\gamma}$ is used to omit the MH correction. SGLD can also be considered as a special case of the stochastic gradient Hamiltonian dynamics by setting the leapfrog step as 1.

Stochastic Subgradient MCMC {#sec:subgradient}
===========================

Central to all the above stochastic sampling methods is the (stochastic) gradient of the log-posterior $\nabla_{\theta} U(\theta; {\mathcal{D}})$. However such gradients might not always be available. In this section, we investigate the more general subgradient-based HMC methods

Generalized Hamiltonian Dynamics
--------------------------------

We first show that the non-stochastic subgradient HMC gives the true posterior samples.

Formally, we consider the generalized Hamiltonian dynamics $\mathcal{H}$ with a non-differentiable but continuous potential energy $U(q)$. In this case, we generalize the Hamiltonian equation (\[eqn:Hamiltonian function\]) as $$\label{eqn:subgradient Hamiltonian function}
    \frac{{\rm d}q}{{\rm d}t} = M^{-1}p ~~~~~~~~  \frac{{\rm d}p}{{\rm d}t} = -\mathcal{G}_{q} U(q),$$ where $\mathcal{G}_{q} U(q)$ is the subgradient of $U$ over $q$. If the potential energy is differentiable, the generalized one Hamiltonian dynamics reduces to the ordinary one.

To analyze the properties of the generalized Hamiltonian dynamics, we restrict ourself to the case where the energy function (i.e., log-posterior) is continuous and *almost differentiable* [@shor1985minimization] which means that the set of non-differentiable points does not have any accumulation points. This is of wide interest and covers the max-margin Bayesian models and the sparse Bayesian methods with a Laplace prior.

In our scheme, we adpot the separable Hamiltonian here for simplicity and the separable Hamiltonian also assumes the following important properties. We show that the generalized Hamiltonian dynamics shares the similar properties [@neal2011mcmc] as the ordinary Hamiltonian dynamics with a differentiable energy function if the Hamiltonian is separable. These properties are significant in constructing the stochastic subgradient sampling methods from the generalized Hamiltonian dynamics.

### Reversibility

Reversibility means that any mapping from one state to another state corresponds to an inverse mapping which turns back to the primary state. For the ordinary Hamiltonian dynamics, the reversibility is proved by simply negating the Hamiltonian equation \[eqn:Hamiltonian function\]. For the generalized Hamiltonian dynamics, we still can show the reversibility by negating the Hamiltonian equation \[eqn:subgradient Hamiltonian function\]. This property indicates the reversibility of the Markov chain transitions of the subgradient HMC.

### Conservation of the Hamiltonian

Even though the potential energy field is non-differentiable, we can still show that the Hamiltonian is kept invariant: $$\begin{aligned}
\label{eqn:InvariantHamiltonian}
  \mathcal{G}_{t}H  =\sum_{i}\left\{ \frac{\partial K}{\partial p_i}\mathcal{G}_{q_{i}} U-\mathcal{G}_{q_{i}}U \frac{\partial K}{\partial p_i} \right\} = 0.\end{aligned}$$

### Volume Preservation

The third significant property of the generalized Hamiltonian dynamics is the volume preservation(also known as Liouville’s theorem for ordinary Hamiltonian dynamics) which means that any phase flow does not change volume.

\[pro:volume preservation\] Volume preservation holds for $\mathcal{H}$, if the potential energy $U$ is continuous and almost differentiable.

See Appendix \[Proof of the volume preservation\].

The basic idea of this proof is to use the phase flow to describe the transformation of the phase-space volume and the separable Hamiltonian assumes the preservation regardless of the non-smoothness [@arnol1989mathematical]. This property is significant for proving detailed balance for the Metropolis update of the subgradient HMC.

Subgradient Hamiltonian Monte Carlo and Its Approximated Detailed Balance
-------------------------------------------------------------------------

When the log-posterior is non-differentiable, gradient-based HMC is not directly applicable. Using the more general subgradients could potentially address this problem, as proven in previous deterministic subgradient descent methods [@shor1985minimization].

Plugging the subgradient of the potential energy $U(q)$, $\mathcal{G}_{q} U(q)$, in the ordinary HMC, we come up with the subgradient HMC with leapfrog method as: $$\begin{aligned}
\label{eqn:subgradient HMC}
\left\{ \begin{array}{rl}
p_{t+1/2} &= p_{t}-\frac{\epsilon}{2} \mathcal{G} U(\theta_{t}|{\mathcal{D}})   \\
\theta_{t+1} &= \theta_{t}-\epsilon M^{-1}p_{t+1/2}\\
p_{t+1} &= p_{t+1/2}-\frac{\epsilon}{2} \mathcal{G} U(\theta_{t+1}|{\mathcal{D}}),
\end{array}\right.\end{aligned}$$ where $p_{0}$ is initialized as $p_{0}\sim {\mathcal{N}}(0,M)$ and $\epsilon$ is the discretization stepsize.

From a theoretical perspective, we may not be able to readily analyze the volume preservation property of the Hamiltonian dynamics with a non-differentiable potential energy nor the detailed balance of a general subgradient HMC sampler. Instead, we give an approximated analysis based on several practical assumptions of the potential energy.

In practical Bayesian models, the non-smoothness of the posterior often lies in the Laplace prior or the hinge loss induced likelihoods which are mainly considered in this paper. Both of these cases are *almost differentiable* which means they are only non-differentiable on a zero mass set, a direct implication of which being that when doing HMC sampling with such dynamics, the sampler will only hit those non-differentiable states with probability $0$. For such dynamics $H_{0}$, we may construct a series of ordinary Hamiltonian dynamics $H_{\epsilon}$ via polynomial interpolation. The resulting dynamics series $H_{\epsilon}$ has the same potential energy function except for a $\epsilon$ mass set. We defer the construction details to Appendix \[appsec:polynomial construction\] to save space.

When drawing at most a countable number of samples from $H_{0}$, the subgradient HMC can be thought of as drawing samples from an ordinary $H_{\epsilon}$ instead. With this approximation, the sampler satisfies detailed balance and is thus valid for generating approximate samples from $H_{0}$. We can make such approximation as accurate as possible by simply letting $\epsilon$ go to zero in the sense of the converging posterior.

With the properties of the generalized Hamiltonian dynamics, we can show the detailed balance of the subgradient HMC as follows:

The subgradient HMC with an MH test and a discretization integrator (either Euler or leapfrog) leaves the canonical distribution invariant.

See Appendix \[Proof of the detailed balance\].

The detailed balance of the subgradient HMC means that the Metropolis update of the subgradient HMC proposal leaves the canonical distribution invariant.

Stochastic Subgradient MCMC in Practice
---------------------------------------

We can readily obtain the straightforward version of stochastic subgradient Langevin dynamics (SSGLD) by simply replacing the gradient of the log-posterior with its subgradient. More formally, SSGLD generates samples by simulating the following dynamics: $$\begin{aligned}
\label{eqn:SSGLD}
\left\{ \begin{array}{rl}
\theta_{t+1} &= \theta_{t}-\frac{\epsilon^{2}}{2}\widetilde{\mathcal{G}} U(\theta_{t+1}|{\mathcal{D}})+\epsilon \eta_t \\
\eta_{t} &\sim  {\mathcal{N}}(0,I),
\end{array}\right.\end{aligned}$$ where $\widetilde{\mathcal{G}}U(\theta;{\mathcal{D}})$ is the stochastic noisy estimate of the true subgradient information $\mathcal{G} U(\theta;{\mathcal{D}})$, $$\label{eqn:subgradient estimate}
\widetilde{\mathcal{G}} U(\theta;{\mathcal{D}}) =  -\mathcal{G}\log P(\theta)-\frac{N}{\widetilde{N}}\mathcal{G}\log P(\widehat{{\mathcal{D}}}|\theta).$$ In existing SGLD methods [@welling2011bayesian], it is recommended to use a polynomial decaying stepsize to save the MH correction step of the Langevin proposals. When the stepsize properly decays, the Markov chain would gradually converge to the target posterior. One subtle part of the method is thus on tuning the discretization stepsize. A pre-specified annealing scheme (if not chosen properly) would make the chain either miss or oscillate around the target. More recent work [@TehThiVol2014a] recommends some relatively optimal scheme for SGLD. Inspired by adaptive stepsizes for (sub)gradient descent methods [@duchi2011adaptive], we in this paper adopt the same adaptive stepsize setting for our SSGLD methods. As we shall see in the experiments, such scheme is indeed beneficial to yielding faster mixing speeds.

data $X$, batchsize $\widetilde{N}$, stepsize $\epsilon$ Initialize $\theta_{0}$. $\theta_{t+1} = \theta_{t}-\frac{\epsilon^{2}}{2}\widetilde{\mathcal{G}} U(\theta_{t+1}|X)+\epsilon {\mathcal{N}}(0,I)$

We can derive stochastic subgradient Hamiltonian dynamics likewise. Here, we adopt an improved version of stochastic gradient HMC [@ding2014bayesian] to derive our stochastic subgradient Nose-Hoover thermostat (SSGNHT), which generates samples via $m$ steps of the following iterations: $$\begin{aligned}
\label{eqn:SSGNHT}
\left\{ \begin{array}{rl}
p_{t+1} &= p_{t} - \xi_t p_{t}\epsilon - \epsilon\widetilde{\mathcal{G}} U(\theta_{t}|X)+\sqrt{2A} {\mathcal{N}}(0,\epsilon)\\
\theta_{t+1} &= \theta_{t}+\epsilon p_{t+1} \\
\xi_{t+1} &= \xi_{t}+\epsilon(\frac{1}{n}p^{\top}_{t}p_{t}-1).
\end{array}\right.\end{aligned}$$ Again we omit the MH correction step by adopting adaptive annealing stepsizes. With the properly chosen stepsizes and thermostats initialization, the SSGNHT simulation would generate posterior samples more efficiently.

Stochastic Subgradient MCMC for Binary Linear SVMs {#sec:bls}
--------------------------------------------------

The above stochastic subgradient MCMC can be directly applied to performing approximate posterior inference in Bayesian linear SVMs. Let ${\mathcal{D}}= \{(x_i, y_i)\}_{i=1}^N$ be the given training dataset, where $x_i \in \mathbb{R}^d$ denotes the $d$-dimensional feature vector of the $i$th instance and $y_i \in \{+1, -1\}$ the binary label. We consider linear classifiers with a weight vector $\eta\in\mathbb{R}^d$ and the decision rule is naturally $\hat{y} = \textrm{sign}(\eta^\top x)$. Then in Bayesian linear SVMs, we are interested in learning a posterior distribution $p(\eta | {\mathcal{D}})$ of the classifier weight. The prior is chosen to be a standard normal distribution $p_{0}(\eta) = {\mathcal{N}}(0, I)$ while the per-datum unnormalized likelihood is set to $p(y_i | x_i, \eta) = \exp(-c \cdot \max(0,1-y_{i}\eta^{T}x_{i}))$ ($c$ is the nonnegative regularization parameter), each of which respectively corresponds to the 2-norm regularizer and hinge loss in linear SVMs.

Taking the subgradient of the log-posterior, or more specifically, the non-differentiable log-likelihood, we obtain $$\begin{aligned}
\mathcal{G}\log P(y_i | x_i, \eta) =
\begin{cases}
- c y_{i}x_{i} & 1-y_{i}\eta^{\top}x_{i} > 0  \\
0 & 1-y_{i}\eta^{\top}x_{i} \leq 0,
\end{cases}\end{aligned}$$ with which we can directly adopt the above stochastic subgradient MCMC samplers Eqn. \[eqn:SSGLD\] and Eqn. \[eqn:SSGNHT\] for approximate inference.

data $X$, batchsize $\widetilde{N}$, stepsize $\epsilon$, step number $m$, $\theta_{0}$. Initialize $p_{0}={\mathcal{N}}(0,I)$, $\xi_{0}=A$ $p_{t+1} = p_{t}-\epsilon\xi_{t} p_{t}-\epsilon\widetilde{\mathcal{G}} U(\theta_{t}|X)+\sqrt{2A} {\mathcal{N}}(0,\epsilon)$ $\theta_{t+1} = \theta_{t}+\epsilon p_{t+1}$ $\xi_{t+1} = \xi_{t}+\epsilon(\frac{1}{n}p^{\top}_{t}p_{t}-1)$

Fast Sampling for Mixture of SVMs {#sec:sampling_mixture_of_svms}
=================================

We now show how to leverage the above stochastic subgradient MCMC methods to derive a fast sampling algorithm for Bayesian mixture of SVMs, a nonlinear extension to Bayesian linear SVMs.

Mixture of SVMs
---------------

Real world data often have some latent clustering structures in them. Therefore mixture of SVMs becomes a natural extension, since mixture-of-experts models are generally capable of capturing these local structures and consequently learn a non-linear model instead of simply a linear one. Here we give an outlined description of Bayesian parametric mixture of SVMs and its nonparametric extension, infinite SVM [@zhu2011infinite].

In the mixture model, we associate each datum with a latent component assignment variable $z$ and then, based on the value of $z$, the input feature $x$ is generated from a corresponding Gaussian distribution $\mathcal{N}(\mu_{z},\Sigma_{z})$. The weight $\eta$ of the corresponding local SVM is generated similarly as the above Bayesian linear SVMs.

Here, the Gibbs classifiers works as follows: given the posterior distribution $q(Z,\eta)$, the Gibbs classifier draws a component assignment $z_{i}$ and a classifier model $\eta_{z_{i}}$ for each data point $i$ and takes the following discriminative function, $$\hat{y}(x_{i})=\mathop{\rm sgn}(\eta_{z_{i}}^{\top}x_{i}).$$ In the Gibbs classifier, the expected hinge loss $\mathcal{R}'=\sum_{i}\mathbb{E}_{q}[\max(0,l-y_{i}\eta_{z_{i}}x_{i})]$ is adopted as an alternative loss function. Then we can build a regularized Bayesian model for mixture model with a Gibbs classifier $$\min_{q(Z,\eta,\gamma)\in\mathcal{P}}\mathcal{L}(q(Z,\eta,\gamma)+2c\cdot\mathcal{R}'(q(Z,\eta,\gamma)))$$ where $\gamma=(\mu,\Sigma)$ is the mean and variance parameter for each Gaussian component, $\mathcal{L}=\mathrm{KL}(q||p_{0}(X,\eta,\gamma))-\mathbb{E}(\log(p(X|Z,\gamma))$ is the objective function when doing standard Bayesian inference. Existing Gibbs sampling methods adopt a Gaussian-Inverse-Wishart conjugate prior for $\gamma$ and use the data augmentation technique to perform a Gibbs sampler for the following posterior distribution $$q(Z,\eta,\gamma)\propto p_{0}(\eta,\gamma,Z)p(X|Z,\gamma)\phi(y|Z,\eta).$$ In infinite SVM, we usually choose a Gaussian prior for $\eta$ and a *Chinese Restaurant Process* (CRP) prior for $Z$. We also fix the number of components to simplify the probabilistic inference and call it *parametric mixture of SVMs*. We defer readers to some previous literature [@zhu2011infinite; @zhang2014max] to see the details.

HMC within Gibbs and Doubly Stochastic Sampling
-----------------------------------------------

[R]{}[0.6]{}

data $(x_{i},y_{i}),i=1,\cdots,N$, batchsize $\widetilde{N}$. Initialization sample $z$ given $\eta$ sample $\eta$ given $z$ using stochastic subgradient HMC

HMC within Gibbs is a sampling scheme that takes the structure of the partial probability and introduces the HMC to some of the Gibbs sampling steps. After collapsing the Gaussian component parameters, the stochastic subgradient HMC within Gibbs can be developed for the parametric mixture of SVMs. The sampling step of associated classifier coefficients are sampled via a stochastic subgradient MCMC. Using the CRP formulation, the method can similarly be used for the infinite SVM model.

[R]{}[0.6]{}

data $(x_{i},y_{i})$, batchsize $\widetilde{N}$¡£ Initialization Randomly draw $\widetilde{N}$ data points calculate $p(z_{i}=k|x_{i},y_{i},\gamma,\eta)$ for $\widetilde{N}$ data points sample $\eta$,$\gamma$ for each component with $\widetilde{N}$ data points

In the Gibbs sampling structure, the step of sampling component assignment $z$ is difficult to take the stochastic subsampling technique. Moreover, the subchains in the Gibbs classifier also slow down the mixing of the Gibbs chain. Inspired by the doubly stochastic method for variational method [@titsias2014doubly], we can develop a doubly stochastic sampler for the mixture model by taking the expectation of assignment $Z$ when taking stochastic gradients over classifier weights $\eta$. In doubly stochastic sampler, we sample both discrete assignment $z$ (in an expectation form) and continuous parameters with stochastic subsampling. Formally, we take the stochastic subgradient of $q(\eta,\gamma)=\sum_{Z}q(Z,\eta,\gamma)$ over $\eta$ as $$\frac{\partial}{\partial \eta_{k}} \log q(\eta,\gamma) = \frac{\partial}{\partial \eta_{k}} \log p_{0}(\eta_{k}) +  \sum_{i}p(z_{i}=k|x_{i},y_{i},\gamma,\eta)\mathcal{G}_{\eta_{k}}\log\phi_{i}(y_{i},x_{i}|z,\eta_{k}).$$ We put the partial gradients over $\gamma$ in the Appendix. \[appsec:doubly stochastic\] to save space.

When evaluating the probability of the component assignment $p(z_{i}=k|x_{i},y_{i},\gamma,\eta)$, we only need to deal with one tiny batchsize which would greatly improves the sampling efficiency (see Alg. \[alg:doubly stochastic\]). Note that, currently we only apply this doubly stochastic for the parametric mixture of SVMs and leave a future extension for infinite SVM.

Experiments {#sec:experiment}
===========

We now empirically show how our stochastic subgradient sampling methods work on two Bayesian max-margin models, namely Bayesian linear SVMs [@vapnik2000nature] and Bayesian mixture of SVMs [@zhu2011infinite]. Our results demonstrate that stochastic subgradient MCMC can achieve great improvement on time efficiency while still getting accurate posterior samples.

Bayesian Linear SVMs
--------------------

We first implement our stochastic subgradient MCMC for Bayesian linear SVMs. The work [@polson2011data] presents a Gibbs sampler by using data augmentation techniques, but it involves inversions of the covariance matrices that are of dimension $n \times n$, which is not scalable for high-dimensional data. We treat this method as a strong baseline. We also compare with random walk Metropolis [@chib1995understanding] on one synthetic dataset and two real datasets. More specifically, we consider stochastic MH test [@korattikara2014austerity] and hence call it *stochastic random walk Metropolis* (SRWM).

### Results on Synthetic Data

We first test on a two dimensional synthetic dataset to show that our methods give the correct samples from the posterior distribution. Specifically, we generate 1,000 input observations from a uniform distribution $x_i \stackrel{i.i.d} \sim U(0,1)$ and the two-dimensional coefficient vector as $\eta \sim {\mathcal{N}}(0,\lambda^{-1} I)$, where $\lambda = 3$. Given the input observations and coefficients, the labels are generated from the Bernoulli distribution as detailed above in Sec. \[sec:bls\].

We compare the samples obtained from SSGLD and SSGNHT with those from the data augmentation method, which is known as an accurate sampler for Bayesian SVMs. We set the batchsize $\widetilde{N}=10$ for stochastic subgradient MCMC methods. As the data is rather simple, we just set the stepsize as small as 0.0001 for SSGLD and 0.0005 for SSGNHT and omit MH corrections. Besides, for SSGNHT we set the diffusion parameter $A=1$ and the step number $m=20$. We take 5,000 samples for each method after a sufficiently long burn-in stage and give the comparison in Figure \[fig:posterior contrast\], where the posterior sample mean and the principal directions of the samples are also marked. Such visual evidence illustrates that our stochastic subgradient MCMC methods are very accurate and that the difference between samples obtained from our methods and those from the data augmentation sampler is almost indistinguishable.

### Results on Real Data

We then test our methods on two real world datasets, namely, the low dimensional UCI Higgs dataset and the high dimensional Realsim dataset [^1].

The Higgs dataset contains $1.1\times 10^7$ samples in a 28-dimensional feature space. We randomly choose $10^7$ samples as the training set and use the rest as the testing set. We select the stochastic batchsize $\widetilde{N}$ to be $1,000$ for both SSGLD and stochastic random walk Metropolis. For SSGHNT, we set $\widetilde{N}=100$, which is a good setting to be justified below in Section \[sec:sensitivity\]. We use the adaptive stepsize for SSGLD, which has been successfully applied in the stochastic (sub)gradient descent [@duchi2011adaptive], and use the polynomial decaying stepsize $0.0001\times t^{-0.2}$ for SSGNHT. We choose the diffusion factor $A=1$ and step number $m=20$ for SSGNHT. For the stochastic random walk Metropolis, the variance parameter is $0.01$. The tuning of these parameters is carried out with cross validation and we discuss in more details how the batchsize parameter is chosen in Section \[sec:sensitivity\].

Figure \[fig:svmhiggs\] shows the convergence curves of various methods on the Higgs dataset with respect to the running time. We can see that our stochastic subgradient MCMC methods spend less time before reaching the best accuracy obtained by linear SVMs. Moreover, although stochastic subgradient MCMC methods often take more iterations in order to get a high accuracy than the data augmentation method, our methods are overall more efficient because at each iteration they draw posterior samples using only a very small mini-batch of the large dataset. Furthermore, compared with stochastic random walk Metropolis, the subgradient information used in our subgradient MCMC methods provides the right direction to the true posterior. Finally, we observe that the two stochastic subgradient MCMC methods have comparable convergence speeds, although in theory SSGNHT would mix faster than SSGLD by using the momenta information [@neal2011mcmc; @chen2014stochastic].

We then test different methods on the high-dimensional Realsim dataset, which contains $92,309$ observations in a $20,958$ dimension feature space. We randomly draw $50,000$ observations as the training data and use the rest as the testing data. We also take stochastic random walk Metropolis and the data augmentation sampler as baseline methods. We take $10$ as the stochastic batchsize for the three stochastic sampling methods. For stochastic random walk Metropolis, we set the variance parameter to $0.05$. For SSGNHT, we choose diffusion factor $A=1$, decaying stepsize $0.001\times t^{-0.2}$ and the step number $m=10$. We use adaptive stepsizes for SSGLD. Figure \[fig:svmrealsim\] shows the convergence of the above methods with respect to the running time in the log scale. We can see that our methods are about ten times faster than the two baseline methods, which might have missed those more efficient sampling directions particularly due to the curse of dimensionality. In contrast, our stochastic subgradient MCMC methods give the more practical solution to the challenging inference task in a high dimensional space.

In our experiments, SSGLD with adaptive stepsizes seem to have a better mixing than the polynomial decaying stepsizes. Better mixing may be resulted from the flexible stepsize decaying at different dimension. We also give an empirical analysis on adaptive stepsizes for SSGLD in Appendix \[adagrad\].

### Sensitivity Analysis {#sec:sensitivity}

Figure \[fig:svmbatchsize\] presents the sensitivity analysis of the batchsize $\widetilde{N}$ for the two stochastic subgradient MCMC methods on both the Higgs and Realsim datasets. We can see that tuning the batchsize $\widetilde{N}$ reflects an accuracy-efficiency trade-off, analogous to the bias-variance tradeoff in stochastic Monte Carlo sampling [@korattikara2014austerity]. In general, using a smaller batchsize often leads to a larger injected noise, but the computation cost at each iteration is reduced, which is linear to the batchsize (i.e., $O(\widetilde{N}$). When doing cross validation to select parameters, both the accuracy and time efficiency are key factors that should be taken into considerations. Take the results of SSGNHT on the Higgs data as an example (i.e., the upper right subgraph in Figure \[fig:svmbatchsize\]). We can see that using the batchsize of 1,000 (i.e., the circled curve) leads to a slightly more stable curve than using a batchsize of 100 (i.e., the dotted curve). When we take both the accuracy and efficiency into consideration, the mild fluctuation of the green line is not that serious, especially when we consider the computation benefits resulting from the smaller batchsize. Therefore, we choose the batchsize as 100 for the results shown in Figure \[fig:svmhiggs\].

A typical example, which shows the balance of the accuracy-efficiency tradeoff, lies in the batchsize parameter choosing of the SSGNHT for binary SVM model on Higgs dataset. In the upper right subgraph of the Figure \[fig:svmbatchsize\], the red line which represents the batchsize 1000, is a bit more stable than the green line which represents batchsize 100. When we take the both accuracy and efficiency into consideration, the mild fluctuation of the green line are counteracted by the computation benefit resulted from the smaller stochastic batchsize. Therefore, we choose the batchsize as 100 in the results shown in Figure\[fig:svmhiggs\]. .

Mixture of SVMs
---------------

For the latent mixture model, we first consider the parametric mixture of SVMs with a fixed number of components. The doubly stochastic subgradient HMC for parametric mixture of SVM and Gibbs samping [@zhang2014max] are implemented for a comparison.

  Model               Accu(%)      Training Time    Time for sampling $\eta$    Time for sample z
  --------------- --------------- ---------------- -------------------------- -------------------
  F-GiSVM+Gibbs    $91.7\pm0.0$    $ 39.86\pm1.6$        $26.02\pm1.21$            $13.85\pm0.52$
  F-GiSVM+HMC      $92.5\pm 0.1$   $0.77\pm0.02$               -                                -

  : Efficiency (in minutes) and accuracy of various methods for parametric mixture of SVMs on the IJCNN dataset.[]{data-label="table:mixture of svm for ijcnn"}

We choose the IJCNN dataset [^2] which contains $49,990$ training samples and $91,701$ test samples in a 22-dimensional space. We use $5,000$ as the batchsize and use adaptive stepsize for each method.

Table \[table:mixture of svm for ijcnn\] shows the final classification performance and the time efficiency of each method. From the results, we can see that with the doubly stochastic HMC sampler, the time efficiency is greatly improved.

Infinite SVMs
-------------

Finally, we implement the Gibbs sampler and the stochastic HMC within Gibbs for the infinite SVM model using the IJCNN dataset. We also use the batchsize $5,000$ and adaptive stepsize for HMC sampling. From the results in Table \[table:isvm for ijcnn\], we can see that the time for sampling $\eta$ is significantly reduced.

We should also note that the stochastic HMC within Gibbs would still be trapped when sampling $z$. Currently, the doubly stochastic sampler for infinite SVMs still remains untested and we leave it as a future extension.

  Model             Accu(%)      Training Time    Time for sampling $\eta$    Time for sample z
  ------------- --------------- ---------------- -------------------------- -------------------
  GiSVM+Gibbs    $94.2\pm0.9$    $68.22\pm1.13$        $24.56\pm0.93$            $43.57\pm1.68$
  GiSVM+HMC      $94.1\pm 0.9$   $48.95\pm3.87$        $1.10\pm0.06$             $47.75\pm3.83$

  : Efficiency (in minutes) and accuracy of various methods for Infinite SVMs on the IJCNN dataset.[]{data-label="table:isvm for ijcnn"}

Conclusions {#sec:conclusion}
===========

We systematically investigate fast subgradient-based sampling methods for Bayesian max-margin models. To this end we first study general stochastic subgradient HMC sampling and several variants. We then apply the stochastic subgradient sampler to Bayesian linear SVMs and mixture of SVMs and further develop a doubly stochastic MCMC for the parametric mixtures of SVMs. However we note that our stochastic subgradient MCMC can also be used for other Bayesian models with non-differentiable posteriors, e.g, sparse models with a Laplace prior (See a detailed analysis in Appendix \[appsec:sparse model\]). Extensive empirical studies demonstrate the effectiveness of the stochastic subgradient MCMC methods on improving time efficiency while remaining a high accuracy of the samples.

We systematically investigate the subgradient MCMC methods for a large class of Bayesian models whose posteriors involve some non-differentiable terms in the log space, such as the sparse Bayesian models with a Laplace prior and the Bayesian SVMs with a piecewise linear hinge loss. We show that by using subgradients the common properties (e.g., volume preservation and detail balance) of an ordinary Hamiltonian Monte Carlo sampler still hold. We further propose to do stochastic subgradient MCMC for dealing with large-scale applications, with extensive empirical results on both low-dimensional and high-dimensional data. Our results demonstrate the effectiveness of the stochastic subgradient MCMC methods on improving the time efficiency and the high accuracy on drawing posterior samples.

In this paper, we have derived a class of the stochastic subgradient MCMC methods for the scalable sampling of the non-differentiable posteriors. This class can be seen as a generalization of the previous stochastic gradient LMC and HMC. The stochastic subgradient MCMC methods realized the stochastic sampling for the Bayesian models which has a non-differentiable posteriors including Bayesian max-margin learning and sparse learning models.

The theoretical analysis is investigated and the detailed balance of the subgradient HMC is shown after proving the several significant properties of the generalized Hamiltonian dynamics.

The method has been tested on the Bayesian max-margin learning and the sparse learning cases. The experimental results demonstrate the efficiency and the correctness of the stochastic subgradient MCMC methods. More over, the adaptive stepsize is also tested for the SSGLD method for a better mixing. The future work for the stochastic subgradient MCMC includes the theoretical analysis of the adaptive stepsize for the SSGLD.

### Acknowledgments {#acknowledgments .unnumbered}

Use unnumbered third level headings for the acknowledgments. All acknowledgments go at the end of the paper. Do not include acknowledgments in the anonymized submission, only in the final paper.

### References {#references .unnumbered}

Supplementary Material {#supplementary-material .unnumbered}
======================

Polynomial smooth of Hamiltonian Dynamics with Almost Differentiable Potential Energy {#appsec:polynomial construction}
=====================================================================================

For a given Hamiltonian dynamics $H_0$ with an almost differentiable potential energy function (we call it generalized Hamiltonian dynamics), a sequence of ordinary Hamiltonian dynamics $H_\epsilon$ can be constructed to approximate $H_{0}$. In this section, we construct the smoothness of the dynamics whose posterior gradient has at most a countable infinite number of non-differentiable points.

For any cuspidal point $q_{0}\in D_{0}$, we give a smoothness construction of the potential energy $U_{\epsilon}(q)$ using the polynomial smooth: $$\label{eqn:polynomial smooth}
U_{\epsilon}(q)= \left\{ \begin{array}{ll}
   U_0(q), & \|q-q_{0}\| > \epsilon\\
   \mathcal{P}(q), & \|q-q_{0}\| \leq \epsilon \end{array} \right.$$ where the polynomial function $\mathcal{P}$ satisfies the following conditions: $$\begin{aligned}
\label{eqn:polynomial construction}
\left\{ \begin{array}{ll}
U_{\epsilon}(q) = U_{0}(q),q=q_{0}\pm\epsilon,\\
\frac{\mathrm{d}U_{\epsilon}(q)}{\mathrm{d}q} = \mathcal{G}_{q}U_{0}(q),q=q_{0}\pm\epsilon.
\end{array} \right.\end{aligned}$$ For a given $U_0(q)$, the multi-dimensional polynomial interpolation ensures the existence of the $U_{\epsilon}(q)$ as defined in Eqn. (\[eqn:polynomial construction\]).

If $\epsilon>0$, $H_{\epsilon}$ is ordinary Hamiltonian dynamics. Specially, $H_\epsilon$ would be generalized Hamiltonian dynamics if $\epsilon=0$. A two-dimensional approximation demo can be seen in Fig. \[fig:approximation demo\].

From the construction of the polynomial smooth in Eqn. (\[eqn:polynomial smooth\]), it can simply be derived when $\epsilon\rightarrow0$, $$U_{\epsilon}\rightarrow U_{0}~,~\frac{\partial}{\partial q}U_{\epsilon} \overset{a.e.}{\rightarrow}  \mathcal{G}_{q}U_{0}$$

Proof of the Proposition 1 {#Proof of the volume preservation}
==========================

We prove the Proposition \[pro:volume preservation\] by using the phase flow to describe the transformation of the volume.

Let $g^t$ be the *phase flow* which is the one parameter group of the transformations of the phase space $$g^{t}:(p(0),q(0))\xrightarrow{f} (p(t),qï¼ˆt))$$ where $p(\cdot)$ and $q(\cdot)$ are the solutions of the generalized Hamiltonian equation Eqn. (\[eqn:subgradient Hamiltonian function\]). The phase flow the corresponding region $D_0$ to $g^{t}D_{0}=D_{t}$.

Using the Hamiltonian equation Eqn. (\[eqn:subgradient Hamiltonian function\]), the phase flow $g_{\epsilon}^{t}$ can be written as $$g^{t}(p,q)=(p,q)+\left(-\mathcal{G}_{q}U(q),M^{-1}q\right)*t+O(t^2),~(t\rightarrow 0)$$

After introducing the phase flow, we can write the phase-space volume of the certain generalized Hamiltonian dynamics after the flow as $$\label{eqn:flow change}
\rho(t)=\int_{D_{0}}\det\frac{\partial g^{t}}{\partial(p,q)}\mathrm{d}p\mathrm{d}q=(\int_{D_{0}\setminus D^{*}}+\int_{D_{*}})  \det\frac{\partial g^{t}}{\partial(p,q)}\mathrm{d}p\mathrm{d}q$$ where $D^{*}$ is the zero-mass set where the potential energy is non-differentiable.

Using the existing result in [@arnol1989mathematical], $$\begin{aligned}
\int_{D_{0}\setminus D^{*}}\det\frac{\partial g^{t}}{\partial(p,q)}\mathrm{d}p\mathrm{d}q &=& \int_{D_{0}\setminus D^{*}} 1+O(t^2)\mathrm{d}p\mathrm{d}q \\
&=& \int_{D_{0}} 1+O(t^2)\mathrm{d}p\mathrm{d}q = \rho(0),~~t\rightarrow0.\end{aligned}$$

To prove the volume preservation, i.e., Eqn. \[eqn:flow change\] equals $\rho(0)$, we only need show $$\label{eqn:zero mass volume}
\int_{D^{*}}\det\frac{\partial g^{t}}{\partial(p,q)}\mathrm{d}p\mathrm{d}q=0,$$ and hereinafter we finish the proof of Proposition. \[pro:volume preservation\] by showing Eqn. \[eqn:zero mass volume\].

We can write the determinant of the Jacobian matrix in Eq as $$\det\frac{\partial g^{t}}{\partial(p,q)}=\begin{vmatrix}
\nabla_{p}{G}_{q}U(q) & M^{-1} \\
{G}_{q}^{2} U(q) & \nabla_{q}M^{-1}p.
\end{vmatrix}$$

The element in the determinant are all finite except the ${G}_{q}^{2} U(q)$ which is a delta function. Then the Eqn. \[eqn:zero mass volume\] can be written as $$\begin{aligned}
&& \int_{D^{*}}\det\frac{\partial g^{t}}{\partial(p,q)}\mathrm{d}p\mathrm{d}q \\
&=& \int_{D^{*}}\begin{vmatrix}
\nabla_{p}{G}_{q}U(q) & M^{-1} \\
{G}_{q}^{2} U(q) & \nabla_{q}M^{-1}p
\end{vmatrix}\mathrm{d}p\mathrm{d}q \\
&=& \int_{D^{*}} {p}{G}_{q}U(q)\cdot-M^{-1}\cdot {G}_{q}^{2} U(q) \mathrm{d}p\mathrm{d}q \\
&<& C_{1}\int_{D^{*}}\end{aligned}$$

Proof of Theorem 2 {#Proof of the detailed balance}
==================

With the proof of the properties of the generalized Hamiltonian dynamics, reversibility, conservation of the Hamiltonian and the volume preservation, the detailed balance of the HMC method can be similarly proved as in the proof for the Hamiltonian dynamics [@neal2011mcmc].

Another significant issue for applying the subgradient HMC is the volume preservation for the commonly used discretization integrators, e.g., the Euler method and the leapfrog method. Similarity, the volume preservation and the reversibility holds in the subgradient HMC [@neal2011mcmc].

Generalization of the Synthetic Dataset {#appsec:synthetic dataset}
=======================================

Specifically, we generate 1,000 input observations from a uniform distribution $x_i \stackrel{i.i.d} \sim U(0,1)$ and the two-dimensional coefficient vector as $w \sim {\mathcal{N}}(0,\lambda^{-1} I)$, where $\lambda = 3$. Given the input observations and the coefficients, the labels are generated from a Bernoulli distribution $B(\alpha, 1-\alpha)$, where the parameter is $$\begin{aligned}
\nonumber \alpha = \frac{ P(y=1 | x_{i}, w) }{P(y=1 | x_{i}, w) + P(y=-1 | x_{i}, w)}. $$

Better Mixing of the Adaptive stepsizes for SSGLD {#adagrad}
=================================================

We give the empirical investigation of the adaptive stepsize for the SSGLD compared with the polynomial stepsizes used in the SGLD [@welling2011bayesian]. For the Bayesian linear SVM model on the Higgs dataset, it can be seen in the Figure \[fig:adagrad\] that adaptive stepsize has a better performance compared with the tuned polynomial decaying stepsizes. This improve might be due to the flexibility decaying at different feature dimension.

0.2in

![Performance of the SSGLD method with adaptive stepsizes on the Higgs dataset, including ([**Left**]{}) running mean plot; and ([**Right**]{}) the convergence curve of testing accuracy.[]{data-label="fig:meanplot"}](meanplot.eps){width="3in"}

0.2in

Derivations of Doubly Stochastic Sampling for Parametric Mixture of SVMs {#appsec:doubly stochastic}
========================================================================

In this section, we give the detailed derivations of the expected gradients of the parametric mixture of SVMs over the Gaussian component parameters $\gamma=(\mu,\Sigma)$. For the Gaussian mean, the expected gradient is $$\frac{\partial}{\partial \mu_{k}}\log q(\eta,\mu,\Sigma)=\sum_{i}p(z_{i}=k|x_{i},y_{i},\gamma,\eta) \Sigma_{k}^{-1}(x_{i}-\mu_{k}).$$ For the Gaussian covariance matrices, the expected gradient is $$\frac{\partial}{\partial \Sigma_{k}}\log q(\eta,\mu,\Sigma)=\sum_{i}p(z_{i}=k|x_{i},y_{i},\gamma,\eta) \frac{1}{2}\left(-\Sigma_{k}^{-1}+\Sigma_{k}^{-1}(x_{i}-\mu_{k})(x_{i}-\mu_{k})^{\top} \Sigma_{k}^{-1}\right)$$ Considering the PSD constrain of the covariance matrix $\Sigma_{k}$, we take the decomposition form $\Sigma_{k}=L_{k}^{\top}L_{k}$ and obtain the following gradients $$\frac{\partial}{\partial L_{k}}\log q(\eta,\mu,L) = \sum_{i}p(z_{i}=k|x_{i},y_{i},\gamma,\eta) \left(-L_{k}^{-\top}+2L_{k}^{-\top}(x_{i}-\mu_{k})(x_{i}-\mu_{k})^{\top }\right).$$

The probability of component assignment would be: $$\begin{aligned}
p(z_{i}=k|x_{i},y_{i},\gamma,\eta) &\propto& p(x_{i},y_{i}|z_{i}=k,\gamma,\eta)p(z_{i}=k) \\
&=& \mathcal{N}(x_{i}|\mu_{k},\Sigma_{k})\phi(y_{i}|\eta_{k},x_{i})\cdot p(z_{i}=k)\end{aligned}$$

Stochastic Subgradient MCMC for Bayesian Sparse Models {#appsec:sparse model}
======================================================

Another important application of stochastic subgradient HMC is the sparse Bayesian models with a non-differentiable Laplace prior. In this section, we use sparse Bayesian logistic regression as an example. We still consider the binary classification setting with inputs $x_i \in \mathbb{R}^n$ and class labels $y_{i}\in\{+1,-1\}$. Using a sigmoid function $\sigma(z)=\frac{1}{1+\exp(-z)}$, the posterior distribution of the model is: $$P(\eta|{\mathcal{D}})=P_{0}(\eta)\prod_{i}P(y_{i}|x_{i}, \eta)=P_{0}(\eta)\prod_{i}\sigma(y_{i}\eta^{\top}x_{i}). \nonumber$$ Here the prior of the vector $\eta$ is the Laplace distribution $P_{0}(\eta)=\exp\left(-\lambda^{-1}|\eta|_{1}\right)$ with which we can do the sparse feature selection. The posterior computing is also difficult due to the non-conjugacy between the Laplace prior and the sigmoid likelihood. In [@korattikara2014austerity], the reverse-jump MCMC method is used to do the stochastic sampling for this model with the Metropolis Hastings test. In RJ-MCMC, the MH iterations are finished by adding a binary vector to the corresponding feature vector $\eta$ and the RJ-MCMC proposal gives samples through a mixture of 3 types of moves (i.e., birth, death and update). This method has already been successfully applied in the sampling of Bayesian lasso [@chen2011bayesian]. For the stochastic subgradient MCMC methods, the subgradient of the log likelihood is $$\begin{aligned}
\nonumber \mathcal{G}_{\eta}\log P(y_{i}|x_{i}, \eta) = \sigma(y_{i}\eta^{\top}x_{i})y_{i}x_{i},\end{aligned}$$ and the subgradient of the log Laplace prior is just taken as $\mathcal{G}_{\eta}\log P_{0}(\eta)=-\lambda^{-1}\mathop{\textrm{sign}}(\eta)$. As for the stochastic subgradient MCMC, [@welling2011bayesian] ever give a description of SGLD to solve this problem using the subgradient information, but without careful justification. Here, we give a systematic study on both stochastic subgradient LMC and HMC. We first test our methods on the MiniBooNE dataset from the UCI machine learning repository [@korattikara2014austerity]. The MiniBooNE dataset contains $130,064$ data samples in a 50-dimensional feature space. We randomly choose $100,000$ samples as the training data and use the rest as the testing data. We compared our stochastic subgradient MCMC methods with the stochastic RJ-MCMC [@korattikara2014austerity] and the stochastic random walk Metropolis [@chib1995understanding]. We choose 0.01 as the variance parameter for the stochastic random walk Metropolis and the same variance setting with [@korattikara2014austerity] for RJ-MCMC. For SSGNHT, the diffusion factor $A$ is set as $1$, the stepsize is $0.0001\times t^{-0.6}$ and the step number is 30. We take 1,000 as the batchsize for the four stochastic sampling methods. We also use the adaptive stepsizes in the SSGLD method for better mixing. We start all the sampling methods with only the first dimension of the coefficient being 1 and all the others being 0.

0.2in

![Feature analysis of the sparse Bayesian logistic regression on the MiniBooNE dataset. **First row**: feature coefficients; and **Second row**: feature frequency histogram. []{data-label="fig:sparseboonefeaturehist"}](sparsefeature.eps){width="3.5in"}

-0.2in

0.15in

  Method          Selected feature indexes               
  -------- --------------------------------------- -- -- --
  SRWM      1, 4, 10, 11, 13, 20, 21, 27, 41, 42         
  SSGLD     1, 7, 13, 20, 21, 30, 31, 32, 37, 49         
  SSGNHT    1, 13, 20, 21, 26, 31, 32, 37, 44, 49        

  : Selected feature indexes of sparse Bayesian logistic regression by various samplers on the MiniBooNE dataset[]{data-label="table:feature"}

Figure \[fig:sparseboone\] presents the results of various methods with respect to the running time (in log-scale). For the stochastic random walk Metropolis and the stochastic RJ-MCMC, it takes a long time for them in the high dimension space to get a meaningful move, and the zigzag pattern of their convergence curves reveals the low sampling efficiency. For the stochastic subgradient MCMC methods, the accuracy of the sparse Bayesian logistic regression converges within a few steps. The accuracy of the RJ-MCMC method converges with about $0.85$ with $9$ selected features. Our stochastic subgradient MCMC methods get the ten times faster convergence compared with the stochastic random-walk Metropolis and the stochastic augmented RJ-MCMC method.

Again, we test the methods on the high dimensional Realsim dataset. We set the batchsize to 10 for all the stochastic sampling methods. For the RJ-MCMC method, the two sigma variances are set as the same value of 1, and for the stochastic random walk Metropolis, the variance parameter is also set as 1. For SSGNHT, the stepsize is taken as $0.01\times t^{-0.4}$ and the step number is set as 30. For SSGLD, the adaptive stepsize is used. Figure \[fig:sparserealsim\] presents the results of various methods with respect to the running time (in log-scale). We can see that the stochastic random walk Metropolis method has a low sampling efficiency, again due to the high dimensionality of the sample space. The stochastic RJ-MCMC suffers from the local minima, which are caused by introducing auxiliary binary variables [@korattikara2014austerity], and we select a relatively good convergence among the local minima to report. We analyze the features selected by the stochastic subgradient MCMC methods and the stochastic random-walk Metropolis. Figure \[fig:sparseboonefeaturehist\] shows the feature coefficient and the standardized feature frequency histogram. We also show the top 10 features with the largest absolute weights in Table \[table:feature\].

We can see that the two stochastic subgradient MCMC methods produce very similar feature ranks (e.g., 8 of the top-10 features are shared), while the similarity to the random-walk method is smaller (e.g., only 4 of the top-10 features are shared between SRWM and SSGLD).

Overall, in the high dimensional sparse learning case, our stochastic subgradient MCMC methods are faster than the chosen baseline methods by several orders of magnitude in terms of running time.

[^1]: http://www.csie.ntu.edu.tw/\~cjlin/libsvmtools/datasets/binary.html

[^2]: http://www.csie.ntu.edu.tw/\~cjlin/libsvmtools/datasets/binary.html
