Stochastic Quasi-Newton Langevin Monte Carlo

Introduction

Markov Chain Monte Carlo (MCMC) methods are one of the most important family of algorithms in Bayesian machine learning. These methods lie in the core of various applications such as the estimation of Bayesian predictive densities and Bayesian model selection. They also provide important advantages over optimization-based point estimation methods, such as having better predictive accuracy and being more robust to over-fitting (Chen et al., 2014; Ahn et al., 2015). Despite their well-known benefits, during the last decade, conventional MCMC methods have lost their charm since they are often criticized as being computationally very demanding. Indeed, classical approaches based on batch Metropolis Hastings would require passing over the whole data set at each iteration and the acceptance-rejection test makes the methods even more impractical for modern data sets.

Recently, alternative approaches have been proposed to scale-up MCMC inference to large-scale regime. An important attempt was made by Welling & Teh (2011), where the authors combined the ideas from Langevin Monte Carlo (LMC) (Rossky et al., 1978; Roberts & Stramer, 2002; Neal, 2010) and stochastic gradient descent (SGD) (Robbins & Monro, 1951; Kushner & Yin, 2003), and developed a scalable MCMC framework referred to as stochastic gradient Langevin dynamics (SGLD). Unlike conventional batch MCMC methods, SGLD uses subsamples of the data per iteration similar to SGD. With this manner, SGLD is able to scale up to large datasets while at the same time being a valid MCMC method that forms a Markov chain asymptotically sampling from the target density. Several extensions of SGLD have been proposed (Ahn et al., 2012; Patterson & Teh, 2013; Ahn et al., 2014; Chen et al., 2014; Ding et al., 2014; Ma et al., 2015; Chen et al., 2015; Shang et al., 2015; Li et al., 2016), which are coined under the term Stochastic Gradient MCMC (SG-MCMC).

One criticism that is often directed at SGLD is that it suffers from poor mixing rates when the target densities exhibit strong couplings and scale differences across dimensions. This problem is caused by the fact that SGLD is based on an isotropic Langevin diffusion, which implicitly assumes that different components of the latent variable are uncorrelated and have the same scale, a situation which is hardly encountered in practical applications.

In order to be able to generate samples from non-isotropic target densities in an efficient manner, SGLD has been extended in several ways. The common theme in these methods is to incorporate the geometry of the target density to the sampling algorithm. Towards this direction, a first attempt was made by (Ahn et al., 2012), where the authors proposed Stochastic Gradient Fisher Scoring (SGFS), a preconditioning schema that incorporates the curvature via Fisher scoring. Later, Patterson & Teh (2013) presented Stochastic Gradient Riemannian Langevin Dynamics (SGRLD), where they defined the sampler on the Riemannian manifold by borrowing ideas from (Girolami & Calderhead, 2011). SGRLD incorporates the local curvature information through the expected Fisher information matrix (FIM), which defines a natural Riemannian metric tensor for probability distributions. SGRLD has shown significant performance improvements over SGLD; however, its computational complexity is very intensive since it requires storing and inverting huge matrices, computing Cholesky factors, and expensive matrix products for generating each sample. It further requires computing the third order derivatives and the expected FIM to be analytically tractable, which limit the applicability of the method to a narrow variety of statistical models (Calderhead & Sustik, 2012).

In a very recent study, Li et al. (2016) proposed Preconditioned SGLD (PSGLD), that aims to handle the scale differences in the target density by making use of an adaptive preconditioning matrix that is chosen to be diagonal for computational purposes. Even though this method is computationally less intensive than SGRLD, it still cannot handle highly correlated target densities since the preconditioning matrix is forced to be diagonal. Besides, in order to reduce the computational burden, the authors discard a correction term in practical applications, which introduces permanent bias that does not vanish asymptotically.

In this study, we propose Hessian Approximated MCMC (HAMCMC), a novel SG-MCMC method that computes the local curvature of the target density in an accurate yet computationally efficient manner. Our method is built up on similar ideas from the Riemannian approaches; however, instead of using the expected FIM, we consider the local Hessian of the negative log posterior, whose expectation coincides with the expected FIM. Instead of computing the local Hessian and inverting it for each sample, we borrow ideas from Quasi-Newton optimization methods and directly approximate the inverse Hessian by using a limited history of samples and their stochastic gradients. By this means, HAMCMC is (i) computationally more efficient than SGRLD since its time and memory complexities are linear with the dimension of the latent variable, (ii) it can be applied to a wide variety of statistical models without requiring the expected FIM to be analytically tractable, and (iii) it is more powerful than diagonal preconditioning approaches such as PSGLD, since it uses dense approximations of the inverse Hessian, hence is able to deal with correlations along with scale differences.

We provide rigorous theoretical analysis for HAMCMC, where we show that HAMCMC is asymptotically unbiased and consistent with posterior expectations. We evaluate HAMCMC on a synthetic and two challenging real datasets. In particular, we apply HAMCMC on two different matrix factorization problems and we evaluate its performance on a speech enhancement and a distributed link prediction application. Our experiments demonstrate that HAMCMC achieves fast convergence rates similar to SGRLD while at the same time having low computational requirements similar to SGLD and PSGLD.

We note that SGLD has also been extended by making use of higher order dynamics (Chen et al., 2014; Ding et al., 2014; Ma et al., 2015; Shang et al., 2015). These extensions are rather orthogonal to our contributions and HAMCMC can be easily extended in the same directions.

Preliminaries

We define an unbiased estimator of $U(\theta)$ as follows:

By combining ideas from LMC and SGD, Welling & Teh (2011) presented SGLD that asymptotically generates a sample $\theta_{t}$ from the posterior distribution by iteratively applying the following update equation:

where $\epsilon_{t}$ is the step-size, and $\eta_{t}$ is Gaussian noise: $\eta_{t}\sim{\cal N}(0,2\epsilon_{t}I)$, and ${I}$ stands for the identity matrix.

This method can be seen as the Euler discretization of the Langevin dynamics that is described by the following stochastic differential equation (SDE):

2 The L-BFGS algorithm

In this study, we consider a limited-memory Quasi-Newton (QN) method, namely the L-BFGS algorithm (Nocedal & Wright, 2006) that iteratively applies the following equation in order to find the maximum a-posteriori estimate:

where $H_{t}$ is an approximation to the inverse Hessian at $\theta_{t-1}$. The L-BFGS algorithm directly approximates the inverse of the Hessian by using the $M$ most recent values of the past iterates. At each iteration, the inverse Hessian is approximated by applying the following recursion: (using $H_{t}=H_{t}^{M-1}$)

where $s_{t}\triangleq\theta_{t}-\theta_{t-1}$, $y_{t}\triangleq\nabla{U}(\theta_{t})-\nabla{U}(\theta_{t-1})$, and $\tau=t-M+m$ and the initial approximation is chosen as $H_{t}^{1}=\gamma I$ for some $\gamma>0$. The matrix-vector product $H_{t}\nabla U(\theta_{t-1})$ is often implemented either by using the two-loop recursion (Nocedal & Wright, 2006) or by using the compact form (Byrd et al., 1994), which both have linear time complexity ${\cal O}(MD)$. Besides, since L-BFGS needs to store only the latest $M-1$ values of $s_{t}$ and $y_{t}$, the space complexity is also ${\cal O}(MD)$, as opposed to classical QN methods that have quadratic memory complexity.

Surprisingly, QN methods have not attracted much attention from the MCMC literature. Zhang & Sutton (2011) combined classical Hamiltonian Monte Carlo with L-BFGS for achieving better mixing rates for small- and medium-sized problems. This method considers a batch scenario with full gradient computations which are then appended with costly Metropolis acceptance steps. Recently, Dahlin et al. (2015) developed a particle Metropolis-Hastings schema based on (Zhang & Sutton, 2011).

Stochastic Quasi-Newton LMC

Recent studies have shown that Langevin and Hamiltonian Monte Carlo methods have strong connections with optimization techniques. Therefore, one might hope for developing stronger MCMC algorithms by borrowing ideas from the optimization literature (Qi & Minka, 2002; Zhang & Sutton, 2011; Bui-Thanh & Ghattas, 2012; Pereyra, 2013; Chaari et al., 2015; Bubeck et al., 2015; Li et al., 2016). However, incorporating ideas from optimization methods to an MCMC framework requires careful design and straightforward extensions often do not yield proper algorithms (Chen et al., 2014; Ma et al., 2015). In this section, we will first show that a naïve way of developing an SG-MCMC algorithm based on L-BFGS would not target the correct distribution. Afterwards, we will present our proposed algorithm HAMCMC and show that it targets the correct distribution and it is asymptotically consistent with posterior expectations.

A direct way of using SQN ideas in an LMC framework would be achieved by considering $H_{t}$ as a preconditioning matrix in an SGLD context (Welling & Teh, 2011; Ahn et al., 2012) by combining Eq. 3 and Eq. 5, which would yield the following update equation:

where $H_{t}(\cdot)$ is the approximate inverse Hessian computed via L-BFGS, $M$ denotes the memory size in L-BFGS and \eta_{t}\sim{\cal N}\bigl{(}0,2\epsilon_{t}H_{t}(\theta_{t-M:t-1})\bigr{)}. Here, we use a slightly different notation and denote the L-BFGS approximation via $H_{t}(\theta_{t-M:t-1})$ instead of $H_{t}$ in order to explicitly illustrate the samples that are used in the L-BFGS calculations.

Despite the fact that such an approach would introduce several challenges, which will be described in the next section, for now let us assume computing Eq. 7 is feasible. Even so, this approach does not result in a proper MCMC schema.

The Markov chain described in Eq. 7 does not target the posterior distribution $p(\theta|x)$ unless $\sum_{j=1}^{D}\frac{\partial}{\partial\theta_{j}}H_{ij}(\theta)=0$ for all $i\in\{1,\dots,D\}$.

Eq. 7 can be formulated as a discretization of a continuous-time diffusion that is given as follows:

Due to the hypothesis, Eq. 8 clearly does not target the posterior distribution, therefore Eq. 7 does not target the posterior distribution either.∎

In fact, this result is neither surprising nor new. It has been shown that in order to ensure targeting the correct distribution, we should instead consider the following SDE (Girolami & Calderhead, 2011; Xifara et al., 2014; Patterson & Teh, 2013; Ma et al., 2015):

where $\bar{g}_{\Omega}(\theta)\triangleq(1/N_{\Omega})\sum_{n\in\Omega}\nabla\log p(x_{n}|\theta)$, the symbols $\cdot/\cdot$ and $\cdot\circ\cdot$ respectively denote element-wise division and multiplication, and $\alpha\in$ and $\lambda\in\mathds{R}_{+}$ are the parameters to be tuned.

In classical Metropolis-Hastings settings, where Langevin dynamics is only used in the proposal distribution, the correction term $\Gamma(\theta_{t})$ can be discarded without violating the convergence properties thanks to the acceptance-rejection step (Qi & Minka, 2002; Girolami & Calderhead, 2011; Zhang & Sutton, 2011; Calderhead & Sustik, 2012). However, such an approach would result in slower convergence (Girolami & Calderhead, 2011). On the other hand, in SG-MCMC settings which do not include an acceptance-rejection step, this term can be problematic and should be handled with special attention.

In (Ahn et al., 2012), the authors discarded the correction term in a heuristic manner. Recently, Li et al. (2016) showed that discarding the correction term induces permanent bias which deprives the methods of their asymptotic consistency and unbiasedness. In the sequel, we will show that the computation of $\Gamma(\theta_{t})$ can be avoided without inducing permanent bias by using a special construction that exploits the limited memory structure of L-BFGS.

2 Hessian Approximated MCMC

In this section, we present our proposed method HAMCMC. HAMCMC applies the following update rules for generating samples from the posterior distribution:

where ${\theta}_{a:b}^{\neg c}\equiv(\theta_{a:b}\setminus\theta_{c})$, \eta_{t}\sim{\cal N}\bigl{(}0,2\epsilon_{t}H_{t}(\cdot)\bigr{)}, and we assume that the initial $2M+1$ samples, i.e. $\theta_{0},\dots,\theta_{2M}$ are already provided.

An important property of HAMCMC is that the correction term $\Gamma(\theta)$ vanishes due to the construction of our algorithm, which will be formally proven in the next section. Informally, since $H_{t}$ is independent of $\theta_{t-M}$ conditioned on ${\theta}_{t-2M+1:t-1}^{\neg(t-M)}$, the change in $\theta_{t-M}$ does not affect $H_{t}$ and therefore all the partial derivatives in Eq. 9 become zero. This property saves us from expensive higher-order derivative computations or inducing permanent bias due to the negligence of the correction term. On the other hand, geometrically, our approach corresponds to choosing a flat Riemannian manifold specified by the metric tensor $H_{t}$, whose curvature is constant (i.e. $\Gamma(\theta)=0$).

In order to have a valid sampler, we are required to have positive definite $H_{t}(\cdot)$. However, even if $U(\theta)$ was strongly convex, due to data subsampling L-BFGS is no longer guaranteed to produce positive definite approximations. One way to address this problem is to use variance reduction techniques as presented in (Byrd et al., 2014); however, such an approach would introduce further computational burden. In this study, we follow a simpler approach and address this problem by approximating the inverse Hessian in a trust region that is parametrized by $\lambda$, i.e. we use L-BFGS to approximate $(\nabla^{2}{U}(\theta)+\lambda I)^{-1}$. For large enough $\lambda$ this approach will ensure positive definiteness. The L-BFGS updates for this problem can be obtained by simply modifying the gradient differences as: $y_{t}\leftarrow y_{t}+\lambda s_{t}$ (Schraudolph et al., 2007).

The final challenge is to generate $\eta_{t}$ whose covariance is a scaled version of the L-BFGS output. This requires to compute the square root of $H_{t}$, i.e. $S_{t}S_{t}^{\top}=H_{t}$, so that we can generate $\eta_{t}$ as $\sqrt{2\epsilon_{t}}S_{t}z_{t}$, where $z_{t}\sim{\cal N}(0,I)$. Fortunately, we can directly compute the product $S_{t}z_{t}$ within the L-BFGS framework by using the following recursion (Brodlie et al., 1973; Zhang & Sutton, 2011):

where $S_{t}\equiv S_{t}^{M-1}$ and $B_{t}^{m}=(H_{t}^{m})^{-1}$. By this approach, $\eta_{t}$ can be computed in ${\cal O}(M^{2}D)$ time, which is still linear in $D$. We illustrate HAMCMC in Algorithm 1.

Note that large $M$ would result in a large gap between two consecutive samples and therefore might require larger $\lambda$ for ensuring positive definite L-BFGS approximations. Such a circumstance would be undesired since $H_{t}$ would get closer to the identity matrix and therefore result in HAMCMC being similar to SGLD. In our computational studies, we have observed that choosing $M$ between $2$ and $5$ is often adequate in several applications.

3 Convergence Analysis

We are interested in approximating the posterior expectations by using sample averages:

where $\pi(\theta)\triangleq p(\theta|x)$, $W_{T}\triangleq\sum_{t=1}^{T}\epsilon_{t}$, and $\theta_{t}$ denotes the samples generated by HAMCMC. In this section we will analyze the asymptotic properties of this estimator. We consider a theoretical framework similar to the one presented in (Chen et al., 2015). In particular, we are interested in analyzing the bias \bigl{|}\mathds{E}[\hat{h}]-\bar{h}\bigr{|} and the mean squared-error \mathds{E}\bigl{[}(\hat{h}-\bar{h})^{2}\bigr{]} of our estimator.

A useful way of analyzing dynamics-based MCMC approaches is to first analyze the underlying continuous dynamics, then consider the MCMC schema as a discrete-time approximation (Roberts & Stramer, 2002; Sato & Nakagawa, 2014; Chen et al., 2015). However, this approach is not directly applicable in our case. Inspired by (Roberts & Gilks, 1994) and (Zhang & Sutton, 2011), we convert Eq. 10 to a first order Markov chain that is defined on the product space $\mathds{R}^{D(2M-1)}$ such that $\Theta_{t}\equiv\{\theta_{t-2M+2},\dots,\theta_{t}\}$. In this way, we can show that HAMCMC uses a transition kernel which modifies only one element of $\Theta_{t}$ per iteration and rearranges the samples that will be used in the L-BFGS computations. Therefore, the overall HAMCMC kernel can be interpreted as a composition of $M$ different preconditioned SGLD kernels whose volatilities are independent of their current states. By considering this property and assuming certain conditions that are provided in the Supplement hold, we can bound the bias and MSE as shown in Theorem 1.

\bigl{|}\mathds{E}[\hat{h}]-\bar{h}\bigr{|}={\cal O}\bigl{(}\frac{1}{L_{K}}+\frac{Y_{K}}{L_{K}}\bigr{)}

\mathds{E}\bigl{[}(\hat{h}-\bar{h})^{2}\bigr{]}={\cal O}\Bigl{(}\sum\limits_{k=1}^{K}\frac{\epsilon_{kM}^{2}}{L_{K}^{2}}\mathds{E}\|\Delta V_{k^{\star}}\|^{2}+\frac{1}{L_{K}}+\frac{Y_{K}^{2}}{L_{K}^{2}}\Bigr{)}

where $h$ is smooth and $\Delta V_{k^{\star}}=\Delta V_{m_{k}^{*}+kM}$, where $m_{k}^{*}=\operatorname*{arg\max}_{1\leq m\leq M}\mathds{E}\|\Delta V_{m+kM}\|^{2}$, and $\|\Delta V_{t}\|$ denotes the operator norm.

The proof is given in the Supplement. The theorem implies that as $T$ goes to infinity, the bias and the MSE of the estimator vanish. Therefore, HAMCMC is asymptotically unbiased and consistent.

Experiments

We conduct our first set of experiments on synthetic data where we consider a rather simple Gaussian model whose posterior distribution is analytically available. The model is given as follows:

for all $n$. Here, we assume $\{a_{n}\}_{n=1}^{N}$ and $\sigma_{x}^{2}$ are known and we aim to draw samples from the posterior distribution $p(\theta|x)$. We will compare the performance of HAMCMC against SGLD, PSGLD, and SGRLD. In these experiments, we determine the step size as $\epsilon_{t}=(a_{\epsilon}/t)^{0.51}$. In all our experiments, for each method, we have tried several values for the hyper-parameters with the same rigor and we report the best results. We also report all of the hyper-parameters used in the experiments (including Sections 4.2-4.3) in the Supplement. These experiments are conducted on a standard laptop computer with $2.5$GHz Quad-core Intel Core i$7$ CPU and all the methods are implemented in Matlab.

In the first experiment, we first generate the true $\theta$ and the observations $x$ by using the generative model given in Eq. 13. We set $D=2$ for visualization purposes and $\sigma_{x}^{2}=10$, then we randomly generate the vectors $\{a_{n}\}_{n}$ in such a way that the posterior distribution will be correlated. Then we generate $T=20000$ samples by using the SG-MCMC methods, where the size of the data subsamples is selected as $N_{\Omega}=N/100$. We discard the first half of the samples as burn-in. Figure 1 shows the typical results of this experiment. We can clearly observe that SGLD cannot explore the posterior efficiently and gets stuck around the mode of the distribution. PSGLD captures the shape of the distribution in a more accurate way than SGLD since it is able to take the scale differences into account; however, it still cannot explore lower probability areas efficiently. Besides, we can see that SGRLD performs much better than SGLD and PSGLD. This is due to the fact that it uses the inverse expected FIM of the full problem for generating each sample, where this matrix is simply $(\sum_{n}a_{n}a_{n}^{\top}/\sigma_{x}^{2}+I)^{-1}$ for this model. Therefore, SGRLD requires much heavier computations as it needs to store a $D\times D$ matrix, computes the Cholesky factors of this matrix, and performs matrix-vector products at each iteration. Furthermore, in more realistic probabilistic models, the expected FIM almost always depends on $\theta$; hence, the computation burden of SGRLD is further increased by the computation of the inverse of the $D\times D$ matrix at each iteration. On the other hand, Fig. 1 shows that HAMCMC takes the best of both worlds: it is able to achieve the quality of SGRLD since it makes use of dense approximations to the inverse Hessian and it has low computational requirements similar to SGLD and PSGLD. This observation becomes more clear in our next experiment.

In our second experiment, we again generate the true $\theta$ and the observations by using Eq. 13. Then, we approximate the posterior mean $\bar{\theta}=\int\theta p(\theta|x)d\theta$ by using the aforementioned MCMC methods. For comparison, we monitor $\|\bar{\theta}-\hat{\theta}\|^{2}$ where $\hat{\theta}$ is the estimate that is obtained from MCMC. Figure 2 shows the results of this experiment for $D=10$ and $D=100$. In the left column, we show the error as a function of iterations. These results show that the convergence rate of SGRLD (with respect to iterations) is much faster than SGLD and PSGLD. As expected, HAMCMC yields better results than SGLD and PSGLD even in the simplest configuration ($M=2$), and it performs very similarly to SGRLD as we set $M=3$. The performance difference becomes more prominent as we increase $D$.

An important advantage of HAMCMC is revealed when we take into account the total running time of these methods. These results are depicted in the right column of Fig. 2, where we show the error as a function of wall-clock time The inverse of the FIM for this model can be precomputed. However, this will not be the case in more general scenarios. Therefore, for fair comparison we perform the matrix inversion and Cholesky decomposition operations at each iteration.. We can observe that, since SGLD, PSGLD, and HAMCMC have similar computational complexities, the shapes of their corresponding plots do not differ significantly. However, SGRLD appears to be much slower than the other methods as we increase $D$.

2 Alpha-Stable Matrix Factorization

In our second set of experiments, we consider a recently proposed probabilistic matrix factorization model, namely alpha-stable matrix factorization ($\alpha$MF), that is given as follows (Şimşekli et al., 2015):

where $X\in\mathds{C}^{I\times J}$ is the observed matrix with complex entries, and $W\in\mathds{R}_{+}^{I\times K}$ and $H\in\mathds{R}_{+}^{K\times J}$ are the latent factors. Here ${\cal GG}$ denotes the generalized gamma distribution (Stacy, 1962) and ${\cal S}\alpha{\cal S}_{c}$ denotes the complex symmetric $\alpha$-stable distribution (Samoradnitsky & Taqqu, 1994). Stable distributions are heavy-tailed distributions and they are the limiting distributions in the generalized central limit theorem. They are often characterized by four parameters: ${\cal S}(\alpha,\beta,\sigma,\mu)$, where $\alpha\in(0,2]$ is the characteristic exponent, $\beta\in$ is the skewness parameter, $\sigma\in(0,\infty)$ is the dispersion parameter, and $\mu\in(-\infty,\infty)$ is the location parameter. The distribution is called symmetric (${\cal S}\alpha{\cal S}$) if $\beta=0$. The location parameter is also chosen as in $\alpha$MF.

The probability density function of the stable distributions cannot be written in closed-form except for certain special cases; however it is easy to draw samples from them. Therefore, MCMC has become the de facto inference tool for $\alpha$-stable models like $\alpha$MF (Godsill & Kuruoglu, 1999). By using data augmentation and the product property of stable distributions, Şimşekli et al. (2015) proposed a partially collapsed Gibbs sampler for making inference in this model. The performance of this approach was then illustrated on a speech enhancement problem.

In this section, we derive a new inference algorithm for $\alpha$MF that is based on SG-MCMC. By using the product property of stable models (Kuruoglu, 1999), we express the model as conditionally Gaussian, given as follows:

where ${\cal N}_{c}$ denotes the complex Gaussian distribution. Here, by assuming $\alpha$ is already provided, we propose a Metropolis and SG-MCMC within Gibbs approach, where we generate samples from the full conditional distributions of the latent variables in such a way that we use a Metropolis step for sampling $\Phi$ where the proposal is the prior and we use SG-MCMC for sampling $W$ and $H$ from their full conditionals. Since all the elements of $W$ and $H$ must be non-negative, we apply mirroring at the end of each iteration where we replace negative elements with their absolute values (Patterson & Teh, 2013).

We evaluate HAMCMC on the same speech enhancement application described in (Şimşekli et al., 2015) by using the same settings. The aim here is to recover clean speech signals that are corrupted by different types of real-life noise signals (Hu & Loizou, 2007). We use the Matlab code that can be downloaded from the authors’ websites in order to generate the same experimental setting. We first train $W$ on a clean speech corpus, then we fix $W$ at denoising and sample the rest of the variables. For each test signal (out of $80$) we have $I=256$ and $J=140$. The rank is set $K=105$ where the first $100$ columns of $W$ is reserved for clean speech and for speech $\alpha=1.2$ and for noise $\alpha=1.89$. Note that the model used for denoising is slightly different from Eq. 14 and we refer the readers to the original paper for further details of the experimental setup.

In this experiment, we compare HAMCMC against SGLD and PSGLD. SGRLD cannot be applied to this model since the expected FIM is not tractable. For each method, we set $N_{\Omega}=IJ/10$ and we generate $2500$ samples where we compute the Monte Carlo averages by using the last $200$ samples. For HAMCMC, we set $M=5$. Fig. 3 shows the denoising results on the test set in terms of signal-to-distortion ratio (SDR) when the mixture signal-to-noise ratio is $5$dB. The results show that SGLD and PSGLD perform similarly, where the quality of the outputs of PSGLD is slightly higher than the ones obtained from SGLD. We can observe that HAMCMC provides a significant performance improvement over SGLD and PSGLD, thanks to the usage of local curvature information.

3 Distributed Matrix Factorization

In our final set of experiments, we evaluate HAMCMC on a distributed matrix factorization problem, where we consider the following probabilistic model: $W_{ik}\sim{\cal N}(0,\sigma_{w}^{2})$, $H_{kj}\sim{\cal N}(0,\sigma_{h}^{2})$, X_{ij}|\cdot\sim{\cal N}\bigl{(}\sum\nolimits_{k}W_{ik}H_{kj},\sigma^{2}_{x}\bigr{)}. Here, $W\in\mathds{R}^{I\times K}$, $H\in\mathds{R}^{K\times J}$, and $X\in\mathds{R}^{I\times J}$. This model is similar to (Salakhutdinov & Mnih, 2008) and is often used in distributed matrix factorization problems (Gemulla et al., 2011; Şimşekli et al., 2015b).

Inspired by Distributed SGLD (DSGLD) for matrix factorization (Ahn et al., 2015) and Parallel SGLD (Şimşekli et al., 2015a), we partition $X$ into disjoint subsets where each subset is further partitioned into disjoint blocks. This schema is illustrated in Fig. 4, where $X$ is partitioned into $4$ disjoint subsets $\Omega^{(1)},\dots,\Omega^{(4)}$ and each subset contains $4$ different blocks that will be distributed among different computation nodes. By using this approach, we extend PSGLD and HAMCMC to the distributed setting similarly to DSGLD, where at each iteration the data subsample $\Omega_{t}$ is selected as $\Omega^{(i)}$ with probability $|\Omega^{(i)}|/\sum_{j}|\Omega^{(j)}|$. At the end of each iteration the algorithms transfer a small block of $H$ to a neighboring node depending on $\Omega_{t+1}$. PSGLD and HAMCMC further need to communicate other variables that are required for computing $H_{t}$, where the communication complexity is still linear with $D$. Finally, HAMCMC requires to compute six vector dot products at each iteration that involves all the computation nodes, such as $s_{t}^{\top}y_{t}$. These computations can be easily implemented by using a reduce operation which has a communication complexity of ${\cal O}(M^{2})$. By using this approach, the methods have the same theoretical properties except that the stochastic gradients would be expected to have larger variances.

For this experiment, we have implemented all the algorithms in C by using a low-level message passing protocol, namely the OpenMPI library. In order to keep the implementation simple, for HAMCMC we set $M=2$. We conduct our experiments on a small-sized cluster with $4$ interconnected computers each of them with $4$ Intel Xeon $2.93$GHz CPUs and $192$ GB of memory.

We compare distributed HAMCMC against DSGLD and distributed PSGLD on a large movie ratings dataset, namely the MovieLens $1$M (grouplens.org). This dataset contains about $1$ million ratings applied to $I=3883$ movies by $J=6040$ users, resulting in a sparse observed matrix $X$ with $4.3\%$ non-zero entries. We randomly select $10\%$ of the data as the test set and partition the rest of the data as illustrated in Fig. 4. The rank of the factorization is chosen as $K=10$. We set $\sigma_{w}^{2}=\sigma_{h}^{2}=\sigma_{x}^{2}=1$. In these experiments, we use a constant step size for each method and discard the first $50$ samples as burn-in. Note that for constant step size, we no longer have asymptotic unbiasedness; however, the bias and MSE are still bounded.

Fig. 5 shows the comparison of the algorithms in terms of the root mean squared-error (RMSE) that are obtained on the test set. We can observe that DSGLD and PSGLD have similar converges rates. On the other hand, HAMCMC converges faster than these methods while having almost the same computational requirements.

Conclusion

We presented HAMCMC, a novel SG-MCMC algorithm that achieves high convergence rates by taking the local geometry into account via the local Hessian that is efficiently computed by the L-BFGS algorithm. HAMCMC is both efficient and powerful since it uses dense approximations of the inverse Hessian while having linear time and memory complexity. We provided formal theoretical analysis and showed that HAMCMC is asymptotically consistent. Our experiments showed that HAMCMC achieves better convergence rates than the state-of-the-art while keeping the computational cost almost the same. As a next step, we will explore the use of HAMCMC in model selection applications (Şimşekli et al., 2016).

Acknowledgments

The authors would like to thank to Ş. İlker Birbil and Figen Öztoprak for fruitful discussions on Quasi-Newton methods. The authors would also like to thank to Charles Sutton for his explanations on the Ensemble Chain Adaptation framework. This work is partly supported by the French National Research Agency (ANR) as a part of the EDISON 3D project (ANR-13-CORD-0008-02). A. Taylan Cemgil is supported by TÜBİTAK 113M492 (Pavera) and Boğaziçi University BAP 10360-15A01P1.

References