Unifying the Stochastic Spectral Descent for Restricted Boltzmann Machines with Bernoulli or Gaussian Inputs

TL;DR

This paper extends stochastic spectral descent (SSD) to Gaussian RBMs, enabling efficient training on continuous data by deriving new bounds and demonstrating improved convergence over traditional stochastic gradient descent.

Contribution

The paper generalizes SSD to Gaussian RBMs, providing theoretical bounds and empirical evidence of faster convergence for continuous data modeling.

Findings

SSD outperforms SGD in convergence speed

Derived upper bounds for the partition function using Schatten-infinity norm

Extended SSD applicability from binary to Gaussian RBMs

Abstract

Stochastic gradient descent based algorithms are typically used as the general optimization tools for most deep learning models. A Restricted Boltzmann Machine (RBM) is a probabilistic generative model that can be stacked to construct deep architectures. For RBM with Bernoulli inputs, non-Euclidean algorithm such as stochastic spectral descent (SSD) has been specifically designed to speed up the convergence with improved use of the gradient estimation by sampling methods. However, the existing algorithm and corresponding theoretical justification depend on the assumption that the possible configurations of inputs are finite, like binary variables. The purpose of this paper is to generalize SSD for Gaussian RBM being capable of mod- eling continuous data, regardless of the previous assumption. We propose the gradient descent methods in non-Euclidean space of parameters, via de- riving the upper bounds of logarithmic partition function for RBMs based on Schatten-infinity norm. We empirically show that the advantage and improvement of SSD over stochastic gradient descent (SGD).

Figures

Tables

Table 1. Table 1: Running time per 1k iterations

	SSD	SGD	Nesterov SGD
GRBM	67s	54s	58s
BRBM	48s	44s	45.5s

Equations

$$E(\mathbf{v},\mathbf{h};\mathbf{\theta })=-{\mathbf{v}}^{\top }W\mathbf{h}-{\mathbf{v}}^{\top }\mathbf{b}-{\mathbf{h}}^{\top }\mathbf{a}$$

$$p(\mathbf{v};\mathbf{\theta })$$

$$Z(\mathbf{\theta })$$

$$\mathcal{L}(\mathbf{\theta })=-\frac{1}{N}\sum _{n=1}^N\log p({\mathbf{v}}_n;\mathbf{\theta })=f(\mathbf{\theta })+g(\mathbf{\theta })$$

$${\nabla }_W\mathcal{L}(\mathbf{\theta })={\mathbb{E}}_p[\mathbf{v}{\mathbf{h}}^{\top }]-{\mathbb{E}}_e[\mathbf{v}{\mathbf{h}}^{\top }]$$

$$\mathcal{L}(W+\Delta W)\le \mathcal{L}(W)+⟨{\nabla }_W\mathcal{L},\Delta W⟩+\frac{N_v{N}_h\Vert \Delta W{\Vert }_{S_{\infty }}^2}{2},$$

$$E_G(\mathbf{v},\mathbf{h};\mathbf{\theta })=-{\mathbf{v}}^{\top }{C}^{-1}W\mathbf{h}+\frac{(\mathbf{v}-\mathbf{b}{)}^{\top }{C}^{-1}(\mathbf{v}-\mathbf{b})}{2}-{\mathbf{h}}^{\top }\mathbf{a}$$

$$f(\mathbf{\theta })=\log \left(\int \sum _{\mathbf{h}}\exp \left(-E_G(\mathbf{v},\mathbf{h};\mathbf{\theta })\right)\mathrm{d}\mathbf{v}\right)$$

$$p(\mathbf{v}|\mathbf{h},\mathbf{\theta })$$

$$p(h_k=1|\mathbf{v},\mathbf{\theta })$$

$$f(\mathbf{a}+\Delta \mathbf{a})$$

$$f(\mathbf{b}+\Delta \mathbf{b})$$

$$f(\mathbf{\theta })=$$

$$f(\mathbf{\theta })=\log \left(\sum _{i=1}^{2^{N_h}}{\omega }_i\exp \left(\frac{1}{2}s(W{)}^{\top }s(W)\right)\right)$$

$$f(W+U)$$

$$+\left(\frac{1}{2}+\frac{3r(R,C,\mathbf{b}{)}^2}{4}\right)N_v{N}_h\Vert U{\Vert }_{S_{\infty }}^2$$

$$g(\mathbf{\theta })=$$

$$-\sum _{k=1}^{N_h}\log \left(1+\exp \left([{\mathbf{v}}_n^{\top }{C}^{-1}W+\mathbf{a}{]}_k\right)\right)).$$

$$\mathcal{L}({\mathbf{\theta }}_i+\Delta {\mathbf{\theta }}_i)\le \mathcal{L}({\mathbf{\theta }}_i)+⟨{\nabla }_{{\mathbf{\theta }}_i}\mathcal{L}({\mathbf{\theta }}_i),\Delta {\mathbf{\theta }}_i⟩+c_i\Vert \Delta {\mathbf{\theta }}_i{\Vert }_{\infty }^2$$

$$\frac{\partial \mathcal{L}}{\partial W}$$

$$\frac{\partial \mathcal{L}}{\partial C^{-1}}$$

$$\frac{\partial \mathcal{L}}{\partial \mathbf{b}}$$

$$\frac{\partial \mathcal{L}}{\partial \mathbf{a}}$$

$$\frac{\partial \mathcal{L}}{\partial C^{-1}}$$

$$\frac{\partial \mathcal{L}}{\partial \mathbf{b}}$$

$$\frac{\partial \mathcal{L}}{\partial \ln (C)}$$

$$\frac{\partial \mathcal{L}}{\partial \ln (C^{-1})}$$

$$g(\mathbf{\theta })$$

$$=-\frac{1}{N}\sum _{n=1}^N\log \left(\exp \left(-\frac{1}{2}({\mathbf{v}}_n-\mathbf{b}{)}^{\top }{C}^{-1}({\mathbf{v}}_n-\mathbf{b})\right)\prod _{k=1}^{N_h}\left(1+\exp \left({\mathbf{v}}_n^{\top }{C}^{-1}{W}_{\cdot ,k}+a_k\right)\right)\right)$$

$$=\frac{1}{N}\sum _{n=1}^N\left(\frac{1}{2}({\mathbf{v}}_n-\mathbf{b}{)}^{\top }{C}^{-1}({\mathbf{v}}_n-\mathbf{b})-\sum _{k=1}^{N_h}\log \left(1+\exp \left({\mathbf{v}}_n^{\top }{C}^{-1}{W}_{\cdot ,k}+a_k\right)\right)\right)$$

$$g(\mathbf{b}+\Delta \mathbf{b})$$

$$f(\mathbf{a}+\Delta \mathbf{a})$$

$$f(\mathbf{\theta })$$

$$=\log \left(\sum _{i=1}^{2^{N_h}}{\omega }_i\exp ({\mathbf{h}}_i^{\top }\mathbf{a})\right)$$

$${\omega }_i$$

$$=(2\pi {)}^{\frac{N_v}{2}}|C{|}^{\frac{1}{2}}\exp \left({\mathbf{b}}^{\top }{C}^{-1}W{\mathbf{h}}_i+\frac{(W{\mathbf{h}}_i{)}^{\top }{C}^{-1}W{\mathbf{h}}_i}{2}\right)$$

Full text

1 Introduction

A crucial challenge to many deep learning models involves optimizing non-convex objective functions over continuous high dimensional matrix or tensor spaces, with large scale training dataset. The stochastic version of gradient descent or quasi-Newton method is ubiquitously applied to optimize the parameters of these models, intending to find an acceptable local optima. However, exploiting the inhomogeneous curvature of loss functions or the corresponding information in their gradients more effectively can greatly speed up the overall inference. [3] developed stochastic spectral descent (SSD) method by changing the space in which the gradient descent performs to better utilize the natural geometry implied in the Bernoulli Restricted Boltzmann Machine (RBM) [11] by operating in the $l_{\infty}$ and Schatten-$\infty$ (S-$\infty$) norms [1]. The advent of such customized decent method encourages more specific designed algorithms that allow for optimizing discrete deep modeling, [4, 2]. Figure 1 illustrates that SSD tends to find the shortcut path toward the local optimal.

The previous research of SSD is however restrictive and it is theoretically correct for discrete data with finite domain. A common trick to modeling real-valued data by RBM is normalizing each input variable to $[0,1]$ interval and considering it as Bernoulli distributed variable [12]. This approach has its own limitation as well since it usually requires the data input to be bounded. This means the generalization of RBM for Gaussian distributed input (typically denoted as Gaussian RBM or abbreviated as GRBM) is important. Alternatively, we can further generalize the SSD algoirthm that can be used to optimize GRBM.

Along the lines of [3], it is crucially required to derive an upper bound for loss function of GRBM in the sense of spectral norm or infinity norm. Since the GRBM will basically introduce another parameter $C^{-1}$ to represent the inverse of covariance matrix for observed variable and modify the energy function to quadratic accordingly, to bound the loss function is not as trivial for GRBM with real valued visible inputs. To achieve this desideratum, we make an additional but reasonable assumption: the norm of weight matrix is bounded by some constant. This assumption is in fact implicitly employed in real deep learning programming. Most of weights regularization tricks are forcing a smaller or bounded norm, such as renormalization or weights decay [13]. Under this assumption, we can successfully derive an upper bound of loss function by a local approximation with first order Taylor expansion plusing the error term in $l_{\infty}$ or Schattern-$\infty$ norm. To the end, we can unify the SSD algorithm for RBM with either binary or real valued inputs, being a substitute of vanilla SGD for RBM and able to apply various adaptive preconditioners.

1.0.1 Related Works

[3] is the most related work of this paper, since our work is its generalization, and a demonstration of its correctness. In addition, many stochastic gradient descent algorithms are developed for deep learning problems in recent years, including Adagrad [7], Rmsprop [25], AdaDelta [26], Adam [15]. These methods are all trying to approximate the inverse of Hessian matrix in an efficient way and use it as a preconditioner. However, the way they derive the update rule for parameter still falls into the Taylor approximation in the sense of $l_{2}$ or Frobenius norm, significantly differing from SSD. Another work [6] mentioned that the pathological curvature can make the algorithms get stuck more iterations for escaping the saddle point, and proposed another upper bound of loss function by adjusting quadratic term of Taylor expansion. It basically set the negative eigenvalue of Hessian matrix to be its absolute value and plugged it back to the quadratic term, which indicated it still considers the $l_{2}$ norm. AdaMax [15] is one algorithm considering the $l_{\infty}$ norm for vector space in a moving average way, which is also different from the S-$\infty$ in matrix space.

As to Gaussian RBM, a typical definition with diagonal covariance matrix is defined in [5], and this work focused on the parallel tempering algorithms. Our analysis of SSD w.r.t. GRM is not limited to diagonal covariance matrix, although it is more efficient and works well in real implementation. In [10]’s work, it introduced the model called Factorized Third-Order RBM intending to model pixel mean and covariances. This model can be considered as one variant of GRBM, while the covariance matrix is presumingly defined as some factorized structure. The parameterization of such a covariance matrix is complicated, and difficult to generalize with SSD in the theoretical aspect.

2 Preliminaries

2.1 RBM with Bernoulli Input

We briefly review the traditional RBM based on [22]. The RBM is a two-layer binary Markov Random Field, where the model includes observed layer units $\mathbf{v}\in\{0,1\}^{N_{v}}$ and its pairwise connected hidden units $\mathbf{h}\in\{0,1\}^{N_{h}}$. The “restricted” means there exists a matrix $W\in\mathbb{R}^{N_{v}\times N_{h}}$ to capture the symmetric interaction cross two layers but no self interaction within the layer, thus defining the energy function of the state $\{\mathbf{v},\mathbf{h}\}$ as

[TABLE]

where all parameters $\bm{\theta}=\{W,\mathbf{b},\mathbf{a}\}$. The probability density function of observed variable is

[TABLE]

where $Z(\bm{\theta})$ is the partition function. The objective is to minimize the negative log-likelihood over the dataset $\{\mathbf{v}_{n=1}^{N}\}$,

[TABLE]

Notice that we can decompose the objective function into the summation of two parts: the data independent log-partition function $f(\bm{\theta})=\log Z(\bm{\theta})$, and the unnormalized data negative log-likelihood $g(\bm{\theta})=-\frac{1}{N}\sum_{n=1}^{N}\sum_{\mathbf{h}}\exp\left(E(\mathbf{v}_{n},\mathbf{h};\bm{\theta})\right)$.

The derivative of the objective function w.r.t. $W$ is:

[TABLE]

where $\mathbb{E}_{p}[\cdot]$ represents the expectation w.r.t. the model defined distribution, and $\mathbb{E}_{e}[\cdot]$ denotes the empirical distribution defined on dataset. However, $\mathbb{E}_{p}[\cdot]$ involves expectation or summation that cannot be analytically computed. In practice parameter learning is performed by following an approximation to gradient of a different objective function, which means $\mathbb{E}_{p}[\cdot]$ is approximated by the distribution of samples obtained from running the Gibbs sampling of RBM for $k$ steps. This is called the Contrastive Divergence (CD-k) learning [9]. $k\rightarrow\infty$ recovers the $\mathbb{E}_{p}[\cdot]$, however, $k$ is set to 1 or small number in practice. Alternatively, Persistent CD [24] is more effective.

2.2 Schatten-$p$ Norm

For $p\geq 1$, the $l_{p}$ norm in finite dimensions for a vector $\mathbf{x}$ is defined as $\|\mathbf{x}\|_{p}=\left(\sum_{i}|x_{i}|^{p}\right)^{1/p}$. As $p\rightarrow\infty$, we obtain the $l_{\infty}$ norm, $\|\mathbf{x}\|_{\infty}=\max_{i}\{|x_{i}|\}$. Suppose $\bm{\lambda}(X)$ is the vector of singular values of matrix $X\in\mathbb{R}^{m\times n}$, we can similarly define the Schatten-$p$ norm $\|X\|_{S_{p}}=\|\bm{\lambda}(X)\|_{p}$, and $\|X\|_{S_{\infty}}=\|\bm{\lambda}(X)\|_{\infty}$. Notice singular values are non-negative, so we can omit the absolute operator for $\lambda_{i}$.

2.3 SSD for Bernoulli Inputs

Many stochastic gradient descent algorithms have been developed for discrete data modeling [3, 4, 2], due to the fact that

[TABLE]

where inner product $\langle\cdot,\cdot\rangle$ means trace of matrix product. In general, if we consider the error term above in $l_{2}$ norm, the local minimization of this objective function indicates $\Delta W\propto\nabla_{W}\mathcal{L}$, i.e., in practice the gradient estimated by CD-$k$ can be directly used for traditional gradient descent. However, in the sense of $S_{\infty}$ norm, the induced minimization is no longer to maximize in the direction of the gradient. It indicates $\Delta W\propto\|\bm{\lambda}(W)\|_{1}UV^{\top}$, where $U,V$ are unitary matrix obtained from proceeding SVD to $\nabla_{W}\mathcal{L}$. Analogously, SSD performs parameter update as $W\leftarrow W-\epsilon\Delta W$, where $\epsilon$ is the step size.

3 Gaussian RBM

In this section, we will generalize the SSD to handle Gaussian RBM (GRBM). Unlike Bernoulli RBM, the observed variable $\mathbf{v}\in\mathbb{R}^{N_{v}}$ is assumed to be a real-value vector as Gaussian distributed variables, while the binary latent variable $\mathbf{h}$ remains unchanged. The Gaussian RBM we discussed here is defined with the following energy function.

[TABLE]

where $\bm{\theta}=\{W,\mathbf{b},\mathbf{a},C^{-1}\}$. Notice that if we also set the hidden variable as Gaussian latent variable, it means that the resulted model would be an undirected version of factor analysis and however, it turns out to be identical to the traditional directed version [18]. If we only change the binary hidden variable of Bernoulli RBM to become following Gaussian distribution, it means that the induced model would be a special case of exponential PCA [19].

The parameters of energy function include an extra covariance matrix of observed variable, and if $C$ is set to be identity matrix, it will recover to the GRBM defined in [20]. Identity matrix is so simple that nothing to learn, thus diagonal covariance matrix [5] is usually preferred in practice due to its efficient computation. More complicated covariance matrix design (e.g., [10] parametrizes $C$ as a structured factorization) exists as well but beyonds our research area.

Since $\mathbf{v}$ is a real-valued vector, the summation in partition function $Z(\bm{\theta})$ (Eq. (3)) over $\mathbf{v}$ is replaced by integration, thus the log-partition function can be written as

[TABLE]

While the definitions of $p(\mathbf{v};\bm{\theta})$, $\mathcal{L}(\theta)$ and data negative likelihood function $g(\bm{\theta})$ will remain the same formula as Eq. (2, 4) in previous section, although their internal computation is completely different. According to energy function (6), the conditional distribution can be represented as,

[TABLE]

where $[\cdot]_{k}$ means the $k$-th element of vector. The tractable conditional distributions enable the possibility of Gibbs sampling for gradient computation (see supplementary materials for the detailed gradient derivation), and it means that CD-$k$ algorithm for gradients estimation can be applied to GRBM directly. However, we want to develop more efficient usage of gradients for parameter updates in the framework of SSD. It essentially requires that we can upper bound the loss function in the sense of $l_{\infty}$ or $S_{\infty}$ norm. Our method is to bound $f(\bm{\theta})$ and $g(\bm{\theta})$ separately.

3.1 Upper Bound of $f(\bm{\theta})$

In this section, we will bound the log-partition function w.r.t. each of $\bm{\theta}=\{W,\mathbf{b},\mathbf{a},C^{-1}\}$ when the other three are fixed or $\bm{\theta}^{old}$ is given. First, we handle the simplest cases for $\mathbf{a}$. By integrating out continuous variable $\mathbf{v}$, we can rewrite $f(\bm{\theta})$ as $f(\bm{\theta})=\log\left(\sum_{i=1}^{2^{N_{h}}}\omega_{i}\exp(\mathbf{h}_{i}^{\top}\mathbf{a})\right)$, where $\omega_{i}=(2\pi)^{\frac{N_{v}}{2}}|C|^{\frac{1}{2}}\exp\left(\mathbf{b}^{\top}C^{-1}W\mathbf{h}_{i}+\frac{(W\mathbf{h}_{i})^{\top}C^{-1}W\mathbf{h}_{i}}{2}\right)$. Notice the bianary vector $\mathbf{h}$ can have $2^{N_{h}}$ configurations. To bound $f(\bm{\theta})$ w.r.t. $\mathbf{a}$, we require the following Lemma and use it to obtain the upper bound shown in Theorem 10 (See proof in supplementary materials).

Lemma 1**.**

[3]** Consider log-sum-exp function $lse_{\bm{\omega}}(\mathbf{x})=\log\sum_{i}\omega_{i}\exp(x_{i})$ where $\bm{\omega}$ is independent of $\mathbf{x}$, the following inequality $lse_{\bm{\omega}}(\mathbf{x}+\Delta\mathbf{x})\leq lse_{\bm{\omega}}(\mathbf{x})+\langle\nabla lse_{\bm{\omega}}(\mathbf{x}),\Delta\mathbf{x}\rangle+\frac{1}{2}\|\Delta\mathbf{x}\|_{\infty}$ holds.

Theorem 2**.**

When fixing other parameters, we have

[TABLE]

In order to bound $f(\bm{\theta})$ of w.r.t. other parameters, we make an assumption that the weight matrix $W$ has a finite norm, i.e., $\|W\|_{2}\leq R$. This assumption is theoretically reasonable and practical in most deep learning implementation ($R$ is usually small, e.g. $R=\sqrt{15}$ in [13]), since either the weights decay or weights renormalization is to force the weights matrix norm to decrease. Thus we can readily bound $f(\bm{\theta})$ w.r.t. $\mathbf{b}$ by applying Lemma 9. We rewrite $f$ as $f(\bm{\theta})=\log\left(\sum_{i=1}^{2^{N_{h}}}\omega_{i}\exp(\mathbf{b}^{\top}C^{-1}W\mathbf{h}_{i})\right)$, where $\omega_{i}=(2\pi)^{\frac{N_{v}}{2}}|C|^{\frac{1}{2}}\exp\left(\mathbf{h}_{i}^{\top}\mathbf{a}+\frac{1}{2}(W\mathbf{h}_{i})^{\top}C^{-1}W\mathbf{h}_{i}\right)$

Theorem 3**.**

When fixing other parameters, we have

[TABLE]

where $r(C)$ is a local constant dependent on the parameter in previous iteration.

$r(C)$ actually depends on the eigenvalues of $C$, where $1/r(C)$ is between the minimum and maximum of the eigenvalues (See discussion in Supplementary materials). For inverse covariance matrix $C^{-1}$, we can derive the following representation by integrating out $\mathbf{v}$.

[TABLE]

where $\omega_{i}=(2\pi)^{\frac{N_{v}}{2}}\exp\left(\mathbf{h}_{i}^{\top}\mathbf{a}\right)$. However, this function involves log determinant, we need another lemma.

Proposition 4**.**

$-\log|C^{-1}|=\log|C|$.

Identity in Proposition 12 allows to consider $f(\bm{\theta})$ w.r.t. $C^{-1}$. Lamma 9 with assumption on $W$ indicates the upper bound for $f$.

Theorem 5**.**

Given parameters in previous iteration, $f(C^{-1}+U)\leq f(C^{-1})+\langle\nabla_{C^{-1}}f(C^{-1}),U\rangle+\left(\frac{N_{h}^{2}R^{2}}{2}+r(C)\right)\|U\|_{S_{\infty}}^{2}$, where $U$ is positive semidefinite, inner product is defined as trace for matrix, and $r(C)$ is a local constant.

For weight matrix $W$, we notice that the log-partition function has a canonical form that does not fall into $lse$ function, where $W$ appears as quadratic.

[TABLE]

where $\omega_{i}=(2\pi)^{\frac{N_{v}}{2}}|C|^{\frac{1}{2}}\exp\left(\mathbf{h}_{i}^{\top}\mathbf{a}-\frac{1}{2}\mathbf{b}^{\top}C^{-1}\mathbf{b}\right)$ and $s(W)=C^{-1/2}(W\mathbf{h}_{i}+\mathbf{b})$. Thus, we prove another lemma.

Lemma 6**.**

Consider log-sum-exp-square function $lse2_{\bm{\omega}}(\mathbf{x})=\log\sum_{i}\omega_{i}\exp(x_{i}^{2}/2)$ where $\bm{\omega}$ is independent of $\mathbf{x}$, if the domain of this function satisfies the bound condition $\|\mathbf{x}\|_{2}\leq r$, then the following inequality $lse2_{\bm{\omega}}(\mathbf{x}+\Delta\mathbf{x})\leq lse2_{\bm{\omega}}(\mathbf{x})+\langle\nabla lse2_{\bm{\omega}}(\mathbf{x}),\Delta\mathbf{x}\rangle+\left(\frac{1}{2}+\frac{3r^{2}}{4}\right)\|\Delta\mathbf{x}\|_{\infty}$ holds.

With the definition $lse2$, we can apply Lemma 14 to bound $f(\bm{\theta})$ w.r.t. $W$ in the sense of Schatten-$\infty$ norm.

Theorem 7**.**

When fixing other parameters, we have

[TABLE]

where $r(R,C,\mathbf{b})$ depends on the parameter $C,\mathbf{b}$ in previous iteration.

Notice that in practice we don’t care about each constant before all error terms in the sense of infinity norm. Since we always need to adjust the learning rate, the learning rate will absorb these constants and it is unnecessary to compute it.

3.2 Upper Bound of $g(\bm{\theta})$

In previous section, we derive the upper bound for log-partition function. We can also rewrite negative likelihood function $g(\bm{\theta})$ as

[TABLE]

However, we still need to derive the upper bound for $g$, if we want to bound the loss function $\mathcal{L}(\bm{\theta})$. Notice that due to the property of negative softplus function $-\log(1+e^{x})$, $g(\bm{\theta})$ is concave w.r.t. $W$ or $\mathbf{a}$, the inequality $g(\mathbf{a}+\Delta\mathbf{a})\leq g(\mathbf{a})+\langle\nabla_{\mathbf{a}}g(\mathbf{a}),\Delta\mathbf{a}\rangle$ naturally holds, so does for $W$. In fact, $g(\bm{\theta})$ is also concave w.r.t. $C^{-1}$, since the first part is linear to $C^{-1}$ (See details in supplementary materials).

As to $\mathbf{b}$, we notice $g$ is quadratic where the gradient that exceeds 2nd order is vainished, i.e., $g(\mathbf{b}+\Delta\mathbf{b})-g(\mathbf{b})=\langle\nabla_{\mathbf{b}}g(\mathbf{b}),\Delta\mathbf{b}\rangle+\frac{1}{2}\Delta\mathbf{b}^{\top}C^{-1}\Delta\mathbf{b}$. In addition, $\Delta\mathbf{b}^{\top}C^{-1}\Delta\mathbf{b}\leq N_{v}r(C)\|\Delta\mathbf{b}\|_{\infty}$, where $r(C)$ can be determined by the smallest eigenvalue of $C$ (Details see supplementary materials).

3.3 Revisiting Stochastic Spectral Descent

To this end, we have successfully bounded $f(\bm{\theta})$ and $g(\bm{\theta})$ w.r.t. each parameter in $\bm{\theta}$. Thus, the final loss function $\mathcal{L}(\bm{\theta})$ can also be bounded in the sense of infinity norm.

[TABLE]

The SSD update strategy is derived from minimizing the surrogate of loss function from the right-hand side of above inequality. Since the Schattern-$\infty$ norm mainly depends on the largest eigenvalue, the optimal direction to perform gradient descent is different from traditional stochastic gradient descent (SGD). However, if the error term is in $l_{2}$ norm, the update strategy will recover SGD update. Additionally, it indicates the inequality holds for both Bernoulli and Gaussian RBM, and guarantees that SSD can be applied. Thus, the general SSD algorithm is summarized in Algorithm 1.

Theorem 8**.**

Optimizing the local approximated loss function $\arg\min_{\bm{\theta}_{i}}\langle\nabla_{\bm{\theta}_{i}}\mathcal{L}(\bm{\theta}_{i}^{old}),\bm{\theta}_{i}-\bm{\theta}_{i}^{old}\rangle+c_{i}\|\bm{\theta}_{i}-\bm{\theta}_{i}^{old}\|_{\infty}^{2}$ will obtain the updates for parameters in Algorithm 2, 3.

3.3.1 Computation complexity

Another important issue w.r.t. SSD is its efficiency on the singular value decomposition (SVD). In most experiments, $C^{-1}$ is usually set as diagonal [20, 5], thus SVD operation is not necessary. For $N_{v}\times N_{h}$ weight matrix $W$, SVD requires $\mathcal{O}(N_{v}N_{h}\min\{N_{v},N_{h}\})$ algorithmic computation. However, this computation cost may be practically equivalent to SGD if $\min\{N_{v},N_{h}\}\approx Bk$, since the gradient computation by CD-$k$ algorithm requires $\mathcal{O}(N_{v}N_{h}Bk)$ computations, where $B$ is the batch size. Alternatively, we can use the randomized SVD algorithm [8] for speeding up with low rank approximation, reducing the computational complexity to $\mathcal{O}(NM+N_{v}N_{h}\log M)$ where $N=\max\{N_{v},N_{h}\}$ and $M$ is approximated level. Meanwhile, the update rule will be modified as $\Delta X=\|\bm{\lambda}\|_{1}UV^{\top}+\|\bm{\lambda}\|_{1}(X-U\bm{\lambda}V^{\top})/\|X-U\bm{\lambda}V^{\top}\|_{S_{\infty}}$. This approach is fast in practice and will not significantly harm the performance.

3.3.2 Discussion

Additionally, if $\mathcal{L}$ is Lipschitz continuous, the convergence rate is $\mathcal{O}\left(L_{p}\|\bm{\theta}-\bm{\theta}^{*}\|_{p}^{2}/T\right)$, where $\bm{\theta}^{*}$ is local optimal and $L_{p}$ is Lipschitz smoothness constant w.r.t. $p$ norm. This implies that SSD does not substantially increase the computational cost but may potentially encourage faster convergence than SGD (correspond to $L_{2}\|\bm{\theta}-\bm{\theta}^{*}\|_{2}^{2}$) if $L_{\infty}\|\bm{\theta}-\bm{\theta}^{*}\|_{\infty}^{2}$ is smaller. [14] proved the conditions that optimization method favors infinity norm, however, it is difficult to verify for RBM with such a complicated non-convex loss function. This is beyond our research and we will leave it as future work.

Unlike the trend of SGD based algorithms, by applying Adagrad [7], Rmsprop [25], AdaDelta [26], Adam [15] on top of SGD, which are universally suitable for most deep learning problems, SSD is specifically designed for RBM. Additionally, SSD does not focus on how to adapt the learning rate or step-size by leveraging the gradients of previous iterations, but is basically trying to adjust the descent direction by taking advantage of geometry information of current gradient. However, for Gaussian RBM we notice that the covariance matrix $C^{-1/2}$ always plays the role of preconditioner in gradients derivation. This means we can naturally adjust or replace $C^{-1/2}$ with other preconditioner; more efficiently, we can directly apply most of the preconditioners, like Adagrad, Rmsprop, and Adam, on top of SSD. In summary, SSD is parallel to SGD.

4 Experiments

Since SSD is an alternative to SGD for RBM, we mainly compare the performance in vanilla version. As mentioned in discussion, SSD with adaptive step-size can also be implemented readily. However, this may misinterpret the significance of spectral descent, since adaptive step-size methods can also accelerate the practical convergence in general.

4.1 Simulation Study

To quantitatively observe the behaviors or properties of SSD algorithm, our first experiment was run on synthetic dataset, with $N_{v}=100$ and $N_{h}=25$. 4000 training and 1000 held-out testing data points are generated as binary vectors to test both Binary RBM and Gaussian RBM. The weights matrix $W$ was a random Gaussian matrix with a variance 0.5 for each element, and the biases are set as zero for convenience. However, we pretended this is an unknown information and allowed the algorithm to learn biases. In our model setting, we only set the $N_{v}$ and $N_{h}$ as simulated ones.

The reconstruction error results are shown in Figure 2. Figure 2(a) illustrates the training performance evaluated by Binary RBM. To closely see the significance of SSD, we also implemented the hybrid descent algorithms, e.g., SSDW/SGD represents only the weight matrix $W$ is updated according to spectral descent while all the other parameters follow the SGD update. We gradually modified updating all parameters by SGD to become updating by SSD. From Figure 2(a), we can see that updating $W$ is crucial for training. Updating $W$ alone by SGD (i.e., SGDW/SSD) has the similar behavior as updating all parameters by SGD, however, once we changed the updating mode to SSDW/SGD, the behavior will become more similar to pure SSD, no matter for small step-size or large step-size. In addition, even we use the SSD to update all parameters, the performance is just slightly better than SSDW/SGD. For Gaussian RBM, we use a simple model with isotropic covariance matrix $c\mathbf{I}$. Since we use the binary input, the GRBM should model worse than BRBM, as proved in Figure 2(b). However, we can still observe that SSD converges faster than SGD, and this phenomenon can also be clearly visualized in Figure 3, where SSD prefers to find a shortcut to the local optimal.

4.2 MNIST and FreyFace Datasets

In this experiment, we use the vectorized MNIST digit dataset that includes 60,000 training and 10,000 test images of handwritten digits with size $28\times 28$. We test both Gaussian and Binary RBM for this dataset, though it is commonly modeled as BRBM. For this real dataset, we assume the covariance matrix as diagonal, i.e. $\text{diag}(\mathbf{c})$, where $\mathbf{c}$ is a vector with positive elements. We mainly evaluate sum of squared error for pixel-wise reconstruction, since the true normalized log-likelihood is intractable for RBM, which usually requires the annealing importance sampling [22, 16] to obtain an approximate result. Such an approximation may suffer high variance if insufficient samples are used. Additionally, we also evaluate the running time of SSD compared with SGD, when CD-10 is applied on a mini-batch with size 100. Since we found vanilla SGD converges too slow (even with non-fixed but exponentially decay stepsize), we applied the Nesterov acceleration [23] for SGD, but did not applied it to SSD, which seems an unfair comparison.

Figure 2(d) shows the results running on MNIST. We found SSD on GRBM can achieve faster convergence rate than BRBM, though the error of initial hundreds of iterations is relatively large. Unlike the simulation, this is a counter-intuitive phenomenon, since our images has been binarized. However, one possible reason is that the covariance matrix enables more flexible modeling, allowing the more variation for same pixel position. Even we use the Nesterov SGD, there still exists a gap between the SSD and SGD for the first 30,000 iterations, where the performance of SGD has become stable but SSD implicates a trend for further decline. Figure 2(e) shows the different weight matrix learned by SGD and SSD. The SGD learned matrix is more random while the SSD can learn more sparsity or patterns. If we plot the histogram for all the elements of $W$, we can also observe the similar sparsity for SSD (See figures in supplements).

We also compared the running time between SGD and SSD. For fair comparison, we also include the vanilla version of SGD. The results of wall-clock time illustrated in Table 1 make sense since SVD and momentum update will require extra computations.

The FreyFace111http://cs.nyu.edu/~roweis/data/frey_rawface.mat dataset contains almost 2000 images of some person’s face, taken from sequential frames of a small video with size $20\times 28$. Each pixel value of face images in this dataset is real value between 0 and 1, which should fit GRBM well. We can also consider it as the Bernoulli probability for observation and apply BRBM, thus we test both BRBM and GRBM and show the result in Figure 4. Similar to MNIST experiment, Nesterov acceleration is implemented for SGD, however, SGD on GRBM converges too slow even with thoroughly fine-tuning. Thus, we did not include this result in the figure. A counter-intuitive phenomenon also appears for SSD, i.e., BRBM can model the data better than GRBM, even if the data point is real-valued. One possible explanation is that the size of this dataset is small such that the more flexible GRBM is not well trained.

5 Conlusion and Future Work

In this paper, we generalize the stochastic spectral descent algorithm to a broader RBM family, enlarging the domain of SSD. A unifying SSD algorithm is summarized due to the newly proved upper bound. Since SSD empirically shows potential advantages over stochastic gradient descent, the emergence of SSD may imply an alternative or substitute of SGD for RBM related model, or adaptive step-size modification. An important avenue of improvement might explore other gradient estimation method for RBM and integrate with SSD. For example, Hamiltonian Monte Carlo (HMC) method [21] allows efficient sampling from conditional distributions that might be difficult to sample from, like mean-covariance RBM [10, 17]. HMC may also be helpful to improve the overall performance of SSD, and the interaction between HMC and SSD is our future work. Additionally, beyond the discussion of traditional $l_{2}$ norm and infinite norm, $l_{0}$ or $l_{1}$ norm constrain can be a new direction to derive other optimization algorithms which allow partial updates for its sparse property.

6 Appendices

6.1 Gradient Computation

[TABLE]

Consider special case, $C^{-1}$ is diagonal, thus denoting as a vector.

[TABLE]

Furthermore, the element of diagonal should be positive, $C=\text{diag}\{e^{c_{1}},\dots,e^{c_{R}}\}$.

[TABLE]

if $C^{-1}=\text{diag}\{e^{c_{1}},\dots,e^{c_{R}}\}$,

[TABLE]

6.2 Derivation of $g(\bm{\theta})$

$N$ is the size of dataset.

[TABLE]

The first part of $g$ is linear to $C^{-1}$, and the second part is concave w.r.t. $W$, $\mathbf{a}$, and $C^{-1}$, since $-\log(1+e^{x})$ is concave.

[TABLE]

One option for $r(C)$ is $1/\lambda_{\min}(C)$, where $\lambda_{\min}(C)$ is minimum eigenvalue of $C$.

6.3 Theorem Proof

Lemma 9**.**

Consider log-sum-exp function $lse_{\bm{\omega}}(\mathbf{x})=\log\sum_{i}\omega_{i}\exp(x_{i})$ where $\bm{\omega}$ is independent of $\mathbf{x}$, the following inequality $lse_{\bm{\omega}}(\mathbf{x}+\Delta\mathbf{x})\leq lse_{\bm{\omega}}(\mathbf{x})+\langle\nabla lse_{\bm{\omega}}(\mathbf{x}),\Delta\mathbf{x}\rangle+\frac{1}{2}\|\Delta\mathbf{x}\|_{\infty}$ holds.

Proof.

See [3]. ∎

Theorem 10**.**

When fixing other parameters, we have

[TABLE]

Proof.

The log partition function can be written as a sum over only the hidden units to give a similar form to Lemma 9.

[TABLE]

where

[TABLE]

and it does not depend on $\mathbf{a}$. Define $\mathbf{H}\in\{0,1\}^{N_{h},2^{N_{h}}}$ be the matrix with each column taking one possible configuration of $\mathbf{h}$, thus we have

[TABLE]

By Lemma 9, we have

[TABLE]

where we omit other fixed parameters in $\bm{\theta}$. Notice the fact that $\mathbf{H}\nabla_{\mathbf{H}^{\top}\mathbf{a}}lse(\mathbf{H}^{\top}\mathbf{a})=\nabla_{\mathbf{a}}f(\mathbf{a})$, thus the $(\nabla_{\mathbf{H}^{\top}\mathbf{a}}lse(\mathbf{H}^{\top}\mathbf{a}))^{\top}\mathbf{H}^{\top}\Delta\mathbf{a}=(\nabla_{\mathbf{a}}f(\mathbf{a}))^{\top}\Delta\mathbf{a}$ holds.

In addition, $\|\mathbf{H}^{\top}\Delta\mathbf{a}\|_{\infty}^{2}=\max_{i}|\mathbf{h}_{i}^{\top}\Delta\mathbf{a}|\leq N_{h}\|\Delta\mathbf{a}\|_{\infty}^{2}$. ∎

Theorem 11**.**

When fixing other parameters, we have

[TABLE]

where $r(C)$ is a local constant dependent on the parameter in previous iteration.

Proof.

We can also rewrite log partition function as the following formulation.

[TABLE]

where $\omega_{i}=(2\pi)^{\frac{N_{v}}{2}}|C|^{\frac{1}{2}}\exp\left(\mathbf{h}_{i}^{\top}\mathbf{a}+\frac{1}{2}(W\mathbf{h}_{i})^{\top}C^{-1}W\mathbf{h}_{i}\right)$. Following the derivation of proof in Theorem 10, we only need to consider $\max_{i}|\Delta\mathbf{b}^{\top}C^{-1}W\mathbf{h}_{i}|$. Using the assumption $\|W\|_{2}\leq R$, this term is bounded by $N_{h}R\|\Delta\mathbf{b}^{\top}C^{-1}\|_{\infty}^{2}$. Denote $\lambda_{\min}(C)$ and $\lambda_{\max}(C)$ as the maximum and minimum eigen values of covariance matrix $C$, thus there exist a constant $r(C)\in[\lambda_{\min}(C),\lambda_{\max}(C)]$ such that the desired term bounded by $N_{h}Rr(C)\|\Delta\mathbf{b}\|_{\infty}^{2}$. Obviously, $r(C)=1/\lambda_{\min(C)}$ is one option. ∎

Proposition 12**.**

$-\log|C^{-1}|=\log|C|$.

Theorem 13**.**

When fixing other parameters, we have

[TABLE]

where $U$ is positive semidefinite and inner product is defined as trace for matrix.

Proof.

With Proposition 12 we can rewritten log partition function as

[TABLE]

where $\omega_{i}=(2\pi)^{\frac{N_{v}}{2}}\exp\left(\mathbf{h}_{i}^{\top}\mathbf{a}\right)$. We also have

[TABLE]

wher $A$ is non-negative and there exists a constant $r^{\prime}(C)$ such that $A\leq r^{\prime}(C)\|vec(U)\|_{2}^{2}=r^{\prime}(C)\|U\|_{F}^{2}=r^{\prime}(C)\|\bm{\lambda}(U)\|_{2}^{2}$ where $\bm{\lambda}(U)$ is the singular value vector, thus there exists a constant $r(C)=N_{v}^{2}r^{\prime}(C)$ such that $A\leq r(C)\|U\|_{S_{\infty}}^{2}$. ∎

Lemma 14**.**

Consider log-sum-exp-square function $lse2_{\bm{\omega}}(\mathbf{x})=\log\sum_{i}\omega_{i}\exp(x_{i}^{2}/2)$ where $\bm{\omega}$ is independent of $\mathbf{x}$, if the domain of this function satisfies the bound condition $\|\mathbf{x}\|_{2}\leq r$, then the following inequality $lse2_{\bm{\omega}}(\mathbf{x}+\Delta\mathbf{x})\leq lse2_{\bm{\omega}}(\mathbf{x})+\langle\nabla lse2_{\bm{\omega}}(\mathbf{x}),\Delta\mathbf{x}\rangle+\left(\frac{1}{2}+\frac{3r^{2}}{4}\right)\|\Delta\mathbf{x}\|_{\infty}$ holds.

Proof.

Denote $lse2_{\bm{\omega}}(\mathbf{x})$ as $l(\mathbf{x})$ for simplicity.

[TABLE]

where $s(\mathbf{x})=\exp(l(\mathbf{x}))$. Let $\pi_{i}=\frac{\omega_{i}\exp(x_{i}^{2}/2)}{s(\mathbf{x})}$, then $\sum_{i}\pi_{i}=1$, $(1+x_{i}^{2})\pi_{i}-x_{i}^{2}\pi_{i}^{2}\leq\pi_{i}+\frac{x_{i}^{2}}{4}$. Then for $\mathbf{u}$ where $|u_{i}|<|x_{i}|,\forall i$,

[TABLE]

where $\mathbf{H}_{l(\mathbf{x}+t\mathbf{u})}$ is the hessian matrix of $f$ evaluated at $\mathbf{x}+t\mathbf{u}$. The integral remainder can be locally bounded, since the following statement.

[TABLE]

Or we can write

[TABLE]

Thus we conclude the Lemma. ∎

Theorem 15**.**

When fixing other parameters, we have

[TABLE]

where $r(R,C,\mathbf{b})$ depends on the parameter $C,\mathbf{b}$ in previous iteration.

Proof.

Rewrite log partition function as quadratic form.

[TABLE]

where $\omega_{i}=(2\pi)^{\frac{N_{v}}{2}}|C|^{\frac{1}{2}}\exp\left(\mathbf{h}_{i}^{\top}\mathbf{a}-\frac{1}{2}\mathbf{b}^{\top}C^{-1}\mathbf{b}\right)$ and $s(W)=C^{-1/2}(W\mathbf{h}_{i}+\mathbf{b})$. Follow Lemma 14 and previous proof can conclude the result. Notice $r(R,C,\mathbf{b})$ also depends the eigenvalues of $C$. ∎

Theorem 16**.**

Optimizing the local approximated loss function $\arg\min_{\bm{\theta}_{i}}\langle\nabla_{\bm{\theta}_{i}}\mathcal{L}(\bm{\theta}_{i}^{old}),\bm{\theta}_{i}-\bm{\theta}_{i}^{old}\rangle+c_{i}\|\bm{\theta}_{i}-\bm{\theta}_{i}^{old}\|_{\infty}^{2}$ will obtain the updates for parameters in Algorithm 2, 3.

Proof.

If $\bm{\theta}_{i}$ is a vector, we have

[TABLE]

This will conclude the update in Algorithm 3.

If $\bm{\theta}_{i}$ is a matrix, we first apply SVD on $\nabla_{\bm{\theta}_{i}}\mathcal{L}(\bm{\theta}_{i}^{old})=U\Sigma V^{\top}$. Thus we can decompose $\Delta\bm{\theta}_{i}=\bm{\theta}_{i}-\bm{\theta}_{i}^{old}$ as $U\Lambda V^{\top}$.

We further have

[TABLE]

This will conclude the update in Algorithm 2. ∎